Numerical issue in odeint function depending on the initial values

Question

I tried to debug my code (which solves coupled ODEs). I found odeint function has specific limits of taking decimal values because of this issue (may be) my plots are different for different initial conditions.

problem is in line 20: x = np.linspace(-30., -50., 5)
problem is in line 25: state0 = [gi,Li,60.]

When I have initial condition in x (-30) and in state0 (60) code work properly:

But when I change the initial conditions x (-25) and in state0 (50) code does not work:

I solved same code in Fortran it works for both limits. Could anyone suggest how to resolve this problem.

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as si
#function for solving differential equations:
def func1(state, x):
    g = state[0]
    L = state[1]
    u = state[2]
    phi1 = 1.645
    phi2 = 2.* 1.202
    N =  (-2.*g/(6.*np.pi + 5.*g))*(18./(1. - 2.*L) + 5.*np.log(1.- 2.*L) - phi1 + 6. )
    B = (-(2. - N)*L) - ((g/np.pi)* (5.*np.log(1.-2.*L) - phi2 + (5.*N/4.)))
    udx = -2.*(1.- ( (3.* 0.5* L)/(3.* 0.330))) #5.1
    gdx = (udx) * (2. + N)*g #5.1a
    Ldx = (udx) *( B)  #5.1b
    print gdx, Ldx, udx
    return gdx, Ldx, udx
gt = 10.**(-60)
Lt = (1.202*gt)/np.pi
x = np.linspace(-30., -50., 500) #initial conditions
for i in np.arange(len(x)):
    gi= np.float64(gt * np.exp(-4.*x[i]))
    Li = np.float64(Lt * np.cosh(4.*x[i]))
    break   
state0 = [gi,Li,60.] #initial conditions
G= si.odeint(func1, state0, x)
a = np.transpose(G) # a[0]=g; a[1]=lambda (L); a[2]=u 
plt.plot(a[1],a[0])
plt.title('RGIII Tragectory')
plt.xlabel('L ---->')
plt.ylabel('g ----->')
plt.show()

Answer 1

Indeed, the odeint solver is struggling with this ODE system. Apparently, the motion that this ODE describes is very non-uniform in time. A fix is to tell the solver to use tighter error control:

G = si.odeint(func1, state0, x, rtol=1e-10)

This yields a spiraling plot for both of your sets of parameters.

---

Generally, when odeint is acting strange, the first thing to do is to request its infodict:

G, info = si.odeint(func1, state0, x, full_output=True)

The infodict has a lot of information in it, but already a cursory look at info['hu'] (with the settings that caused your problem) shows ridiculous values of time step being used; it's as if the solver decided this solution isn't going anywhere and covered the requested interval in one time step; drawing a tiny line segment.