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′=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 个方程，
$\begin{pmatrix} y'_1 \\ y'_2 \\ \vdots \\ y'_n \end{pmatrix} = \begin{pmatrix} f_1(t,y_1,y_2,...,y_n) \\ f_2(t,y_1,y_2,...,y_n) \\ \vdots \\ f_n(t,y_1,y_2,...,y_n) \end{pmatrix}$
则用于编写该方程组代码的函数将返回一个向量，其中包含 n 个元素，对应于 $y'_1,y'_2,...,y'_n$ 值。例如，考虑以下包含两个方程的方程组
$P(x|pa_x)= \left \{\begin{array}{cc} y′_1=y_2\\ y'_2=y_1y_2-2 \end{array}\right.$
用于编写该方程组代码的函数为

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

讯享网

高阶 ODE

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

$y^{'''} - y^{''} y + 1 = 0 .$

使用代换法

讯享网

$y_1 = y \\ y_2 = y' \\ y_3=y''. \\$

生成等效的一阶方程组

$\begin{cases} y'_1=y_2 \\ y'_2=y3 \\ y'''=y_1y_3-1. \end{cases}$

此方程组的代码则为

讯享网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+iy_2$ 。为解算该方程，需要将实部和虚部分解为不同的解分量，最后重新组合相应的结果。从概念上讲，这类似于

例如，如果 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。

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