我有一个耦合方程组:流体静力平衡方程、质量连续性方程和理想气体的状态方程。这些是,在数学语法中,
\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
是质子质量。
我想使用 Runge-Kutta 数值技术解决这些耦合方程,我在这里展示了我为解决这个问题而设计的 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
在三个步骤中简单地爆炸,所以显然有一些非常错误的流媒体,但它会有助于有另一个眼睛或视角,因为这是我第一次表演Runga-Kutta 方法。