尝试用solve_ivp求解牛顿多面体，其正确性因n的不同值而不同？

Question

我在尝试解牛顿的恒星结构方程，解多向性状态方程（我在解多向性）。我没有用莱恩·埃姆登方程。这应该是一个非常简单的代码，因为它是由两个方程组成的简单线性系统。状态方程也很简单

因为堆栈溢出不接受tex，也不允许我嵌入图像。我不知道如何将牛顿方程整齐地放在这里，因此我只能做到以下几点：

dm/dr=4πρ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值。我已经检查了输出席等于π。我需要在继续GR解决方案之前找到错误，我被卡住了

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

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

我还检查了scipy.odeint

Answer 1

更新：原来我比较的出版物中有一个错误。结果表明，大约1.43是正确的（它与钱德拉塞卡质量有关）