python - 刚性 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-method 可以解决这个棘手的问题 - 但我必须发送整个 t-vector，因此需要设置数量已完成的步骤，这破坏了我的任务重点。

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

或者您能看出我在使用 vode 集成器时遗漏了什么吗？

Answer 1

我看到了类似的东西；使用 'vode' 求解器，在 'adams' 和 'bdf' 之间更改方法不会对步数产生太大影响。 (顺便说一句，使用 order=15 没有意义；'vode' 的 'bdf' 方法的最大阶数solver 是 5('adams' solver 的最大阶数是 12)。如果你省略参数，它应该默认使用最大值。)

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

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

更新:事实证明，至少对于一个棘手的问题，如果您提供一个函数来计算雅可比矩阵。

下面的脚本使用这两种方法运行 '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))