为什么在SciPy中使用integrate.odeint时不调用Dfunc(渐变)？

Question

任何人都可以提供一个integrate.odeint在SciPy中为函数提供jacobian的例子吗？我尝试从SciPy教程odeint示例运行此代码,但似乎永远不会调用Dfunc(渐变).

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

在引擎盖下,scipy.integrate.odeint使用ODEPACK FORTRAN库中的LSODA解算器.为了处理你试图集成的函数是僵硬的情况,LSODA在两种不同的方法之间自适应地切换以计算积分 - Adams的方法,这种方法更快但不适合刚性系统,而BDF则更慢但是更健壮僵硬.

您尝试集成的特定功能是非僵硬的,因此LSODA将在每次迭代时使用Adams.您可以通过返回infodict(...,full_output=True)并检查来检查这一点infodict['mused'].

由于Adams的方法不使用雅可比行列式,因此您的渐变函数永远不会被调用.但是如果你给出odeint一个强大的功能来整合,比如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])

你应该看到odeint切换到使用BDF,现在调用雅可比函数.

如果您想要更多地控制求解器,您应该研究一下scipy.integrate.ode,这是一个更灵活的面向对象的多个不同集成器的接口.