Main Content

选择 ODE 求解器

常微分方程

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

例如,这是一个二阶 ODE:

y''=9y

初始值问题中,从初始状态开始解算 ODE。利用初始条件 y0 以及要在其中求得答案的时间段 (t0,tf),以迭代方式获取解。在每一步,求解器都对之前各步的结果应用一个特定算法。在第一个这样的时间步,初始条件将提供继续积分所需的必要信息。最终结果是,ODE 求解器返回一个时间步向量 t=[t0,t1,t2,...,tf] 以及在每一步对应的解 y=[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'=M1(t,y)f(t,y)。不过,将质量矩阵直接指定给 ODE 求解器可避免这种既不方便还可能带来大量计算开销的变换操作。

  • 如果 y' 的某些分量缺失,则这些方程称为微分代数方程或 DAE,并且 DAE 方程组会包含一些代数变量。代数变量是导数未出现在方程中的因变量。可通过对方程求导来将 DAE 方程组重写为等效的一阶 ODE 方程组,以消除代数变量。将 DAE 重写为 ODE 所需的求导次数称为微分指数。ode15sode23t 求解器可解算微分指数为 1 的 DAE。

  • f(t,y,y')=0 形式的完全隐式 ODE。完全隐式 ODE 不能重写为显式形式,还可能包含一些代数变量。ode15i 求解器专为完全隐式问题(包括微分指数为 1 的 DAE)而设计。

ODE 方程组

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

(y'1y'2y'n)=(f1(t,y1,y2,...,yn)f2(t,y1,y2,...,yn)fn(t,y1,y2,...,yn)),

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

{y'1=y2y'2=y1y22.

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

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

高阶 ODE

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

y1=yy2=y'y3=y''yn=y(n1).

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

{y'1=y2y'2=y3y'n=f(t,y1,y2,...,yn).

例如,考虑三阶 ODE

y'''y''y+1=0.

使用代换法

y1=yy2=y'y3=y''

生成等效的一阶方程组

{y'1=y2y'2=y3y'3=y1y31.

此方程组的代码则为

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=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 + 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 问题,一般情况下应作为您的首选求解器。但对于精度要求更宽松或更严格的问题而言,ode23ode78ode89ode113 可能比 ode45 更加高效。

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

您可以使用 ode 对象根据问题的属性自动选择求解器。如果您不确定要使用哪个求解器,请参考下表,其中提供关于每个求解器的适用情形的一般指导。

求解器问题类型精度何时使用
ode45非刚性

大多数情况下,您应当首先尝试求解器 ode45

ode23

对于容差较宽松的问题或在刚度适中的情况下,ode23 可能比 ode45 更加高效。

ode113低到高

对于具有严格误差容限的问题或在 ODE 函数需要大量计算开销的情况下,ode113 可能比 ode45 更加高效。

ode78

对于具有高准确度要求的平滑解的问题,ode78 可能比 ode45 更加高效。

ode89

对于非常平滑的问题,当在较长的时间区间内进行积分时,或当容差特别严格时,ode89 可能比 ode78 更加高效。

ode15s刚性低到中

ode45 失败或效率低下并且您怀疑面临刚性问题,请尝试 ode15s。此外,当解算微分代数方程 (DAE) 时,请使用 ode15s

ode23s

对于误差容限较宽松的问题,ode23s 可能比 ode15s 更加高效。它可以解算一些刚性问题,而使用 ode15s 解算这些问题的效率不高。

ode23s 会在每一步计算雅可比,因此通过 odeset 提供雅可比有利于最大限度地提高效率和精度。

如果存在质量矩阵,则它必须为常量矩阵。

ode23t

对于仅仅是刚度适中的问题,并且您需要没有数值阻尼的解,请使用 ode23t

ode23t 可解算微分代数方程 (DAE)。

ode23tb

ode23s 一样,对于误差容限较宽松的问题,ode23tb 求解器可能比 ode15s 更加高效。

ode15i完全隐式

对于完全隐式问题 f(t,y,y’) = 0 和微分指数为 1 的微分代数方程 (DAE),请使用 ode15i

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

ODE 示例和文件摘要

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

odeexamples

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

edit exampleFileName.m

要运行示例,请键入

exampleFileName

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

示例文件使用的求解器指定的选项描述文档链接
amp1daeode23t
  • 'Mass'

刚性 DAE - 包含常量奇异质量矩阵的电路

求解刚性晶体管微分代数方程
ballodeode23
  • 'Events'

  • 'OutputFcn'

  • 'OutputSel'

  • 'Refine'

  • 'InitialStep'

  • 'MaxStep'

简单事件位置 - 弹球

ODE 事件位置
batonodeode45
  • 'Mass'

时间依赖和状态依赖的质量矩阵的 ODE - 短棒的移动

求解抛向空中的短棒的运动方程
brussodeode15s
  • 'JPattern'

  • 'Vectorized'

刚性大问题 - 化学反应中的扩散 (Brusselator)

解算刚性 ODE
burgersodeode15s
  • 'Mass'

  • 'MStateDependence'

  • 'JPattern'

  • 'MvPattern'

  • 'RelTol'

  • 'AbsTol'

强烈依赖于状态的质量矩阵的 ODE - 使用移动网格方法解算伯格斯方程

求解具有强状态依赖质量矩阵的 ODE
fem1odeode15s
  • 'Mass'

  • 'MStateDependence'

  • 'Jacobian'

时间依赖的质量矩阵的刚性问题 - 有限元方法

fem2odeode23s
  • 'Mass'

常量质量矩阵的刚性问题 - 有限元方法

hb1odeode15s

在非常长的区间内解算的刚性 ODE 问题 - 罗伯逊化学反应

hb1daeode15s
  • 'Mass'

  • 'RelTol'

  • 'AbsTol'

  • 'Vectorized'

基于守恒定律的刚性线性隐式 DAE - 罗伯逊化学反应

将 Robertson 问题作为半显式微分代数方程 (DAE) 求解
ihb1daeode15i
  • 'RelTol'

  • 'AbsTol'

  • 'Jacobian'

刚性完全隐式 DAE - 罗伯逊化学反应

将 Robertson 问题作为隐式微分代数方程 (DAE) 求解
iburgersodeode15i
  • 'RelTol'

  • 'AbsTol'

  • 'Jacobian'

  • 'JPattern'

隐式 ODE 方程组 - 伯格斯方程

kneeodeode15s
  • 'NonNegative'

具有非负约束的“膝盖问题”

非负 ODE 解
orbitodeode45
  • 'RelTol'

  • 'AbsTol'

  • 'Events'

  • 'OutputFcn'

高级事件位置 - 限制性三体问题

ODE 事件位置
rigidodeode45

非刚性问题 - 刚体在不受外力作用时的欧拉方程

求解非刚性 ODE
vdpodeode15s
  • 'Jacobian'

可参数化的 van der Pol 方程(对于较大 μ,为刚性问题)

解算刚性 ODE

参考

[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.

另请参阅

| |

相关主题

外部网站