常微分方程
常微分方程 (ODE) 包含与一个自变量 t(通常称为时间)相关的因变量 y 的一个或多个导数。此处用于表示 y 相对于 t 的导数的表示法对于一阶导数为 y′,对于二阶导数为 y′′,依此类推。ODE 的阶数等于 y 在方程中出现的最高阶导数。
例如,这是一个二阶 ODE:
y′′=9y
在初始值问题中,从初始状态开始解算 ODE。利用初始条件 y 0 y_0 y0以及要在其中求得答案的时间段 ( t 0 , t f ) (t_0, t_f) (t0,tf),以迭代方式获取解。在每一步,求解器都对之前各步的结果应用一个特定算法。在第一个这样的时间步,初始条件将提供继续积分所需的必要信息。最终结果是,ODE 求解器返回一个时间步向量 t = [ t 0 , t 1 , t 2 , . . . , t f ] t=[t_0,t_1, t_2, ..., t_f] t=[t0,t1,t2,...,tf] 以及在每一步对应的解 t = [ y 0 , y 1 , y 2 , . . . , y f ] t=[y_0,y_1, y_2, ..., y_f] t=[y0,y1,y2,...,yf]。
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 ) y′=M^{-1}(t,y)f(t,y) 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 个方程,
(
y
1
′
y
2
′
⋮
y
n
′
)
=
(
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
)
)
\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}
⎝⎜⎜⎜⎛y1′y2′⋮yn′⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛f1(t,y1,y2,...,yn)f2(t,y1,y2,...,yn)⋮fn(t,y1,y2,...,yn)⎠⎟⎟⎟⎞
则用于编写该方程组代码的函数将返回一个向量,其中包含 n 个元素,对应于
y
1
′
,
y
2
′
,
.
.
.
,
y
n
′
y'_1,y'_2,...,y'_n
y1′,y2′,...,yn′ 值。例如,考虑以下包含两个方程的方程组
P
(
x
∣
p
a
x
)
=
{
y
′
1
=
y
2
y
2
′
=
y
1
y
2
−
2
P(x|pa_x)= \left \{\begin{array}{cc} y′_1=y_2\\ y'_2=y_1y_2-2 \end{array}\right.
P(x∣pax)={y′1=y2y2′=y1y2−2
用于编写该方程组代码的函数为
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
′
′
⋮
y
n
′
=
y
n
−
1
.
y_1=y \\ y_2=y' \\ y_2=y'' \\ \vdots \\ y'_n=y^{n-1}.
y1=yy2=y′y2=y′′⋮yn′=yn−1.
这些代换将生成一个包含 n 个一阶方程的方程组
{
y
1
′
=
y
2
y
2
′
=
y
3
⋮
y
n
′
=
f
(
t
,
y
1
,
y
−
2
,
.
.
.
,
y
n
)
.
\begin{cases} y'_1=y_2 \\ y'_2=y3 \\ \vdots \\ y'_n=f(t,y_1,y-2,...,y_n). \end{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧y1′=y2y2′=y3⋮yn′=f(t,y1,y−2,...,yn).
例如,考虑三阶 ODE
y ′ ′ ′ − y ′ ′ y + 1 = 0. y'''−y''y+1=0. y′′′−y′′y+1=0.
使用代换法
y 1 = y y 2 = y ′ y 3 = y ′ ′ . y_1 = y \\ y_2 = y' \\ y_3=y''. \\ y1=yy2=y′y3=y′′.
生成等效的一阶方程组
{ y 1 ′ = y 2 y 2 ′ = y 3 y ′ ′ ′ = y 1 y 3 − 1. \begin{cases} y'_1=y_2 \\ y'_2=y3 \\ y'''=y_1y_3-1. \end{cases} ⎩⎪⎨⎪⎧y1′=y2y2′=y3y′′′=y1y3−1.
此方程组的代码则为
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 y=y_1+iy_2 y=y1+iy2。为解算该方程,需要将实部和虚部分解为不同的解分量,最后重新组合相应的结果。从概念上讲,这类似于
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。 |