odeint程序中的一个错误

Question

我写了一个程序，用odeint求解微分方程。但有个问题。当我将Cosmopara设置为np.array([70.0,0.3,0,-1.0,0])时，它给出了一个警告，即invalid value encountered in sqrt和{}位于{}。但我检查了那条线，没有发现任何错误。如果Cosmopara=np.array([70.0,0.3,0.0,-1.0,0.0])，Y不应该改变，但是它改变了。另外，如果我选择Cosmopara=np.array([70.0,0.3,0.1,-1.0,0.1])，这个程序可以给出正确的结果。在

我做错什么了？在

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

当sqrt()中的值小于0时导致的错误。并且sqrt()中的值将小于0是由于时间值t突然减小到{}。在

Reason

为了找出原因，我测试了几个例子，我发现函数derivfun的时间值t的时间间隔和全局数组{}的时间间隔是不同的，即使函数odeint工作正常。我想这是一种加速odeint的机制。因为当导数dydt1和{}变化不快时，结果可以看作是具有较大时间间隔的线性函数。在导数变化不快的情况下，该函数将增大下一步的时间间隔，并且可以用较少的步骤求解。在你提供的情况下。导数dydt1和{}始终等于零且不改变。这种情况会导致错误的时间值。在

根据这个结果，我们可以知道当导数不变时，odeint函数会给出错误的时间值，而这不在用户给出的时间数组范围内。如果要避免这种情况，可以使用全局时间变量而不是原始时间值。但是用odeint来求解常微分方程会花费更多的时间。在

Code

下面的代码是您提供的代码，我添加了一些行进行测试。您可以更改参数并检查结果。在

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

下图显示，当函数正常工作时（使用您提供的参数），函数中的时间值和全局时间值具有不同的时间间隔：

下图显示，当导数dydt1和dydt2不变（也使用您提供的参数），并且函数中的时间值突然减小到小于0的值：