刚性 ODE 求解器

Question

我需要一个 ODE 求解器来解决类似于 MATLAB ode15s 的棘手问题。

对于我的问题，我需要检查不同初始值需要多少步（计算），并将其与我自己的 ODE 求解器进行比较。

我尝试过使用

solver = scipy.integrate.ode(f)
solver.set_integrator('vode', method='bdf', order=15, nsteps=3000)
solver.set_initial_value(u0, t0)

然后与：

i = 0
while solver.successful() and solver.t<tf:
    solver.integrate(tf, step=True)
    i += 1
print(i)

tf 是我的时间间隔的结束。

使用的函数定义为：

def func(self, t, u):
    u1 = u[1]
    u2 = mu * (1-numpy.dot(u[0], u[0]))*u[1] - u[0]
    return numpy.array([u1, u2])

初始值u0 = [ 2, 0] 是一个棘手的问题。

这意味着步数不应该取决于我的常量mu。

但确实如此。

我认为odeint-方法可以解决这个棘手的问题-但是我必须发送整个t-vector，因此需要设置已完成的步数，这会破坏重点我的任务。

是否可以在两个t0 和tf 之间使用具有自适应步长的odeint？

或者你能看到我在使用vode-integrator 时遗漏了什么吗？

Answer 1

我看到了类似的东西；使用'vode' 求解器，在'adams' 和'bdf' 之间更改方法不会改变步数。 （顺便说一句，使用order=15 是没有意义的；'vode' 求解器的'bdf' 方法的最大阶数是 5（而 'adams' 求解器的最大阶数是 12）。如果你离开参数，默认情况下应该使用最大值。）

odeint 是 LSODA 的封装。 ode 还提供了 LSODA 的封装： 将 'vode' 更改为 'lsoda'。不幸的是，'lsoda' 求解器忽略了 integrate 方法的 step=True 参数。

'lsoda' 求解器比'vode' 和method='bdf' 做得很多。 你可以得到一个上限 通过初始化 tvals = [] 使用的步数， 在func 中，执行tvals.append(t)。求解器完成后，设置 tvals = np.unique(tvals)。 tvals 的长度告诉你 评估您的函数的时间值的数量。 这并不完全是您想要的，但它确实显示出 巨大 的差异 在使用 'lsoda' 求解器和 'vode' 求解器之间 方法'bdf'。 'lsoda' 求解器使用的步数为 与您在评论中为 matlab 引用的顺序相同。 （我用过mu=10000，tf = 10。）

更新：事实证明，如果你提供一个计算雅可比矩阵的函数，至少对于一个僵硬的问题，它对 'vode' 求解器有很大的不同。

下面的脚本使用这两种方法运行 'vode' 求解器，它 运行'lsoda' 求解器。在每种情况下，它都会运行带有和不带有雅可比函数的求解器。这是它生成的输出：

vode   adams    jac=None  len(tvals) = 517992
vode   adams    jac=jac   len(tvals) = 195
vode   bdf      jac=None  len(tvals) = 516284
vode   bdf      jac=jac   len(tvals) = 55
lsoda           jac=None  len(tvals) = 49
lsoda           jac=jac   len(tvals) = 49

脚本：

from __future__ import print_function

import numpy as np
from scipy.integrate import ode


def func(t, u, mu):
    tvals.append(t)
    u1 = u[1]
    u2 = mu*(1 - u[0]*u[0])*u[1] - u[0]
    return np.array([u1, u2])


def jac(t, u, mu):
    j = np.empty((2, 2))
    j[0, 0] = 0.0
    j[0, 1] = 1.0
    j[1, 0] = -mu*2*u[0]*u[1] - 1
    j[1, 1] = mu*(1 - u[0]*u[0])
    return j


mu = 10000.0
u0 = [2, 0]
t0 = 0.0
tf = 10

for name, kwargs in [('vode', dict(method='adams')),
                     ('vode', dict(method='bdf')),
                     ('lsoda', {})]:
    for j in [None, jac]:
        solver = ode(func, jac=j)
        solver.set_integrator(name, atol=1e-8, rtol=1e-6, **kwargs)
        solver.set_f_params(mu)
        solver.set_jac_params(mu)
        solver.set_initial_value(u0, t0)

        tvals = []
        i = 0
        while solver.successful() and solver.t < tf:
            solver.integrate(tf, step=True)
            i += 1

        print("%-6s %-8s jac=%-5s " %
              (name, kwargs.get('method', ''), j.func_name if j else None),
              end='')

        tvals = np.unique(tvals)
        print("len(tvals) =", len(tvals))