刚性求解器

Question

我需要一个类似于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之间使用具有自适应步长的tf？

或者，在使用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'解算器比使用method='bdf'的'vode'要好得多。您可以获得初始化tvals = []所使用的步骤数的上限，在func中，可以得到do tvals.append(t)。当解决程序完成后，设置tvals = np.unique(tvals)。tvals的长度告诉您计算函数的时间值的数量。这并不完全是您想要的，但它确实显示了使用'lsoda'求解器和使用方法'bdf'的'vode'求解器之间的巨大差异。'lsoda'求解器使用的步骤数与您在评论中为matlab引用的顺序相同。(我用的是mu=10000，tf = 10。)

更新:如果您提供一个函数来计算雅可比矩阵，那么至少对于一个棘手的问题，它会对'vode'求解器产生巨大的影响。

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

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))