当我们有 sqrt 时如何避免数学域错误

How to avoid math domain error when we have sqrt

我写了下面的代码,在打印大约 10 个输出后出现错误

return 1/sqrt(Om*(1+z)**3+omg0*(1+z)**6+(1-Om-omg0))
ValueError: math domain error

前三个函数是正确的,在其他程序中也能正常工作请不要对它们的结构敏感,我认为问题出在for-loop和选择一些不能满足sqrt条件的值。我可以添加一些行来告诉 Python 避免导致 math domain error 的数字吗?如果是,我应该怎么做?我的意思是在 sqrt?

下通过导致负值的步骤

Cov 是一个 31*31 矩阵,xx[n] 是一个包含 31 个数字的列表。

def ant(z,Om,w):
    return 1/sqrt(Om*(1+z)**3+w*(1+z)**6+(1-Om-w))

def dl(n,Om,w,M):  
    q=quad(ant,0,xx[n],args=(Om,w))[0]
    h=5*log10((1+xx[n])*q)
    fn=(yy[n]-M-h)                  
    return fn

def c(Om,w,M):
    f_list = []
    for i in range(31):  # the value '2' reflects matrix size
        f_list.append(dl(i,Om,w,M))
    A=[f_list]
    B=[[f] for f in f_list]
    C=np.dot(A,Cov)
    D=np.dot(C,B)
    F=np.linalg.det(D)*0.000001
    return F

N=100
for i in range (1,N):
    R3=np.random.uniform(0,1)

    Omn[i]=Omo[i-1]+0.05*np.random.normal()
    wn[i]=wo[i-1]+0.05*np.random.normal()
    Mn[i]=Mo[i-1]+0.1*np.random.normal()

    L=exp(-0.5*(c(Omn[i],wn[i],Mn[i])-c(Omo[i-1],wo[i-1],Mo[i-1])))

    if L>R3:
        wo[i]=wn[i]

    else:
        wo[i]=wo[i-1]

    print(wo[i])

输出为:

0.12059556415714912
0.16292726528216397
0.16644447885609648
0.1067588804671105
0.0321446951572128
0.0321446951572128
0.013169965429457382
Traceback (most recent call last):
......
return 1/sqrt(Om*(1+z)**3+omg0*(1+z)**6+(1-Om-omg0))
ValueError: math domain error

您的代码试图计算负数的平方根。

这对于实数值是不允许的(出于明显的数学原因),所以最好的办法是使用 try/except 块:

def ant(z,Om,w):
    try:
        return 1/sqrt(Om*(1+z)**3+omg0*(1+z)**6+(1-Om-omg0))
    except ValueError:
        return None

这里有一些选项:

  1. sqrt(...) 替换为 (...)**0.5。这将产生复杂的数字,这可能会或可能不会被接受。例如 (-1)**0.5 生成 i,在 Python 中将显示为 1j(忽略浮点错误)。

  2. 发现错误并继续。由于您可能想捕获更高级别的错误,我建议将 ValueErrorsqrt 转换为自定义错误:

class SqrtError(ValueError):
    pass

def ant(z, Om, w):
    val = Om * (1 + z) ** 3 + w * (1 + z) ** 6 + (1 - Om - w)
    try:
        sq = sqrt(val)
    except ValueError:
        raise SqrtError
    return 1 / sq

那么听起来您想继续尝试新的随机数,直到找到一个有效的随机数,这可以像这样完成:

for i in range(1, N):
    R3 = np.random.uniform(0, 1)
    while True:

        Omn[i] = Omo[i - 1] + 0.05 * np.random.normal()
        wn[i] = wo[i - 1] + 0.05 * np.random.normal()
        Mn[i] = Mo[i - 1] + 0.1 * np.random.normal()

        try:
            L = exp(-0.5 * (c(Omn[i], wn[i], Mn[i]) - c(Omo[i - 1], wo[i - 1], Mo[i - 1])))
        except SqrtError:
            continue
        else:
            break

    if L > R3:
        wo[i] = wn[i]

    else:
        wo[i] = wo[i - 1]