lumartor - 1 month ago 13

Python Question

Can anyone provide an example of providing a jacobian to a

`integrate.odeint`

I try to run this code from SciPy tutorial odeint example but seems that Dfunc (gradient) 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

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.

Source (Stackoverflow)

Comments