Dlaczego Dfunc (gradient) nie jest wywoływany podczas używania funkcji integate.odeint w SciPy?

Question

ktoś może stanowić przykład zapewniając jakobian do integrate.odeint funkcji w scipy ?. próbuję uruchomić ten kod z scipy tutorialu odeint example ale wydaje się, że Dfunc (gradient) nigdy nie jest wywoływana.

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

Pod maską scipy.integrate.odeint używa solver LSODA z ODEPACK FORTRAN library. W celu radzenia sobie z sytuacjami, w których funkcja staramy się zintegrować jest stiff, LSODA włącza adaptacyjny między dwa różne sposoby obliczania całki - Adams' method, co jest szybsze, ale nie nadają się do systemów sztywne i BDF, który jest wolniejszy, ale wytrzymała do sztywności .

Szczególna funkcja starasz się zintegrować nie jest sztywny, więc LSODA użyje Adamsa w każdej iteracji. Możesz to sprawdzić, zwracając infodict (...,full_output=True) i sprawdzając infodict['mused'].

Ponieważ metoda Adamsa nie używa Jacobiego, czynność gradientu nie jest wywoływana. Jednak jeśli dasz odeint sztywny funkcję integracji, takie jak 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]) 

należy zauważyć, że odeint przełącza się za pomocą BDF, a funkcja Jacobiego teraz jest wywoływana.

Jeśli chcesz mieć większą kontrolę nad solver, powinieneś zajrzeć do scipy.integrate.ode, który jest znacznie bardziej elastycznym interfejsem obiektowym dla wielu różnych integratorów.