Using Runge-Kutta to solve coupled differential equations

Question

I have a system of coupled equations: the hydrostatic equilibrium equation, the mass continuity equation, and an equation of state of the ideal gas. These are, in mathematical grammer,

\frac{dP}{dr}=- \rho*g,

where \rho is the density and g is the gravitational acceleration.

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

and

p=\rho* k_B* T/(\mu *m_p),

where k_B is boltzmann's constant, \mu is the mean molecular weight, and m_p is the proton mass.

I want to solve these coupled equations using the Runge-Kutta numerical technique, and I herein show the python code I have devised to solve this problem:

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)

for intial conditions for the pressure, p_atm, radius, r_atm, and mass, m_atm, I use the values I have shown in h, r_atm, p_atm, m_atm = rk4_1(g,f, 6.991e7, 1e-6*1e5, 1.898e27, 2.0e7,10000000). Notice that I am approaching this boundary-value problem from the upper atmosphere (where the initial conditions are given) and progressing downward in the atmosphere (notice that h is negative). My intention is to evaluate this numerical integration from 10^-1 Pascals to 100000 Pascals. The result I get from running this code is that the pressure simply blows up to ~1e+123 in three steps, so there is obviously something terribly wrong streaming about, but it would help to have another eye or perspective, for this is the first time I am performing Runga-Kutta methodology.

To make sure you have enough precision I would highly recommend using the decimal module so you can set your own rpecision: https://docs.python.org/2/library/decimal.html Looking at the code, I would think the `h = (r1 - r0)/n` line might already be problematic, it might do integer division which quickly kills any precision. — Wolph, Jun 30 '15 at 12:55

score 1 · Accepted Answer · answered Jun 30 '15 at 13:45

As Wolph says, dividing by n might simply give you h=0, depending on which version of Python you're using. If you're using 2.x, you should include from __future__ import division in the beginning, or handle the division in some other way (e.g., divide by float(n)). (Oh, and I guess perhaps you also intended to use n in your loop, rather than hard-coding range(0,10000000)? And there are a couple of indentation errors in the code as it stands, but I guess that's just from posting it here.)

This doesn't seem to be the main problem, though. You say you get a high pressure early; when I run it, it gets really low? Even with proper divisions, I get p_atm[3] = -2.27e+97, and from that, I start getting infinities (inf and -inf) and nans.

It's hard, without knowing the specific problem better, to see if there's an error in your implementation, or if this is simply a matter of numerical instability. It looks right to me, but I may very well have missed something (sort of hard to read.) If this is your first time with Runge–Kutta, I'd strongly suggest using an existing implementation rather than trying to get it right yourself. Numerical computation and avoiding floating-point issues can be quite challenging. You're already using scipy — why not use their implementation of the R–K method, or related numerical integration solutions? have a look at scipy.integrate, for example. If nothing else, if the scipy integrators can't solve your problem, at least you know more about what your challenges are.

thank you for the few minor corrections. However, the problem still remains, and wonder if anyone has additional insight into this. — inquiries, Jun 30 '15 at 15:28

Wolph · Answer 2 · 2015-07-01T13:06:01.663

Here's a version that uses decimals btw, it seems to work slightly better:

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

Thank you for your consideration again. I ran your modification, and this time we are able to get 13 steps before the code stops. There is something terribly wrong still, however, because I expect many steps and the mass itself blows up after the first step. See this by putting in `print 'm_atm', m_atm[i]` — inquiries, Jul 01 '15 at 10:31
Yes, I see. It goes outside of the `Decimal` range which is `999999999` by default. It could be increased but I don't know enough of this specific problem — Wolph, Jul 01 '15 at 13:09

score 0 · Answer 3 · answered Oct 13 '22 at 20:45

As f is the P derivative and g the M derivative function, then the k's are the slopes of M and the l's s the slopes of P. p_atm[i] + 0.5*k0 in k1 is thus using the wrong offset, it should be

p_atm[i] + 0.5*l0,

as is done in the next line for l1.

The effect this has on the result is unpredictable. For sufficiently small step sizes it just reduces the order of the method to one. For larger step sizes it could make the integration unstable (where RK4 is still stable), with chaotic results

Using Runge-Kutta to solve coupled differential equations

3 Answers3