How to avoid math domain error when we have sqrt

Question

I wrote the following code, after printing about 10 outputs it gets error

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

the three first functions are correct and work well in other programs please do not be sensitive about their structures, I think the problem is in for-loop and choosing some values which cannot satisfy the condition of sqrt. Could I add some lines to tell Python avoid numbers which lead to math domain error? If yes, How should I do that? I mean passing the step which leads to negative value under sqrt?

Cov is a 31*31 mtrix, xx[n] is a list of 31 numbers.

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])

The output is:

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

Here are some options:

1. Replace sqrt(...) with (...)**0.5. This will produce complex numbers, which may or may not be acceptable. For example (-1)**0.5 produces i which in Python will appear as 1j (ignoring floating point errors).

2. Catch errors and move on. Since you probably want to catch the errors higher up, I recommend translating the ValueError from sqrt into a custom error:

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

Then it sounds like you want to keep trying new random numbers until you get one that works, which could be done like so:

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]