scipy 用龙格-库塔法求解耦合微分方程

mwg9r5ms 于 2022-11-09 发布在其他

关注(0)|答案(3)|浏览(292)

我有一个耦合方程组：流体静力学平衡方程、质量连续性方程和理想气体的状态方程。用数学语法来说，这些方程是：

\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，所以很明显这里有一些非常不对劲的地方，但是如果有另一只眼睛或者另一个视角会有所帮助，因为这是我第一次使用龙加-库塔方法。

scipy

来源：https://stackoverflow.com/questions/31101172/using-runge-kutta-to-solve-coupled-differential-equations

3条答案

按热度按时间

8ehkhllq1#

正如Wolph所说，除以n可能只给予h=0，这取决于你使用的Python版本。如果你使用的是Python 2.x，你应该在开头包含from __future__ import division，或者以其他方式处理除法（例如，除以float(n)）。（哦，我猜您可能还打算在循环中使用n，而不是硬编码range(0,10000000)？代码中有一些缩进错误，但我想这只是因为在这里发布了它。）
这似乎不是主要的问题，你说你很早就有了高血压;当我运行它的时候，它会变得很低，即使用正确的除法，我也会得到p_atm[3] = -2.27e+97，然后，我开始得到无穷大（inf和-inf）和nan。
在不更好地了解具体问题的情况下，很难看出实现中是否存在错误，或者这仅仅是数值不稳定的问题。如果这是你第一次使用龙格-库塔法，我强烈建议使用现有的实现，而不是尝试自己去做。数值计算和避免浮点问题可能是相当具有挑战性的。如果您已经在使用scipy-为什么不使用他们的R-K方法实现或相关的数值积分解决方案？例如，看看scipy.integrate。如果scipy积分器不能解决您的问题，至少您知道更多关于您的挑战是什么。

赞(0）回复(0）举报 2022-11-09

x0fgdtte2#

这里有一个版本，使用小数顺便说一句，它似乎工作得稍微好：

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

赞(0）回复(0）举报 2022-11-09

ajsxfq5m3#

由于f是P的导数，g是M的导数函数，因此k是M的斜率，l是P的斜率。因此，k1中的p_atm[i] + 0.5*k0使用了错误的偏移量，它应该是