为什么在使用集成.odeint在学校里？

Question

有谁能提供一个例子来给SciPy中的integrate.odeint函数提供雅可比矩阵？。 我试图从SciPy教程odeint example运行这段代码，但似乎从未调用Dfunc（gradient）。在

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

在引擎盖下，^{}使用来自ODEPACK FORTRAN library的LSODA解算器。为了处理您试图积分的函数是stiff的情况，LSODA在计算积分的两种不同方法之间自适应地切换-Adams' method，后者更快，但不适合刚性系统；和{a5}，后者速度较慢，但对刚度鲁棒。在

您要集成的特定函数是非刚性的，因此LSODA将在每次迭代中使用Adams。您可以通过返回infodict（...,full_output=True）并检查infodict['mused']来检查这个问题。在

由于Adams的方法没有使用Jacobian，所以梯度函数永远不会被调用。但是，如果给odeint一个刚性函数来进行积分，例如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])

您应该看到odeint切换到使用BDF，现在调用了Jacobian函数。在

如果你想对解算器有更多的控制，你应该看看^{}，这是一个更灵活的面向对象的接口，可供多个不同的集成器使用。在