Tại sao không phải là Dfunc (gradient) được gọi trong khi sử dụng integration.odeint trong SciPy?

Question

Bất kỳ ai cũng có thể cung cấp một ví dụ về việc cung cấp một jacobian cho một hàm integrate.odeint trong SciPy ?. Tôi cố gắng chạy mã này từ hướng dẫn SciPy odeint example nhưng dường như Dfunc (gradient) không bao giờ được gọi.

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

Dưới mui xe, scipy.integrate.odeint sử dụng bộ giải LSODA từ ODEPACK FORTRAN library. Để xử lý các tình huống mà chức năng bạn đang tích hợp là stiff, LSODA chuyển đổi thích nghi giữa hai phương pháp khác nhau để tính tích phân - Adams' method, nhanh hơn nhưng không phù hợp với hệ thống cứng và BDF. .

Chức năng cụ thể mà bạn đang cố gắng tích hợp là không cứng, vì vậy LSODA sẽ sử dụng Adams trên mọi lần lặp lại. Bạn có thể kiểm tra điều này bằng cách trả lại infodict (...,full_output=True) và kiểm tra infodict['mused'].

Vì phương pháp Adams 'không sử dụng Jacobian, hàm gradient của bạn không bao giờ được gọi. Tuy nhiên nếu bạn cho odeint một hàm cứng tích hợp, chẳng hạn như 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]) 

bạn sẽ thấy rằng odeint chuyển sang sử dụng BDF, và các chức năng Jacobian nay được gọi.

Nếu bạn muốn kiểm soát nhiều hơn đối với bộ giải, bạn nên xem xét scipy.integrate.ode, đây là giao diện hướng đối tượng linh hoạt hơn nhiều cho nhiều trình tích hợp khác nhau.