Perché non viene chiamato il Dfunc (gradiente) durante l'utilizzo di integr.odeint in SciPy?

Question

Qualcuno può fornire un esempio di fornitura di un jacobiano a una funzione integrate.odeint in SciPy ?. Cerco di eseguire questo codice dal tutorial di SciPy odeint example ma sembra che Dfunc (gradiente) non venga mai chiamato.

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

Sotto il cofano, scipy.integrate.odeint utilizza il componente LSODA dal ODEPACK FORTRAN library. Per far fronte a situazioni in cui la funzione che stai cercando di integrare è stiff, LSODA passa in modo adattivo tra due metodi diversi per calcolare l'integrale - Adams' method, che è più veloce ma non adatto per sistemi rigidi, e BDF, che è più lento ma robusto alla rigidità .

La particolare funzione che si sta tentando di integrare non è rigida, quindi LSODA utilizzerà Adams a ogni iterazione. È possibile controllare questo restituendo il infodict (...,full_output=True) e controllando infodict['mused'].

Poiché il metodo di Adams non utilizza Jacobian, la funzione di gradiente non viene mai chiamata. Tuttavia, se si dà odeint una funzione rigida da integrare, come ad esempio il 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]) 

si dovrebbe vedere che odeint passa utilizzando BDF, e la funzione di Jacobiano ora viene chiamata.

Se si desidera un maggiore controllo sul risolutore, è necessario esaminare scipy.integrate.ode, che è un'interfaccia orientata agli oggetti molto più flessibile per più integratori diversi.