python - An error in odeint program -

i wrote program use odeint solve differential equation. had problem. when setted cosmopara np.array([70.0,0.3,0,-1.0,0]), gave warning invalid value encountered in sqrt , invalid value encountered in double_scalars in h = np.sqrt(y1 ** 2 + omega_m * t ** (-3) + omega_de*y2). checked line , didn't find error. if cosmopara = np.array([70.0,0.3,0.0,-1.0,0.0]), y shouldn't change changed. besides, if chose cosmopara = np.array([70.0,0.3,0.1,-1.0,0.1]), program gives right result.

what doing wrong?

import numpy np import matplotlib.pyplot plt 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() 

the error caused situation when value in sqrt() less 0. , value in sqrt() less 0 caused time value t decreasing -0.0000999.

reason

i have tested several cases find out reason, , found time spacing of time value t given function derivfun , time spacing of global array time different when function odeint works fine. guess mechanism speeding odeint. because when derivatives dydt1 , dydt2 not change fast, result can considered linear function larger time spacing. if derivatives not change fast, function increase next step's time spacing , function can solve in less steps. in case provided. derivatives dydt1 , dydt2 equal 0 , not change. situation causes error time value.

based on result, can know function odeint give wrong time value not in range of time array users give when derivatives not change. if want avoid situation, can use global time variables instead of original time value. cost more time when using odeint solve ordinary differential equations.

code

the following code code provided , add lines testing. can change parameters , check out result.

import numpy np import matplotlib.pyplot plt scipy.integrate import odeint  cosmopara = np.array([70.0,0.3,0.1,-1.0,0.1]) = 0 def derivfun(y,t):     global     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)     += 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 time value in function , global time value have different time spacing when function works fine (using parameters provided):

the following picture shows when derivatives dydt1 , dydt2 not change (using parameters provided, too), , time value in function decreases value less 0: