MATLAB 数学应用 微分方程 常微分方程 选择ODE求解器

常微分方程

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

例如，这是一个二阶 ODE：

y′′=9y

在初始值问题中，从初始状态开始解算 ODE。利用初始条件 $y_0$以及要在其中求得答案的时间段 $(t_0,t_f)$，以迭代方式获取解。在每一步，求解器都对之前各步的结果应用一个特定算法。在第一个这样的时间步，初始条件将提供继续积分所需的必要信息。最终结果是，ODE 求解器返回一个时间步向量 $t=[t_0,t_1,t_2,...,t_f]$ 以及在每一步对应的解 $t=[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\prime =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）而设计。

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

ODE 方程组

您可以指定需要解算的任意数量的 ODE 耦合方程，原则上，方程的数量仅受计算机可用内存的限制。如果方程组包含 n 个方程，
$$\left(\begin{array}{c}{y}_1^{\prime }\\ y_2^{\prime }\\ ⋮\\ y_n^{\prime }\end{array}\right)=\left(\begin{array}{c}{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{array}\right)$$
则用于编写该方程组代码的函数将返回一个向量，其中包含 n 个元素，对应于 $y_1^{\prime },y_2^{\prime },...,y_n^{\prime }$ 值。例如，考虑以下包含两个方程的方程组
$$P(x|p{a}_x)=\{\begin{array}{c}y{\prime }_1=y_2\\ y_2^{\prime }=y_1{y}_2-2\end{array}$$
用于编写该方程组代码的函数为

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

高阶 ODE

MATLAB ODE 求解器仅可解算一阶方程。您必须使用常规代换法，将高阶 ODE 重写为等效的一阶方程组
$$\begin{gathered} y_1=y \\ {y}_2=y^{\prime } \\ {y}_2=y^{\prime \prime } \\ ⋮ \\ {y}_n^{\prime }=y^{n-1}. \end{gathered}$$
这些代换将生成一个包含 n 个一阶方程的方程组
$$\{\begin{array}{l}{y}_1^{\prime }=y_2\\ y_2^{\prime }=y3\\ ⋮\\ y_n^{\prime }=f(t,y_1,y-2,...,y_n).\end{array}$$
例如，考虑三阶 ODE

$y^{\prime \prime \prime }-y^{\prime \prime }y+1=0.$

使用代换法

$$\begin{gathered} y_1=y \\ {y}_2=y^{\prime } \\ {y}_3=y^{\prime \prime }. \end{gathered}$$

生成等效的一阶方程组

$$\{\begin{array}{l}{y}_1^{\prime }=y_2\\ y_2^{\prime }=y3\\ y^{\prime \prime \prime }=y_1{y}_3-1.\end{array}$$

此方程组的代码则为

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

复数 ODE

考虑复数 ODE 方程

y′=f(t,y) ,

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

yv=[Real(y)    Imag(y)]
fv=[Real(f(t,y))    Imag(f(t,y))] .

例如，如果 ODE 为 y′=yt+2i，则可以使用函数文件来表示该方程。

function f = complexf(t,y)
f = y.t + 2i;
然后，分解实部和虚部的代码为

function fv = imaginaryODE(t,yv)
y = yv(1) + i*yv(2);            

yp = complexf(t,y);             

fv = [real(yp); imag(yp)];      

在运行求解器以获取解时，初始条件 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 和 ode113 可能比 ode45 更加高效。

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

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

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