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.title('quintessence cosmological model')

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$.

z = [t0]
i = 0;

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

  • $\begingroup$ Thank you for helping. I changed the increments( edited in the ques too) but still end up getting the same result unless I reduce the initial value of v which is 1101 here to a single digit number. I know that we can scale things but am not sure how to do it in this case so that the other equations are not affected by scaling the value of v. And I am relatively new in using python,I understood the concept of using one vector valued ODE function, but am not sure what or how I should write the RK4 argument after that. Any suggestions here? $\endgroup$ – tsv Jul 3 '17 at 5:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.