我有一个耦合方程组：流体静力平衡方程、质量连续性方程和理想气体的状态方程。在数学语法中，这些是
- CCD_ 1
其中rho是密度，g是重力加速度。
frac{dM}{dr}=4*pi* r^2*rho
和
- CCD_ 5
其中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、半径frac{dP}{dr}=- rho*g0和质量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，因此流媒体显然存在严重的错误，但换一种眼光或视角会有所帮助，因为这是我第一次执行Runga-Kutta方法论。
正如Wolph所说，除以n可能只会得到h=0，这取决于您使用的Python版本。如果使用2.x，则应该在开头包含from __future__ import division，或者以其他方式处理除法（例如，除以float(n)）。（哦，我想你可能也打算在循环中使用n，而不是硬编码range(0,10000000)？目前的代码中有几个缩进错误，但我想这只是因为把它发布在这里。）
不过，这似乎不是主要问题。你说你很早就有高压；当我运行它时，它真的很低吗？即使有适当的除法，我也得到p_atm[3] = -2.27e+97，从中，我开始得到无穷大（inf和-inf）和nan s。
如果不更好地了解具体问题，很难判断您的实现中是否存在错误，或者这是否只是数值不稳定的问题。在我看来，这是，但我很可能错过了一些东西（有点难以阅读。）如果这是你第一次使用Runge–Kutta，我强烈建议使用现有的实现，而不是自己尝试正确的实现。数值计算和避免浮点问题可能非常具有挑战性。您已经在使用scipy了——为什么不使用他们的R–K方法或相关的数值积分解决方案呢？例如，看看scipy.integrate。如果没有别的，如果scipy集成商不能解决你的问题，至少你知道你面临的挑战是什么。
这里有一个使用小数btw的版本，它似乎工作得稍微好一点：
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是P导数函数，g是M导数函数，则k是M的斜率，l是P的斜率。k1中的p_atm[i] + 0.5*k0因此使用了错误的偏移量，它应该是
p_atm[i] + 0.5*l0,
如在CCD_ 38的下一行中所做的那样。
这对结果的影响是不可预测的。对于足够小的步长，它只是将方法的阶数减少到一。对于较大的步长，它可能会使积分不稳定（其中RK4仍然稳定），产生混乱的结果