Why isn't the Dfun(gradient) called while using integrate.odeint in SciPy?

Question

Can anyone provide an example of providing a Jacobian to a integrate.odeint function in SciPy?. I try to run this code from SciPy tutorial odeint example but seems that Dfun() (the Jacobian function) is never called.

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):
    return [t*y[1],y[0]]
    
    
def gradient(y,t):
    print 'jacobian'  # added
    return [[0,t],[1,0]]
    
    
x = arange(0,4.0, 0.01)
t = x
ychk = airy(x)[0]
y = odeint(func, y0, t)
y2 = odeint(func, y0, t, Dfun=gradient)
print y2 # added

Answer 1

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 infodict (...,full_output=True) and checking infodict['mused'].

Since Adams' method does not use the Jacobian, your gradient function never gets called. However if you give odeint a stiff function to integrate, such as the Van der Pol equation:

def vanderpol(y, t, mu=1000.):
    return [y[1], mu*(1. - y[0]**2)*y[1] - y[0]]

def vanderpol_jac(y, t, mu=1000.):
    return [[0, 1], [-2*y[0]*y[1]*mu - 1, mu*(1 - y[0]**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.