Dfun (Jacobian)在odeint Python中的使用

Question

这是一个相当普遍的问题，其中一部分可能适用于任何耦合常微分方程的数值模拟，一部分可能仅适用于scipy.integrate Python库中的odeint方法。

首先，odeint如何使用手动输入的雅可比矩阵( Dfunc参数)，为什么它会使大型ODE系统的速度如此之快？

其次，更确切地说，对于我的具体问题，如果雅可比函数略有错误，odeint会产生不正确的解，还是只会减慢它的速度？通过肉眼(模拟结果的动画)我看不出有什么不同；我希望这是因为雅可比是正确的，但我不能完全确定。

Answer 1

odeint调用的ODEPACK方法使用隐式多步法来求解y'=f(t,y)。这意味着它必须在每一步中求解一个非线性方程组，本质上是不动点方程。

y[i+1] = h*b[0]*f(t[i+1],y[i+1]) + C

其中C是之前的点y[i-k]和f(t[i-k],y[i-k])，k=0,..,q的值的线性组合，这是计算y[i+1]的当前步骤中的常量。

现在在任何不动点方程中，如果它在感兴趣的区域收缩，你可以迭代它，找到一个不动点y=g(y)作为极限。收敛将是线性的，其中收缩因子主要由雅可比g'(yfix)的范数/谱半径决定。现在假设g的分裂已知为线性部分和非线性余数(可能具有小的线性部分)，

g(y) = A*y + r(y) = yfix + A*(y-yfix) + [r(y)-r(yfix)]

其中“小”与最后分解中的最后一项有关。将不动点方程中的线性部分转移到另一端，从而得到一个新的不动点方程

y = (I-A)^(-1) * r(y)

所获得的数据。它的收敛因子现在取决于到固定点的距离，越接近，收敛越快。这可能有一个非零的下界，并且不是二次的，但它将比原始迭代更快。