An error in odeint program

Question

I wrote a program to use odeint to solve a differential equation. But it had a problem. When I setted Cosmopara as np.array([70.0,0.3,0,-1.0,0]), it gave a warning that invalid value encountered in sqrt and invalid value encountered in double_scalars in h = np.sqrt(y1 ** 2 + Omega_M * t ** (-3) + Omega_DE*y2). But I checked that line and didn't find any error. If Cosmopara = np.array([70.0,0.3,0.0,-1.0,0.0]), Y shouldn't change but it changed. Besides, If I chose Cosmopara = np.array([70.0,0.3,0.1,-1.0,0.1]), this program could gives the right result.

What am I doing wrong?

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

global CosmoPara
Cosmopara = np.array([70.0,0.3,0.0,-1.0,0.0])

def derivfun(Y,t):
    Omega_M = Cosmopara[1]
    Sigma_0 = Cosmopara[2]
    omega = Cosmopara[3]
    delta = Cosmopara[4]
    Omega_DE = 1-Omega_M-Sigma_0**2
    y1 = Y[0]
    y2 = Y[1]
    h = np.sqrt(y1**2 + Omega_M*t**(-3) + Omega_DE*y2)
    dy1dt = -3.0*y1/t + (delta*Omega_DE*y2)/(t*h)
    dy2dt = -(3.0*(1+omega) + 2.0*delta*y1/h)*y2/t
    return np.array([dy1dt,dy2dt])

z = np.linspace(1,2.5,15001)
time = 1.0/z
Omega_M = Cosmopara[1]
Sigma_0 = Cosmopara[2]
omega = Cosmopara[3]
delta = Cosmopara[4]
Omega_DE = 1-Omega_M-Sigma_0**2
y1init = Sigma_0
y2init = 1.0
Yinit = np.array([y1init,y2init])

Y = odeint(derivfun,Yinit,time)

y1 = Y[:,0]
y2 = Y[:,1]

h = np.sqrt(y1**2 + Omega_M*time**(-3) + Omega_DE*y2)

plt.figure()
plt.plot(z,h)
plt.show()

Answer 1

I had a similar problem, it seems to arrise because as previousely stated, the odeint program seems to want to take a shortcut by increasing the step size. I fixed this for myself by capping the step size. You can do this very easilly:

y = odeint(model, x, t, hmax= maximum step size)

This does make the computing time quite a bit slower, so for big computations and long time ranges, you want to avoid it if possible. I found it also worked to decrease the time range, but this is not always possible. (in my case, I wanted to do a time range of 2000-20000 seconds, but I couldn't go any higher than 200 before values would just be all over the place)

Answer 2

The error caused by the situation when the value in the sqrt() is less than 0. And the value in the sqrt() will be less than 0 is caused by the time value t suddenly decreasing to -0.0000999.

Reason

I have tested several cases to find out the reason, and I found that the time spacing of the time value t given to the function derivfun and the time spacing of the global array time are different even when the function odeint works fine. I guess it is a mechanism for speeding up the odeint. Because when the derivatives dydt1 and dydt2 do not change fast, the result can be considered as a linear function with a larger time spacing. If the derivatives will not change fast, this function will increase next step's time spacing and the function can solve it in less steps. In the case you provided. the derivatives dydt1 and dydt2 always equal to zero and do not change. This situation causes the error time value.

Based on this result, we can know that the function odeint will give the wrong time value which is not in the range of the time array users give when the derivatives do not change. If you want to avoid this situation, you can use global time variables instead of the original time value. But it will cost more time when you using odeint to solve ordinary differential equations.

Code

The following code is the code you provided and I add some lines for testing. You can change the parameters and check out the result.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

Cosmopara = np.array([70.0,0.3,0.1,-1.0,0.1])
i = 0
def derivfun(Y,t):
    global i
    Omega_M = Cosmopara[1]
    Sigma_0 = Cosmopara[2]
    omega = Cosmopara[3]
    delta = Cosmopara[4]
    Omega_DE = 1-Omega_M-Sigma_0**2
    y1 = Y[0]
    y2 = Y[1]
    h = np.sqrt(y1**2 + Omega_M*t**(-3) + Omega_DE*y2)
    dy1dt = -3.0*y1/t + (delta*Omega_DE*y2)/(t*h)
    dy2dt = -(3.0*(1+omega) + 2.0*delta*y1/h)*y2/t
    print "Time: %14.12f / Global time: %14.12f / dy1dt: %8.5f / dy2dt: %8.5f" % (t, time[i], dy1dt, dy2dt)
    i += 1
    return np.array([dy1dt,dy2dt])

z = np.linspace(1,2.5,15001)
time = 1.0/z
Omega_M = Cosmopara[1]
Sigma_0 = Cosmopara[2]
omega = Cosmopara[3]
delta = Cosmopara[4]
Omega_DE = 1-Omega_M-Sigma_0**2
y1init = Sigma_0
y2init = 1.0
Yinit = np.array([y1init,y2init])

Y = odeint(derivfun,Yinit,time)

y1 = Y[:,0]
y2 = Y[:,1]

h = np.sqrt(y1**2 + Omega_M*time**(-3) + Omega_DE*y2)

plt.figure()
plt.plot(z,h)
plt.show()

Testing results

The following picture shows that the time value in the function and the global time value have different time spacing when the function works fine (using the parameters you provided):

The following picture shows that when the derivatives dydt1 and dydt2 do not change (using the parameters you provided, too), and the time value in the function suddenly decreases to a value less than 0: