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

Question

我写了以下代码，在打印了大约10个输出后，它会出错

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

前三个函数是正确的，在其他程序中工作良好。请不要对它们的结构敏感，我认为问题在于for loop和选择一些不能满足sqrt条件的值。我可以添加一些行来告诉Python避免使用导致数学域错误的数字吗？如果是，我应该怎么做？我的意思是通过sqrt下导致负值的步骤？

Cov是一个31*31mtrix，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

Answer 1

您的代码正在尝试查找负数的平方根。

这对于实值数字是不允许的（出于明显的数学原因），所以你最好的选择是使用try/除块：

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

Answer 2

以下是一些选项：

>

将sqrt（...）替换为（...）**0.5。这将产生复数，这可能是可接受的，也可能是不可接受的。例如（-1）**0.5生成i，在Python中显示为1j（忽略浮点错误）。

抓住错误，继续前进。由于您可能希望捕获更高级别的错误，我建议将ValueError从sqrt转换为自定义错误：

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]