dimanche 19 juin 2016

Calculating decay using scipy.integrate.odeint()

I'm trying to calculate the decay of a few species using scipy.integrate.odeint(). One of the decay coefficients is very large, and that seems to be causing an issue and I'm not sure why or how to fix it. Here is my code:

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

def decay(init,t):
    a,b,c,d = init
    f0 = - 6*a
    f1 = 3 * a - 1e12*b
    f2 = 3 * a
    f3 = 1e12 * b
    return [f0,f1,f2,f3]

if __name__ == '__main__':
    init = [1,0,0,0]
    time = np.linspace(0, 1e12, 100001)
    soln  = odeint(decay, init ,time)

which gives me the following error:

lsoda--  at t (=r1) and step size h (=r2), the    
       corrector convergence failed repeatedly       
       or with abs(h) = hmin   ls
      in above,  r1 =  0.0000000000000D+00   r2 =  0.1552206225945D-09
/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/integrate/odepack.py:218: ODEintWarning: Repeated convergence failures (perhaps bad Jacobian or tolerances). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
Jakes-MacBook-Pro:Radiation_Calculator jakehuneau$ python tester.py
 lsoda--  at t (=r1) and step size h (=r2), the    
       corrector convergence failed repeatedly       
       or with abs(h) = hmin   ls
      in above,  r1 =  0.0000000000000D+00   r2 =  0.1552206312178D-09
/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/integrate/odepack.py:218: ODEintWarning: Repeated convergence failures (perhaps bad Jacobian or tolerances). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)

The code does work however if the coefficients on b in the decay function are 1e10 instead of 1e12.

Aucun commentaire:

Enregistrer un commentaire