刚性ODE求解器 [英] Stiff ODE-solver

问题描述

对于类似于MATLAB ode15s的刚性问题，我需要一个ODE求解器.

I need an ODE-solver for a stiff problem similar to MATLAB ode15s.

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

For my problem I need to check how many steps (calculations) is needed for different initial values and compare this to my own ODE-solver.

我尝试使用

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是我的时间间隔的结尾.

Where tf is the end of my time interval.

使用的功能定义为:

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]的问题很棘手.

Which with the initial value u0 = [ 2, 0] is a stiff problem.

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

This means that the number of steps should not depend on my constant mu.

但确实如此.

我认为odeint方法可以解决这个棘手的问题-但随后我必须发送整个t -vector，因此需要设置要执行的步骤数量，这会破坏观点我的任务.

I think the odeint-method can solve this as a stiff problem - but then I have to send in the whole t-vector and therefore need to set the amount of steps that is done and this ruins the point of my assignment.

反正有没有在两个t0和tf之间使用具有自适应步长的odeint?

Is there anyway to use odeint with adaptive stepsize between two t0 and tf?

或者您能看到使用vode -integrator时我错过的任何东西吗?

Or can you see anything I miss in the use of the vode-integrator?

推荐答案

我看到了类似的东西；使用'vode'求解器，在'adams'和'bdf'之间更改方法不会显着改变步骤数. (顺便说一句，使用order=15毫无意义； 'vode'求解器的'bdf'方法的最大顺序为5(而'adams'求解器的最大顺序为12).如果离开参数，默认情况下应该使用最大值.)

I'm seeing something similar; with the 'vode' solver, changing methods between 'adams' and 'bdf' doesn't change the number of steps by very much. (By the way, there is no point in using order=15; the maximum order of the 'bdf' method of the 'vode' solver is 5 (and the maximum order of the 'adams' solver is 12). If you leave the argument out, it should use the maximum by default.)

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

odeint is a wrapper of LSODA. ode also provides a wrapper of LSODA: change 'vode' to 'lsoda'. Unfortunately the 'lsoda' solver ignores the step=True argument of the integrate method.

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

The 'lsoda' solver does much better than 'vode' with method='bdf'. You can get an upper bound on the number of steps that were used by initializing tvals = [], and in func, do tvals.append(t). When the solver completes, set tvals = np.unique(tvals). The length of tvals tells you the number of time values at which your function was evaluated. This is not exactly what you want, but it does show a huge difference between using the 'lsoda' solver and the 'vode' solver with method 'bdf'. The number of steps used by the 'lsoda' solver is on the same order as you quoted for matlab in your comment. (I used mu=10000, tf = 10.)

更新:事实证明，至少对于一个刚性问题，如果提供计算雅可比矩阵的函数，则对'vode'求解器会产生巨大的影响.

Update: It turns out that, at least for a stiff problem, it make a huge difference for the 'vode' solver if you provide a function to compute the Jacobian matrix.

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

The script below runs the 'vode' solver with both methods, and it runs the 'lsoda' solver. In each case, it runs the solver with and without the Jacobian function. Here's the output it generates:

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

这篇关于刚性ODE求解器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！