Использование Рунге-Кутта для решения связанных дифференциальных уравнений
У меня есть система связанных уравнений: уравнение гидростатического равновесия, уравнение неразрывности массы и уравнение состояния идеального газа. Это, в математической грамматике,
\frac{dP}{dr}=- \rho*g,
где \rho это плотность и g это гравитационное ускорение.
\frac{dM}{dr}=4*pi* r^2*\rho
а также
p=\rho* k_B* T/(\mu *m_p),
где k_B постоянная Больцмана, \mu средняя молекулярная масса, и m_p это масса протона.
Я хочу решить эти связанные уравнения, используя числовую технику Рунге-Кутты, и здесь я покажу код Python, который я разработал для решения этой проблемы:
from scipy.constants import m_p,G,k,pi
from pylab import *
#mu may be changed for different molecular composition:
mu=2
def g(r_atm, p_atm):
T=165
return 4*pi*r_atm**2*mu*m_p*p_atm/(k*T)
def f(r_atm,p_atm, m_atm):
T=165
return -mu*m_p*p_atm*G*m_atm/(k*T*r_atm**2)
def rk4_1(g,f, r0, p0,m0, r1, n):
r_atm = [0]*(n + 1)
p_atm = [0]*(n + 1)
m_atm=[0]*(n + 1)
h = (r1 - r0)/n
# h=-20
r_atm[0]=r0
p_atm[0]=p0
m_atm[0]=m0
for i in range(0,10000000):
if p_atm[i]<100000:
k0 = h*g(r_atm[i], p_atm[i])
l0 = h*f(r_atm[i], p_atm[i], m_atm[i])
k1 = h*g(r_atm[i] + 0.5*h, p_atm[i] + 0.5*k0)
l1 = h*f(r_atm[i] + 0.5*h, p_atm[i] + 0.5*l0, m_atm[i]+0.5*k0)
k2 = h*g(r_atm[i] + 0.5*h, p_atm[i] + 0.5*k1)
l2 = h*f(r_atm[i] + 0.5*h, p_atm[i] + 0.5*l1, m_atm[i]+0.5*k1)
k3 = h*g(r_atm[i] + h, p_atm[i] + k2)
l3 = h*f(r_atm[i] + h, p_atm[i] + l2, m_atm[i]+k2)
r_atm[i+1] = r0 + (i+1)*h
p_atm[i+1] = p_atm[i] + (l0 + 2*l1 + 2*l2 + l3)/6
m_atm[i+1] = m_atm[i] + (k0 + 2*k1 + 2*k2 + k3)/6
else:
break
return h, r_atm, p_atm, m_atm
h, r_atm, p_atm, m_atm = rk4_1(g,f, 6.991e7, 1e-6*1e5, 1.898e27, 2.0e7,10000000) #bar to pascals (*1e5)
для начальных условий для давления, p_atm, радиус, r_atmи масса, m_atmЯ использую значения, которые я показал в h, r_atm, p_atm, m_atm = rk4_1(g,f, 6.991e7, 1e-6*1e5, 1.898e27, 2.0e7,10000000), Обратите внимание, что я подхожу к этой краевой задаче из верхней атмосферы (где заданы начальные условия) и продвигаюсь вниз в атмосфере (обратите внимание, что h отрицательно). Мое намерение состоит в том, чтобы оценить эту числовую интеграцию из 10^-1 Паскали в 100000 Па. В результате выполнения этого кода я получаю, что давление просто поднимается до ~1e+123 в трех шагах, поэтому, очевидно, что-то ужасно неправильно течет, но это помогло бы иметь другой взгляд или перспективу, потому что я впервые выполняю методологию Рунга-Кутты.
3 ответа
Как говорит Вольф, деление на n может просто дать вам h=0в зависимости от того, какую версию Python вы используете. Если вы используете 2.x, вы должны включить from __future__ import division в начале, или обрабатывать деление каким-либо другим способом (например, делить на float(n)). (О, и я думаю, возможно, вы также намеревались использовать n в вашем цикле, а не жесткое кодирование range(0,10000000)? И в коде есть пара ошибок отступа, но я думаю, что это просто из-за публикации здесь.)
Это, похоже, не главная проблема. Вы говорите, что рано получаете высокое давление; когда я запускаю его, он становится действительно низким? Даже при правильном делении я получаю p_atm[3] = -2.27e+97и от этого я начинаю получать бесконечности (inf а также -inf) а также nans.
Трудно, не зная конкретной проблемы, лучше понять, есть ли ошибка в вашей реализации или это просто вопрос числовой нестабильности. Мне это кажется правильным, но я вполне мог что-то упустить (что-то трудно читаемое). Если вы впервые используете Runge-Kutta, я бы настоятельно рекомендовал использовать существующую реализацию, а не пытаться сделать это самостоятельно, Численные вычисления и избежание проблем с плавающей точкой могут быть довольно сложными. Вы уже используете scipy - Почему бы не использовать их реализацию метода R–K или связанных с ними решений для численного интегрирования? взгляните на scipy.integrate, например. Если ничего другого, если scipy Интеграторы не могут решить вашу проблему, по крайней мере, вы знаете больше о своих проблемах.
Вот версия, которая использует десятичные знаки между прочим, кажется, работает немного лучше:
from decimal import Decimal as D
from scipy.constants import m_p,G,k,pi
m_p = D(m_p)
G = D(G)
k = D(k)
pi = D(pi)
# mu may be changed for different molecular composition:
mu = D(2)
def g(r_atm, p_atm):
T = D(165)
return D(4) * pi * r_atm ** D(2) * mu * m_p * p_atm/(k * T)
def f(r_atm,p_atm, m_atm):
T = D(165)
return -mu * m_p * p_atm * G * m_atm/(k * T * r_atm ** D(2))
def rk4_1(g,f, r0, p0,m0, r1, n):
r_atm = [D(0)] * (n + 1)
p_atm = [D(0)] * (n + 1)
m_atm = [D(0)] * (n + 1)
h = (r1 - r0) / n
# h = -20
r_atm[0] = r0
p_atm[0] = p0
m_atm[0] = m0
for i in range(0, 10000000):
if p_atm[i] < 100000:
k0 = h * g(r_atm[i], p_atm[i])
l0 = h * f(r_atm[i], p_atm[i], m_atm[i])
k1 = h * g(r_atm[i] + D('0.5') * h, p_atm[i] + D('0.5') * k0)
l1 = h * f(r_atm[i] + D('0.5') * h, p_atm[i] + D('0.5') * l0,
m_atm[i]+D('0.5') * k0)
k2 = h * g(r_atm[i] + D('0.5') * h, p_atm[i] + D('0.5') * k1)
l2 = h * f(r_atm[i] + D('0.5') * h, p_atm[i] + D('0.5') * l1,
m_atm[i]+D('0.5') * k1)
k3 = h * g(r_atm[i] + h, p_atm[i] + k2)
l3 = h * f(r_atm[i] + h, p_atm[i] + l2, m_atm[i]+k2)
r_atm[i + 1] = r0 + (i + 1) * h
p_atm[i + 1] = p_atm[i] + (l0 + D('2') * l1 + D('2') * l2 +
l3)/D('6')
m_atm[i + 1] = m_atm[i] + (k0 + D('2') * k1 + D('2') * k2 + k3)/D('6')
else:
break
return h, r_atm, p_atm, m_atm
h, r_atm, p_atm, m_atm = rk4_1(
g,
f,
D('6.991e7'),
D('1e-6') * D('1e5'),
D('1.898e27'),
D('2.0e7'),
10000000,
) # bar to pascals (*1e5)
print 'h', h
Какfявляется производной иgпроизводная функция, тоkэто склоныMиlна склонахP.p_atm[i] + 0.5*k0вk1таким образом, использует неправильное смещение, оно должно быть
p_atm[i] + 0.5*l0,
как это делается в следующей строке дляl1.
Влияние этого на результат непредсказуемо. Для достаточно малых размеров шага он просто снижает порядок метода до единицы. Для больших размеров шага это может сделать интеграцию нестабильной (где RK4 все еще стабилен) с хаотическими результатами.