处理odeint中极小/极大的值

Question

我正在尝试使用

scipy.integrate.odeint

求解这组耦合微分方程。

我写了这个Python代码：

def odes(x, r):
    #constants
    alpha = 1.473
    beta = 52.46
    gamma = 1.33 
    
    #Assign each ode to a vector element
    p = x[0]
    m = x[1]
    
    #Define each ode
    dp_dt = -alpha*(x[0]**(1/gamma))*x[1]/(r**2)
    dm_dt = beta*(r**2)*(x[0]**(1/gamma))
    
    #Return list
    return [dp_dt, dm_dt]

#initial conditions(Initial p_bar, initial m_bar)
init = [1e-16, 0]

#declare a r vector
r = np.linspace(0, 15000, 1000000)

#Solve the odes
x = odeint(odes, init, r)

print(x)
plt.semilogy(r, x[:,0])
plt.semilogy(r, x[:,1])
plt.show()

然而，Jupyter 抛出了这个

RuntimeWarning

：

C:\conda_temp\ipykernel_14612\3000110852.py:12: RuntimeWarning: invalid value encountered in double_scalars
  dp_dt = -alpha*(x[0]**(1/gamma))*x[1]/(r**2)

这就是

x

数组的样子：

[[1.e-16 0.e+00]
 [   nan    nan]
 [0.e+00 0.e+00]
 ...
 [0.e+00 0.e+00]
 [0.e+00 0.e+00]
 [0.e+00 0.e+00]]

根据我收集的信息，这与 Python/

odeint

处理极小值的方式有关。但是，我不知道如何解决这个问题。我尝试过使用

scipy.special.logsumexp

（错误？）但没有成功。

Answer 1

您的问题不是数字的大小，您的问题是您的第一个

r

值为 0，这意味着您在计算

dp_dt

时除以零。相反，将第一个

r

值设置为较小但非零的值，例如

r = np.linspace(1e-6, 15000, 1000000)