为什么在 SciPy 中使用集成.odeint 时不调用 Dfun(gradient)?

问题描述

谁能提供一个为 SciPy 中的 integrate.odeint 函数提供雅可比行列式的示例?我尝试从 SciPy 教程中运行此代码 odeint 示例 但似乎Dfun()(雅可比函数)从未被调用.

Can anyone provide an example of providing a Jacobian to a integrate.odeint function in SciPy?.I try to run this code from SciPy tutorial odeint example but seems that Dfun() (the Jacobian function) is never called.

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

推荐答案

在幕后，scipy.integrate.odeint 使用来自 ODEPACK FORTRAN 库.为了处理您尝试集成的函数 stiff 的情况，LSODA 在计算积分的两种不同方法之间自适应切换 - Adams 方法，速度更快但不合适对于刚性系统，以及 BDF，其速度较慢，但​​对刚度具有鲁棒性.

Under the hood, scipy.integrate.odeint uses the LSODA solver from the ODEPACK FORTRAN library. In order to deal with situations where the function you are trying to integrate is stiff, LSODA switches adaptively between two different methods for computing the integral - Adams' method, which is faster but unsuitable for stiff systems, and BDF, which is slower but robust to stiffness.

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

The particular function you're trying to integrate is non-stiff, so LSODA will use Adams on every iteration. You can check this by returning the infodict (...,full_output=True) and checking infodict['mused'].

由于 Adams 的方法不使用雅可比矩阵，因此您的梯度函数永远不会被调用.但是，如果您给 odeint 一个刚性函数进行积分，例如 Van der Pol 方程:

Since Adams' method does not use the Jacobian, your gradient function never gets called. However if you give odeint a stiff function to integrate, such as the 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 函数.

you should see that odeint switches to using BDF, and the Jacobian function now gets called.

如果您想更好地控制求解器，您应该查看 scipy.integrate.ode，这是一个更灵活的面向多个不同集成器的面向对象接口.

If you want more control over the solver, you should look into scipy.integrate.ode, which is a much more flexible object-oriented interface to multiple different integrators.

这篇关于为什么在 SciPy 中使用集成.odeint 时不调用 Dfun(gradient)?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！