[Math] Runge-Kutta 4th order

Question

I am not sure if it's a mathematics question or coding part although I am pretty sure that the code I used is right,so do excuse me if I a wrong.

I have a set of three first order ode and I am trying to numerically integrate them in python using RK4 method. The problem arises with the exponential term where python rounds it to $0$ and the values thereafter are returned "NAN". I tried various packages to deal with large values still no use. I am thinking if I modify the equations or scale them so that the value is smaller might help, but I am not sure how to do that. Any help or the reference to that would be really helpful as I am kind of in a really stuck up situation.
Here's the code I have used for this purpose:

    #simple exponential potential
# u = K*phi'/H0; v = omega_matter**(1/3)*(1+z); w = l*K*phi' - ln((K**2)*V0/H0**2),z from 1100(initial) to z=0(final)
# f,g,h are functions of derivation of u,v,w respectively derived w.r.t t*H0 = T

import matplotlib.pyplot as plt
import numpy as np
import math

def f(u,v,w):
    return -3*u*((v**3 + (u**2)/6 + np.exp(-w)/3)**0.5) + l*np.exp(-w)

def g(u,v,w):
    return -v*(v**3 + (u**2)/6 + np.exp(-w)/3)**0.5

def h(u):
    return l*u

z = np.linspace(start=0.0,stop=1.0,num = 10001)
print (z)    
p = 0.1
q = 1.0
n = 10000.0
dh = (q-p)/n
u = [0.0]
v = [1101]
w = [1.8]
l = 1.0


for i in range(0,int(n)):
k1 = f(u[i],v[i],w[i])
r1 = g(u[i],v[i],w[i])
s1 = h(u[i])
k2 = f(u[i] + k1*dh/2,v[i] + r1*dh/2,w[i] + s1*dh/2)
r2 = g(u[i] + k1*dh/2,v[i] + r1*dh/2,w[i] + s1*dh/2)
s2 = h(u[i] + k1*dh/2)
k3 = f(u[i] + k2*dh/2,v[i] + r2*dh/2,w[i] + s2*dh/2)
r3 = g(u[i] + k2*dh/2,v[i] + r2*dh/2,w[i] + s2*dh/2)
s3 = h(u[i] + k2*dh/2)
k4 = f(u[i] + dh*k3,v[i] + dh*r3,w[i] + s3*dh)
r4 = g(u[i] + k3*dh,v[i] + dh*r3,w[i] + s3*dh)
s4 = h(u[i] + dh*k3)
u == u.append(u[i] + (dh/6)*(k1 + 2.0*k2 + 2.0*k3 + k4))
v == v.append(v[i] + (dh/6)*(r1 + 2.0*r2 + 2.0*r3 + r4))
w == w.append(w[i] + (dh/6)*(s1 + 2.0*s2 + 2.0*s3 + s4))

plt.plot(z,u, '-b')
plt.plot(z,v, '-r')
plt.plot(z,w, '-g') 
#plt.plot(u,v)
#plt.plot(u,w)  
plt.title('quintessence cosmological model')
plt.show()

Answer 1

The immediate problem is that you are using the increments of the stages in a wrong way. With

k1 = f(u[i],v[i],w[i])
r1 = g(u[i],v[i],w[i])
s1 = h(u[i])

k1 is always the increment of u, r1 always the increment of v and s1 always the increment of w. The second stage should thus read

k2 = f(u[i] + k1*dh/2,v[i] + r1*dh/2,w[i] + s1*dh/2)
r2 = g(u[i] + k1*dh/2,v[i] + r1*dh/2,w[i] + s1*dh/2)
s2 = h(u[i] + k1*dh/2)

and so on.

As you are using numpy it would be in general better to have one vector valued ODE function so that the RK4 step has only 5 lines.

def derivs(t,y):
    u,v,w = y
    emw = np.exp(-w)
    term1 = (v**3 + (u**2)/6 + emw/3)**0.5
    return np.array([-3*u*(term1) + l*emw, -v*term1, l*u ]);

---

(additional comments jul/24) The initial slope for the $v$ component is about $−v^{5/2}$ which for $v∼10^3$ gives $g∼−10^{7.5}$. If you want qualitatively sensible values (saying nothing about quantitative correctness) for the second stage of RK4, you need $dh⋅10^{7.5}≤10^3$. To get somewhat quantitatively correct values, this suggests $n≥10^5$. With n=1e5 I got a sensible looking graph.

The initially large slope then falling towards a horizontal asymptote also strongly suggests using an adaptive step size. Roughly estimating the local Lipschitz constant as being about $\sim v^{3/2}$ in the initial phase, RK4 in double precision works best if $dh·v^{3/2}\sim 10^{-3}$, so something like dh = h0/(1+abs(v[i])**1.5), $h_0>10^{-3}$ will distribute the numerical error somewhat evenly over the integration interval in about $7/h_0$ integration steps. This even works with $h_0=0.5$.

t0=0.0
tf=1.0
h0=0.5
z = [t0]
i = 0;


while z[i] < tf:
    dh = h0/(1+abs(v[i])**1.5)
    ...

    z.append(z[i]+dh);
    i+=1;