选择 ode 求解器 -凯发k8网页登录

选择 ode 求解器

常微分方程

常微分方程 (ode) 包含与一个自变量 t（通常称为时间）相关的因变量 y 的一个或多个导数。此处用于表示 y 关于 t 的导数的表示法对于一阶导数为 $y'$ ，对于二阶导数为 $y''$ ，依此类推。ode 的阶数等于 y 在方程中出现的最高阶导数。

例如，这是一个二阶 ode：

$y'' = 9 y$

在初始值问题中，从初始状态开始解算 ode。利用初始条件 $y_{0}$ 以及要在其中求得答案的时间段 $(t_{0}, t_{f})$ ，以迭代方式获取解。在每一步，求解器都对之前各步的结果应用一个特定算法。在第一个这样的时间步，初始条件将提供继续积分所需的必要信息。最终结果是，ode 求解器返回一个时间步向量 $t = [t_{0}, t_{1}, t_{2}, ..., t_{f}]$ 以及在每一步对应的解 $y = [y_{0}, y_{1}, y_{2}, ..., y_{f}]$ 。

ode 的类型

matlab^® 中的 ode 求解器可解算以下类型的一阶 ode：

$y' = f (t, y)$ 形式的显式 ode。
$m (t, y) y' = f (t, y)$ 形式的线性隐式 ode，其中 $m (t, y)$ 为非奇异质量矩阵。该质量矩阵可以是时间或状态依赖的矩阵，也可以是常量矩阵。线性隐式 ode 涉及在质量矩阵中编码的一阶 y 导数的线性组合。
线性隐式 ode 可随时变换为显式形式 $y' = m^{- 1} (t, y) f (t, y)$ 。不过，将质量矩阵直接指定给 ode 求解器可避免这种既不方便还可能带来大量计算开销的变换操作。
如果 $y'$ 的某些分量缺失，则这些方程称为微分代数方程或 dae，并且 dae 方程组会包含一些代数变量。代数变量是导数未出现在方程中的因变量。可通过对方程求导来将 dae 方程组重写为等效的一阶 ode 方程组，以消除代数变量。将 dae 重写为 ode 所需的求导次数称为微分指数。ode15s 和 ode23t 求解器可解算微分指数为 1 的 dae。
$f (t, y, y') = 0$ 形式的完全隐式 ode。完全隐式 ode 不能重写为显式形式，还可能包含一些代数变量。ode15i 求解器专为完全隐式问题（包括微分指数为 1 的 dae）而设计。

可通过使用函数创建 options 结构体，来针对某些类型的问题为求解器提供附加信息。

ode 方程组

您可以指定需要解算的任意数量的 ode 耦合方程，原则上，方程的数量仅受计算机可用内存的限制。如果方程组包含 n 个方程，

$(\begin{matrix} y'_{1} \\ y'_{2} \\ ⋮ \\ y'_{n} \end{matrix}) = (\begin{matrix} f_{1} (t, y_{1}, y_{2}, ..., y_{n}) \\ f_{2} (t, y_{1}, y_{2}, ..., y_{n}) \\ ⋮ \\ f_{n} (t, y_{1}, y_{2}, ..., y_{n}) \end{matrix}),$

则用于编写该方程组代码的函数将返回一个向量，其中包含 n 个元素，对应于 $y'_{1}, y'_{2}, \dots, y'_{n}$ 值。例如，考虑以下包含两个方程的方程组

${\begin{cases} y'_{1} = y_{2} \\ y'_{2} = y_{1} y_{2} - 2 . \end{cases}$

用于编写该方程组代码的函数为

function dy = myode(t,y)
  dy(1) = y(2);
  dy(2) = y(1)*y(2)-2;
end

高阶 ode

matlab ode 求解器仅可解算一阶方程。您必须使用常规代换法，将高阶 ode 重写为等效的一阶方程组

$\begin{array}{l} y_{1} = y \\ y_{2} = y' \\ y_{3} = y'' \\ ⋮ \\ y_{n} = y^{(n - 1)} . \end{array}$

这些代换将生成一个包含 n 个一阶方程的方程组

${\begin{cases} y'_{1} = y_{2} \\ y'_{2} = y_{3} \\ ⋮ \\ y'_{n} = f (t, y_{1}, y_{2}, ..., y_{n}) . \end{cases}$

例如，考虑三阶 ode

$y''' - y'' y 1 = 0.$

使用代换法

$\begin{array}{l} y_{1} = y \\ y_{2} = y' \\ y_{3} = y'' \end{array}$

生成等效的一阶方程组

${\begin{cases} y'_{1} = y_{2} \\ y'_{2} = y_{3} \\ y'_{3} = y_{1} y_{3} - 1. \end{cases}$

此方程组的代码则为

function dydt = f(t,y)
  dydt(1) = y(2);
  dydt(2) = y(3);
  dydt(3) = y(1)*y(3)-1;
end

复数 ode

考虑复数 ode 方程

$y' = f (t, y),$

其中 $y = y_{1} i y_{2}$ 。为解算该方程，需要将实部和虚部分解为不同的解分量，最后重新组合相应的结果。从概念上讲，这类似于

$\begin{array}{l} y v = [real (y) imag (y)] \\ f v = [real (f (t, y)) imag (f (t, y))] . \end{array}$

例如，如果 ode 为 $y' = y t 2 i$ ，则可以使用函数文件来表示该方程：

function f = complexf(t,y)
  f = y.*t   2*i;
end

然后，分解实部和虚部的代码为

function fv = imaginaryode(t,yv)
  % construct y from the real and imaginary components
  y = yv(1)   i*yv(2);            
  % evaluate the function
  yp = complexf(t,y);             
  % return real and imaginary in separate components
  fv = [real(yp); imag(yp)]; 
end

在运行求解器以获取解时，初始条件 y0 也会分解为实部和虚部，以提供每个解分量的初始条件。

y0 = 1 i;
yv0 = [real(y0); imag(y0)];
tspan = [0 2];
[t,yv] = ode45(@imaginaryode, tspan, yv0);

获得解后，将实部和虚部分量组合到一起可获得最终结果。

y = yv(:,1)   i*yv(:,2);

基本求解器选择

ode45 适用于大多数 ode 问题，一般情况下应作为您的首选求解器。但对于精度要求更宽松或更严格的问题而言，ode23、ode78、ode89 和 ode113 可能比 ode45 更加高效。

一些 ode 问题具有较高的计算刚度或难度。术语“刚度”无法精确定义，但一般而言，当问题的某个位置存在标度差异时，就会出现刚度。例如，如果 ode 包含的两个解分量在时间标度上差异极大，则该方程可能是刚性方程。如果非刚性求解器（例如 ode45）无法解算某个问题或解算速度极慢，则可以将该问题视为刚性问题。如果您观察到非刚性求解器的速度很慢，请尝试改用 ode15s 等刚性求解器。在使用刚性求解器时，可以通过提供 jacobian 矩阵或其稀疏模式来提高可靠性和效率。

下表提供了关于何时使用每种不同求解器的一般指导原则。

问题类型	精度	何时使用
非刚性	中	大多数情况下，您应当首先尝试求解器 `ode45`。
	低	对于容差较宽松的问题或在刚度适中的情况下，`ode23` 可能比 `ode45` 更加高效。
	低到高	对于具有严格误差容限的问题或在 ode 函数需要大量计算开销的情况下，`ode113` 可能比 `ode45` 更加高效。
	高	对于具有高准确度要求的平滑解的问题，`ode78` 可能比 `ode45` 更加高效。
	高	对于非常平滑的问题，当在较长的时间区间内进行积分时，或当容差特别严格时，`ode89` 可能比 `ode78` 更加高效。
刚性	低到中	若 `ode45` 失败或效率低下并且您怀疑面临刚性问题，请尝试 `ode15s`。此外，当解算微分代数方程 (dae) 时，请使用 `ode15s`。
	低	对于误差容限较宽松的问题，`ode23s` 可能比 `ode15s` 更加高效。它可以解算一些刚性问题，而使用 `ode15s` 解算这些问题的效率不高。 `ode23s` 会在每一步计算 jacobian，因此通过 `odeset` 提供 jacobian 有利于最大限度地提高效率和精度。如果存在质量矩阵，则它必须为常量矩阵。
	低	对于仅仅是刚度适中的问题，并且您需要没有数值阻尼的解，请使用 `ode23t`。 `ode23t` 可解算微分代数方程 (dae)。
	低	与 `ode23s` 一样，对于误差容限较宽松的问题，`ode23tb` 求解器可能比 `ode15s` 更加高效。
完全隐式	低	对于完全隐式问题 f(t,y,y’) = 0 和微分指数为 1 的微分代数方程 (dae)，请使用 `ode15i`。

有关何时使用每种求解器的详细信息和更多建议，请参阅 [5]。

ode 示例和文件摘要

有几个示例文件可用作大多数 ode 问题的有用起点。要运行微分方程示例应用，以便轻松浏览和运行示例，请键入

odeexamples

要打开单独的示例文件进行编辑，请键入

edit examplefilename.m

要运行示例，请键入

examplefilename

此表包含可用的 ode 和 dae 示例文件及其使用的求解器和选项的列表。其中包含示例子集的链接，这些示例也已直接发布在文档中。

示例文件	指定的选项	说明	文档链接
`amp1dae`	`'mass'`	刚性 dae - 包含常量奇异质量矩阵的电路
`ballode`	`'events'` `'outputfcn'` `'outputsel'` `'refine'` `'initialstep'` `'maxstep'`	简单事件位置 - 弹球
`batonode`	`'mass'`	时间依赖和状态依赖的质量矩阵的 ode - 短棒的移动
`brussode`	`'jpattern'` `'vectorized'`	刚性大问题 - 化学反应中的扩散 (brusselator)
`burgersode`	`'mass'` `'mstatedependence'` `'jpattern'` `'mvpattern'` `'reltol'` `'abstol'`	强烈依赖于状态的质量矩阵的 ode - 使用移动网格方法解算 burgers 方程
`fem1ode`	`'mass'` `'mstatedependence'` `'jacobian'`	时间依赖的质量矩阵的刚性问题 - 有限元方法	-
`fem2ode`	`'mass'`	常量质量矩阵的刚性问题 - 有限元方法	-
`hb1ode`	-	在非常长的区间内解算的刚性 ode 问题 - robertson 化学反应	-
`hb1dae`	`'mass'` `'reltol'` `'abstol'` `'vectorized'`	基于守恒定律的刚性线性隐式 dae - robertson 化学反应
`ihb1dae`	`'reltol'` `'abstol'` `'jacobian'`	刚性完全隐式 dae - robertson 化学反应
`iburgersode`	`'reltol'` `'abstol'` `'jacobian'` `'jpattern'`	隐式 ode 方程组 - burgers 方程	-
`kneeode`	`'nonnegative'`	具有非负约束的“膝盖问题”
`orbitode`	`'reltol'` `'abstol'` `'events'` `'outputfcn'`	高级事件位置 - 限制性三体问题
`rigidode`	-	非刚性问题 - 刚体在不受外力作用时的欧拉方程
`vdpode`	`'jacobian'`	可参数化的 van der pol 方程（对于较大 μ，为刚性问题）

参考

[1] shampine, l. f. and m. k. gordon, computer solution of ordinary differential equations: the initial value problem, w. h. freeman, san francisco, 1975.

[2] forsythe, g., m. malcolm, and c. moler, computer methods for mathematical computations, prentice-hall, new jersey, 1977.

[3] kahaner, d., c. moler, and s. nash, numerical methods and software, prentice-hall, new jersey, 1989.

[4] shampine, l. f., numerical solution of ordinary differential equations, chapman & hall, new york, 1994.

[5] shampine, l. f. and m. w. reichelt, “the matlab ode suite,” siam journal on scientific computing, vol. 18, 1997, pp. 1–22.

[6] shampine, l. f., gladwell, i. and s. thompson, solving odes with matlab, cambridge university press, cambridge uk, 2003.

选择 ode 求解器 -凯发k8网页登录