python - library - ¿Por qué no se llama a Dfunc(degradado) mientras se usa integrate.odeint en SciPy?

solve differential equation odeint

Question

¿Alguien puede dar un ejemplo de proporcionar un jacobian a una función integra.odeint en SciPy ?. Intento ejecutar este código desde el ejemplo de SciPy tutorial odeint pero parece que Dfunc (gradiente) nunca se llama.

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

Bajo el capó, scipy.integrate.odeint utiliza el solucionador LSODA de la biblioteca ODEPACK FORTRAN . Para lidiar con situaciones en las que la función que intenta integrar es rígida , LSODA conmuta adaptativamente entre dos métodos diferentes para calcular la integral: el método de Adams , que es más rápido pero no apto para sistemas rígidos, y BDF , que es más lento pero robusto a la rigidez.

La función particular que intentas integrar no es rígida, por lo que LSODA usará Adams en cada iteración. Puede verificar esto devolviendo el infodict ( ...,full_output=True ) y revisando infodict[''mused''] .

Como el método de Adams no usa el Jacobiano, nunca se llamará a su función de gradiente. Sin embargo, si le das a odeint una función rígida para integrar, como la ecuación de Van der Pol :

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])

debería ver que odeint cambia a usar BDF, y ahora se llama a la función jacobiana.

Si desea tener más control sobre el solucionador, debe buscar en scipy.integrate.ode , que es una interfaz mucho más flexible orientada a objetos para múltiples integradores diferentes.