RuntimeWarning: overflow encountered in exp

Question

I have three coupled ODEs and I am using RK4 method in python to solve them with the given initial conditions. When I run them I am getting the following error:

RuntimeWarning: overflow encountered in exp

Any help on where I might have done wrong would be highly appreciated.

Here's the code (originally pasted here):

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

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

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

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

def h(T,u):
    return l*u

p = 0.1
q = 1.0
dh = 0.1
n = (q-p)/dh
u = [0.0]
v = [1100]
T = [0.00001]
w = [-44.927]
l = 1.3

for i in range(0,int(n)):
    k1 = f(T[i],u[i],v[i],w[i])
    r1 = g(T[i],u[i],v[i],w[i])
    print k1, r1
    s1 = h(T[i],u[i])
    print s1
    k2 = f(T[i] + 0.5*dh,u[i] + k1*dh,v[i] + k1*dh,w[i] + k1*dh)
    r2 = g(T[i] + 0.5*dh,u[i] + r1*dh,v[i] + r1*dh,w[i] + r1*dh)
    s2 = h(T[i] + 0.5*dh,u[i] + s1*dh)
    print k2,r2,s2
    k3 = f(T[i] + 0.5*dh,u[i] + k2*dh,v[i] + k2*dh,w[i] + k2*dh)
    r3 = g(T[i] + 0.5*dh,u[i] + r2*dh,v[i] + r2*dh,w[i] + r2*dh)
    s3 = h(T[i] + 0.5*dh,u[i] + s2*dh)
    k4 = f(T[i] + dh,u[i] + dh*k3,v[i] + dh*k3,w[i] + k3*dh)
    r4 = g(T[i] + dh,u[i] + r3*dh,v[i] + dh*r3,w[i] + r3*dh)
    s4 = h(T[i] + dh,u[i] + dh*s3)
    T == T.append(T[i] + dh)
    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(T,u, '-b')
plt.plot(T,v, '-r')
plt.plot(T,w, '-g')   
plt.title('quintessence cosmological model')
plt.show()

Welcome! It might be useful to try asking this question in a slightly different way. It's best to provide the *minimal* possible code that reproduces the problem in question and to provide it in the body of the question. Questions that just ask for help in debugging a large section of code aren't really a good fit for SO. You might check out https://stackoverflow.com/help/how-to-ask for a guide on how best to ask questions here. — wphicks, Jun 13 '17 at 18:37
Possible duplicate of [overflow in exp, python](https://stackoverflow.com/questions/43100902/overflow-in-exp-python) — Knowledge Cube, Jun 13 '17 at 19:16
I strongly suspect that you have your equations for `f`, `g` and `h` wrong. Either that, or the problem is _extremely_ ill conditioned. The derivative of `v` at your starting point has value `-3619146406000.8389`, which already seems large, and after one tiny step you're now feeding enormous (negative) values into an `np.exp(-w)`, and getting infinities as a result. (I don't see anything wrong with the RK4 implementation.) Do you have a reference to the physical problem you're solving? — Mark Dickinson, Jun 13 '17 at 19:22
Voting to close: this isn't an issue with the code, which is solving the ODE just fine. It's an issue with the equations. (It may possibly be an issue with the way that the original equations were translated into Python, but we can't help with that unless we see those equations.) — Mark Dickinson, Jun 14 '17 at 05:43

Knowledge Cube · Accepted Answer · 2017-06-13T19:55:39.910

2

Your problem seems like it would originate in your f or g function calls. Watch your value of -w when it is passed to the numpy.exp function in that code. If it evaluates to a large enough positive number, calling the exponential function on that would seem like it could easily cause an overflow. Google tells me that e^44.927 would evaluate to 3.2474927e+19.

You may need to choose a smaller absolute value for w to run this code without getting an overflow error.

(See Python-How to determine largest/smallest int/long/float/complex numbers my system can handle for how you can determine what the largest allowed integer values are in Python.)

edited Jun 13 '17 at 19:55

answered Jun 13 '17 at 19:15

Knowledge Cube

990
12
35

Minimum and maximum values for ints aren't really relevant here; this is all floating-point. – Mark Dickinson Jun 13 '17 at 19:28
@MarkDickinson Good catch. I linked to a better question. – Knowledge Cube Jun 13 '17 at 19:56
Thank you for the help. The code used is correct but the problem is with the differential equations used. They are autonomous equations which would blow up after a certain time if the initial value is huge. Any idea or reference on how to compute autonomous equations cause I am not able to find any with proper explanation in net. @MarkDickinson – anon_particle Jun 16 '17 at 21:34

RuntimeWarning: overflow encountered in exp

1 Answers1