使用Runge-Kutta求解耦合微分方程

Question

我有一个耦合方程系统：流体静力平衡方程，质量连续性方程和理想气体状态方程。这些是，在数学语法，

1. \frac{dP}{dr}=- \rho*g，

其中\rho是密度和g是重力加速度。

\frac{dM}{dr}=4*pi* r^2*\rho

和

1. 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 INTIAL条件半径，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方法。

Answer 1

正如沃尔夫所说，除以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。

很难，不知道具体问题更好，看你的实现是否有错误，或者这只是一个数值不稳定问题。它看起来对我来说，但我可能错过了某些东西（有点难以阅读。）如果这是你第一次与龙格库塔，我强烈建议使用现有的实现，而不是试图得到它对你自己。数值计算和避免浮点问题可能非常具有挑战性。您已经在使用scipy - 为什么不使用它们的R-K方法或相关数值集成解决方案？例如，看看scipy.integrate。如果没有别的，如果scipy集成商不能解决你的问题，至少你更了解你的挑战是什么。

Answer 2

这是一个使用小数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