Пытаетесь решить ньютоновские многогранники с помощью Solve_ivp, и корректность меняется для разных значений n? - PullRequest
Новая
       44

Пытаетесь решить ньютоновские многогранники с помощью Solve_ivp, и корректность меняется для разных значений n?

1 голос
/ 17 апреля 2020

Я пытаюсь решить ньютоновские уравнения звездной структуры для политропи c уравнения состояния (я решаю политропы). Я не использую уравнения Лейна-Эмдена. Это должен быть очень простой код, поскольку это простая линейная система из двух уравнений. Уравнение состояния также простое.

Поскольку переполнение стека не принимает текс и не позволяет мне вставить изображение. Я не уверен, как аккуратно представить ньютоновские уравнения здесь, поэтому вот лучшее, на что я способен:

дм / др = 4 пи ро р ^ 2

дп / dr = - G rho M / r ^ 2

Для n = 1 есть точное решение политропа, и моя кривая давления-радиуса точно совпадает, а числовой радиус находится в пределах 0,0003% ошибки точного радиуса решения.

Однако при n = 3 я получаю массу 1,43 солнечной массы вместо 1,24.

При n = 3/2 все мои рассчитанные радиусы в 3,7 раза превышают опубликованную версию, а массы в 28 раз превышают опубликованные результаты. Я не уверен, что вызвало это.

Я сделал код в геометрических единицах c единиц с безразмерными величинами и со всем в СИ, и результаты, которые я получаю, согласуются. Что говорит мне, что это не из-за ошибок при работе с большими числами. Поэтому я помещу здесь код 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)

Я проверил свое значение K с 3 различными бумагами. Я проверил, что выходной XI равен пи. Мне нужно найти ошибку до того, как я go перейду к решению GR, и я застрял.

Я также проверил меньшие значения r_0 (поскольку вы фактически не можете использовать r = 0), и я обнаружил, что решение устойчиво в этой точке.)

Я также попытался понизить rtol / atol на интеграторе, если он просто накапливал ошибку, но изменение rtol со значения по умолчанию rtol = 10E-3 на 10E-6 ничего не сделало.

Я тоже проверял с scipy.odeint

...