python - 尝试使用 Solve_ivp 求解牛顿多面体，并且正确性因 n 的不同值而异？

问题描述

我正在尝试为多方状态方程求解牛顿恒星结构方程（我正在求解多方体）。我没有使用 Lane-Emden 方程。这应该是一个非常简单的代码，因为它是两个方程的简单线性系统。状态方程也很简单。

由于堆栈溢出不接受 tex 并且它不会让我嵌入图像。我不确定如何以简洁的方式将牛顿方程放在这里，所以这是我能做到的最好的：

dm/dr = 4 pi rho r^2

dp/dr = - G rho M / r^2

对于 n = 1，存在多方体的精确解，并且我的压力半径曲线完全匹配，并且数值半径在精确解半径的 0.0003% 误差范围内。

然而，对于 n=3，我得到的质量是 1.43 个太阳质量而不是 1.24。

对于 n=3/2，我计算的所有半径都是公布版本的 3.7 倍，质量是公布结果的 28 倍。我不确定是什么原因造成的。

我已经用无量纲量的几何单位编写了代码，并且所有东西都是 SI，我得到的结果是一致的。这告诉我，这不是来自处理大量数字的错误。所以我将 SI 代码放在这里，这样就不会因为单位和比例因子的变化而混淆。

多面体计算的代码是这样的：

#for gamma = 4/3

K = ((hbar * c) /(12*np.pi**2.)) * ((3*np.pi**2.)/(m_h))**(4./3.)  
K = K * 0.5**(4./3.)


#for gamma = 5/3 

K_nr = hbar_cgs**2.
K_nr = K_nr /(15. *np.pi**2. * me_cgs) 
K_nr = K_nr * (3. * np.pi**2. )**(5./3.)
K_nr = K_nr * (mh_cgs)**(-5./3.)
K_nr = K_nr * 0.5**(5./3.)

#Equation of State

def EOS(p):

    rho = (p/K)**(1./gamma)

    return rho


def TOV(r,y):
    M = y[0]
    p = y[1]

    rho = EOS(p)
    #print p
    dMdr = 4. * np.pi * rho * r**2.
    dpdr = - G * M * rho /r**2.
    #print dpdr


    return [dMdr,dpdr]

def star_boundary(r,y):
    return y[1]

#Set star boundary at pressure = 0
star_boundary.terminal = True

M_0 = 0.
r_0 = 0.01 #m
r_stop = 20. #km
r_stop = r_stop * 10.**(3.) #SI = m
t_span = (r_0,r_stop)
t_eval = np.linspace(r_0,r_stop,1000)
p0 = 10**33.
y0 = [M_0, p0]

soln = solve_ivp(TOV,t_span,y0,method='RK45', events=star_boundary, t_eval=t_eval, dense_output=True)

r = soln.t
M = soln.y[0]
p = soln.y[1]

计算精确解的代码在这里：

rho0 = EOS(p0)
R = (((1.+1.)*p0 )/(4*np.pi * g_cgs * rho0**2))**0.5 * np.pi  
error = abs((r[-1] - R)/R)
print "percent error is ", error


R_s = (((1.+1.)*p0 )/(4*np.pi * g_cgs * rho0**2))**0.5
xi = r / R_s
theta = np.sin(xi) / xi
P = p0 * theta**(n+1)

我用 3 篇不同的论文检查了我的 K 值。我检查了输出 xi 是否等于 pi。在继续使用 GR 解决方案之前，我需要找到错误，但我很困惑。

我还检查了 r_0 的较小值（因为您实际上不能使用 r=0），我发现解决方案在这一点附近是稳定的。）

我还尝试降低积分器上的 rtol/atol，以防它只是累积错误，但将 rtol 从默认的 rtol=10E-3 更改为 10E-6 没有任何效果。

我还检查了 scipy.odeint

标签： pythonscipysolver

python - 尝试使用 Solve_ivp 求解牛顿多面体，并且正确性因 n 的不同值而异？

问题描述

解决方案

推荐阅读