我试图用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（不正确？）无济于事。