Can anyone provide an example of providing a jacobian to a
from numpy import * # added
from scipy.integrate import odeint
from scipy.special import gamma, airy
y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)
y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)
y0 = [y0_0, y1_0]
def func(y, t):
print 'jacobian' # added
x = arange(0,4.0, 0.01)
t = x
ychk = airy(x)
y = odeint(func, y0, t)
y2 = odeint(func, y0, t, Dfun=gradient)
print y2 # added
Under the hood,
scipy.integrate.odeint uses the LSODA solver from the ODEPACK FORTRAN library. In order to deal with situations where the function you are trying to integrate is stiff, LSODA switches adaptively between two different methods for computing the integral - Adams' method, which is faster but unsuitable for stiff systems, and BDF, which is slower but robust to stiffness.
The particular function you're trying to integrate is non-stiff, so LSODA will use Adams on every iteration. You can check this by returning the
...,full_output=True) and checking
def vanderpol(y, t, mu=1000.): return [y, mu*(1. - y**2)*y - y] def vanderpol_jac(y, t, mu=1000.): return [[0, 1], [-2*y*y*mu - 1, mu*(1 - y**2)]] y0 = [2, 0] t = arange(0, 5000, 1) y,info = odeint(vanderpol, y0, t, Dfun=vanderpol_jac, full_output=True) print info['mused'] # method used (1=adams, 2=bdf) print info['nje'] # cumulative number of jacobian evaluations plot(t, y[:,0])
you should see that
odeint switches to using BDF, and the Jacobian function now gets called.
If you want more control over the solver, you should look into
scipy.integrate.ode, which is a much more flexible object-oriented interface to multiple different integrators.