计算科学 - Scipy odeint 意外结果 - 吾爱随笔录

我正在尝试对方程进行数值积分

\frac{d P}{d r} = - (P + ρ (r)) \frac{m (r) + 4 π r^{3} P}{r [r - 2 m (r)]}

$\frac{\mathrm{dP} }{\mathrm{d} r}=-\left ( P+\rho\left ( r \right ) \right )\frac{m\left ( r \right )+4\pi r^{3}P}{r\left [ r-2m\left ( r \right ) \right ]}$ 去寻找

P (r)

$P\left ( r \right )$ . 我已经定义了以下内容：

ρ (r) = ρ_{c} [1 - {(\frac{r}{R})}^{2}]

$\rho\left ( r \right )=\rho_{c}\left [ 1-\left ( \frac{r}{R} \right )^{2} \right ]$

m (r) = \frac{4 π ρ_{c} (5 r^{3} R^{2} - 3 r^{5})}{15 R^{2}} .

$m\left ( r \right )=\frac{4\pi \rho_{c}\left ( 5r^{3}R^{2}-3r^{5} \right )}{15R^{2}}.$ 我的代码如下：

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

def dP_dr(x, r):
    return -(x+density(r, p_c, R))*(mass(r, p_c, R)+4*np.pi*r**3*x)/(r*(r- 
        2*mass(r, p_c, R)))

def mass(r, p_c, R):
    return 4*np.pi*p_c*(5*r**3*R**2-3*r**5)/(15*R**2)

def density(r, p_c, R):
    return p_c*(1-(r/R)**2)

p_c = 7*10**17      #central density
P_c = 1.6*10**34    #central pressure
R = 10              #radius

deltaR = 0.01
r = np.linspace(deltaR, R, int(R/deltaR))

P = odeint(dP_dr, P_c, r)

plt.figure()
plt.plot(r, P, label='Pressure')
plt.xlabel('Radius (km)')
plt.title('Neutron Star')
plt.legend(loc=0)
plt.show()

我不确定出了什么问题，但我收到了溢出错误。图表结果也是垃圾。压力， $P\left ( r \right )$ 也应该减少 $r$ ，但它会以非常尖锐的尖峰跳跃。对我做错的任何帮助将不胜感激。

import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint def dP_dr(P, r, p_c, R, G, c): return - G * (P + density(r, p_c, R) / c**2) * (mass(r, p_c, R) + 4 * np.pi * r**3 * P / c**2) / (r * (r - 2 * G* mass(r, p_c, R) / c**2)) def mass(r, p_c, R): return 4 * np.pi * p_c * (5 * r**3 * R**2 - 3 * r**5) / (15 * R**2) def density(r, p_c, R): return p_c * (1 - (r/R)**2) p_c = 7 * 10**10 #central density P_c = 1.6 * 10**15 #central pressure R = 10 #radius c = 3 * 10**8 G = 6.674 * 10**(-11) deltaR = 0.0001 r = np.linspace(deltaR, R, int(R/deltaR)) P, dic = odeint(dP_dr, P_c, r, args=(p_c, R, G, c), full_output=True, printmessg=True) num = 10000 plt.figure() plt.plot(r[:num], P[:num], label='Pressure') plt.xlabel('Radius (km)') plt.title('Neutron Star') plt.legend() plt.show()