Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

ode23tb

求解刚性微分方程 - 梯形法则 + 后向差分公式

说明

示例

[t,y] = ode23tb(odefun,tspan,y0)(其中 tspan = [t0 tf])求微分方程组 y'=f(t,y)t0tf 的积分,初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。

所有 MATLAB® ODE 求解器都可以解算 y'=f(t,y) 形式的方程组,或涉及质量矩阵 M(t,y)y'=f(t,y) 的问题。求解器都使用类似的语法。ode23s 求解器只能解算质量矩阵为常量的问题。ode15sode23t 可以解算具有奇异质量矩阵的问题,称为微分代数方程 (DAE)。使用 odesetMass 选项指定质量矩阵。

示例

[t,y] = ode23tb(odefun,tspan,y0,options) 还使用由 options(使用 odeset 函数创建的参数)定义的积分设置。例如,使用 AbsTolRelTol 选项指定绝对误差容限和相对误差容限,或者使用 Mass 选项提供质量矩阵。

[t,y,te,ye,ie] = ode23tb(odefun,tspan,y0,options) 还求 (t,y) 的函数(称为事件函数)在何处为零。在输出中,te 是事件的时间,ye 是事件发生时的解,ie 是触发的事件的索引。

对于每个事件函数,应指定积分是否在零点处终止以及过零方向是否重要。为此,请将 'Events' 属性设置为函数(例如 myEventFcn@myEventFcn),并创建一个对应的函数:[value,isterminal,direction] = myEventFcn(t,y)。有关详细信息,请参阅 ODE 事件位置

sol = ode23tb(___) 返回一个结构体,您可以将该结构体与 deval 结合使用来计算区间 [t0 tf] 中任意点位置的解。您可以使用上述语法中的任何输入参数组合。

示例

全部折叠

在对求解器的调用中,可将只有一个解分量的简单 ODE 指定为匿名函数。该匿名函数必须同时接受两个输入 (t,y),即使其中一个输入未使用也是如此。

解算 ODE

y=-10t.

使用时间区间 [0,2] 和初始条件 y0 = 1

tspan = [0 2];
y0 = 1;
[t,y] = ode23tb(@(t,y) -10*t, tspan, y0);

对解绘图。

plot(t,y,'-o')

Figure contains an axes. The axes contains an object of type line.

表示张弛振荡的 van der Pol 方程组是一个刚性方程组示例。范围环的某些区域中的解分量变化缓慢且问题呈刚性,这些区域与变化显著的非刚性区域交替出现。

方程组为:

$$\begin{array}{cl} y_1' &= y_2\\y_2' &= 1000(1-y_1^2)y_2-y_1\end{array}$$

初始条件为 $y_1(0)=2$$y_2(0)=0$。函数 vdp1000 随 MATLAB® 一起提供,用于对方程进行编码。

function dydt = vdp1000(t,y)
%VDP1000  Evaluate the van der Pol ODEs for mu = 1000.
%
%   See also ODE15S, ODE23S, ODE23T, ODE23TB.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

使用 ode45 以及默认的相对和绝对误差容限(分别为 1e-31e-6)解算此方程组特别慢,需要好几分钟才能完成解算和绘制。ode45 需要几百万个时间步才能完成积分,因为它要在刚度区域满足容差要求。

下面是通过 ode45 获得的解的绘图,需要很长时间来计算。请注意,通过刚度区域需要大量时间步。

使用 ode23tb 求解器对刚性方程组求解,然后绘制解 y 的第一列对时间点 t 的图。ode23tb 求解器通过刚性区域所需的步数远远少于 ode45

[t,y] = ode23tb(@vdp1000,[0 3000],[2 0]);
plot(t,y(:,1),'-o')

ode23s 仅适用于使用两个输入参数(ty)的函数。但是,通过在函数外部定义参数并在指定函数句柄时传递这些参数,可以传入额外参数。

解算 ODE

$$y'' = \frac{A}{B} t y.$$

将该方程重写为一阶方程组可以得到

$$\begin{array}{cl} y'_1 &= y_2\\ y'_2 &= \frac{A}{B} t y_1.
\end{array}$$

odefcn.m 将此方程组表示为接受四个输入参数(tyAB)的函数。

function dydt = odefcn(t,y,A,B)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = (A/B)*t.*y(1);

使用 ode23tb 解算 ODE。指定函数句柄,使其将 AB 的预定义值传递给 odefcn

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode23tb(@(t,y) odefcn(t,y,A,B), tspan, y0);

绘制结果。

plot(t,y(:,1),'-o',t,y(:,2),'-.')

ode15s 求解器是适用于大多数刚性问题的首选求解器。但是,对于特定类型的问题,其他刚性求解器可能更高效。本示例使用所有四个刚性 ODE 求解器解算刚性测试方程。

请考虑以下测试方程:

y=-λy.

随着 λ 的量级增加,方程的刚性逐渐增强。使用 λ=1×109、初始条件 y(0)=1 和时间区间 [0 0.5]。这些值使该问题足够刚性,从而使 ode45ode23 需要对该方程进行积分。此外,还使用 odeset 传入常量 Jacobian J=fy=-λ 并打开求解器统计信息的显示。

lambda = 1e9;
y0 = 1;
tspan = [0 0.5];
opts = odeset('Jacobian',-lambda,'Stats','on');

使用 ode15sode23sode23tode23tb 对方程求解。生成子图进行比较。

subplot(2,2,1)
tic, ode15s(@(t,y) -lambda*y, tspan, y0, opts), toc
104 successful steps
1 failed attempts
212 function evaluations
0 partial derivatives
21 LU decompositions
210 solutions of linear systems
Elapsed time is 1.785332 seconds.
title('ode15s')
subplot(2,2,2)
tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc
63 successful steps
0 failed attempts
191 function evaluations
0 partial derivatives
63 LU decompositions
189 solutions of linear systems
Elapsed time is 0.495668 seconds.
title('ode23s')
subplot(2,2,3)
tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc
95 successful steps
0 failed attempts
125 function evaluations
0 partial derivatives
28 LU decompositions
123 solutions of linear systems
Elapsed time is 0.658186 seconds.
title('ode23t')
subplot(2,2,4)
tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc
71 successful steps
0 failed attempts
167 function evaluations
0 partial derivatives
23 LU decompositions
236 solutions of linear systems
Elapsed time is 0.777485 seconds.
title('ode23tb')

Figure contains 4 axes. Axes 1 with title ode15s contains 2 objects of type line. Axes 2 with title ode23s contains 2 objects of type line. Axes 3 with title ode23t contains 2 objects of type line. Axes 4 with title ode23tb contains 2 objects of type line.

这些刚性求解器都表现不错,但对于这个特定问题,ode23s 完成积分所用的步长最少,运行速度最快。由于指定了常量 Jacobian,任何求解器都不需要计算偏导数即可计算出解。指定 Jacobian 对 ode23s 最有利,因为它通常需要在每个步长中计算 Jacobian。

对于一般的刚性问题,刚性求解器的性能因问题的形式和指定的选项而异。提供 Jacobian 矩阵或稀疏模式始终可以提高求解器解决刚性问题的效率。但是,由于刚性求解器使用 Jacobian 的方式不同,计算效率的提高程度也有很大差异。实际上,如果方程组非常大,或者要解算很多次,则非常值得研究一下不同求解器的性能,以最大程度地缩短执行时间。

输入参数

全部折叠

要求解的函数,指定为指向待积分函数的句柄。

对于标量 t 和列向量 y 来说,函数 dydt = odefun(t,y) 必须返回数据类型为 singledouble 的列向量 dydt,该列向量对应于 f(t,y)odefun 必须同时接受输入参数 ty,即使其中一个参数未在函数中使用也是如此。

例如,要解算 y'=5y3,请使用此函数:

function dydt = odefun(t,y)
dydt = 5*y-3;

对于方程组,odefun 的输出为向量。向量中的每个元素是一个方程的解。例如,要求解

y'1=y1+2y2y'2=3y1+2y2

使用函数:

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);

有关如何为函数 odefun 提供其他参数的信息,请参阅参数化函数

示例: @myFcn

数据类型: function_handle

积分区间,指定为向量。tspan 必须至少是一个二元素向量 [t0 tf],用于指定初始时间和最终时间。要获取 t0tf 之间的特定时间的解,请使用 [t0,t1,t2,...,tf] 形式的长向量。tspan 中的元素必须单调递增或单调递减。

求解器在初始时间 tspan(1) 施加由 y0 给出的初始条件,然后求 tspan(1)tspan(end) 的积分:

  • 如果 tspan 有两个元素,[t0 tf],求解器将返回在该区间内的每个内部积分步计算的解。

  • 如果 tspan 包含两个以上的元素,[t0,t1,t2,...,tf],求解器将返回在给定点处计算的解。但是,求解器不会精确步进到 tspan 中指定的每个点。此时,求解器使用自己的内部积分步来计算解,然后在 tspan 中请求的各点处计算解。在指定点处生成的解与在每个内部积分步计算的解具有相同的准确度级别。

    指定多个中间点对计算效率影响甚微,但在大型系统中可能会影响内存管理。

求解器使用 tspan 的值计算 InitialStepMaxStep 的合适值:

  • 如果 tspan 包含多个中间点 [t0,t1,t2,...,tf],则指定的点表示了问题的规模,这可能影响求解器使用的 InitialStep 的值。因此,根据您是将 tspan 指定为二元素向量还是包含中间点的向量,求解器获得的解可能有所不同。

  • tspan 中的初始值和最终值用于计算最大步长 MaxStep。因此,更改 tspan 中的初始值或最终值可能导致求解器使用不同步长序列,从而可能会更改解。

示例: [1 10]

示例: [1 3 5 7 9 10]

数据类型: single | double

初始条件,指定为向量。y0 的长度必须与 odefun 的向量输出相同,使 y0odefun 中定义的每个方程包含一个初始条件。

数据类型: single | double

options 结构体,指定为结构体数组。使用 odeset 函数创建或修改 options 结构体。有关与每个求解器兼容的选项列表,请参阅 ODE 选项摘要

示例: options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) 指定 1e-5 的相对误差容限、打开求解器统计信息的显示,并指定输出函数 @odeplot 在计算时绘制解。

数据类型: struct

输出参数

全部折叠

计算点,以列向量形式返回。

  • 如果 tspan 包含两个元素,[t0 tf],则 t 包含用于执行积分的内部计算点。

  • 如果 tspan 包含两个以上元素,则 ttspan 相同。

解,以数组形式返回。y 中的每一行都与 t 的相应行中返回的值处的解相对应。

事件的时间,以列向量形式返回。te 中的事件时间对应于 ye 中返回的解,而 ie 指定发生了哪个事件。

事件时间的解,以数组形式返回。te 中的事件时间对应于 ye 中返回的解,而 ie 指定发生了哪个事件。

触发的事件函数的索引,以列向量形式返回。te 中的事件时间对应于 ye 中返回的解,而 ie 指定发生了哪个事件。

用于计算的结构体,以结构体数组形式返回。此结构体与 deval 函数一起使用,用于计算区间 [t0 tf] 内任何点的解。sol 结构体数组始终包括下列字段:

结构体字段说明

sol.x

求解器选择的步的行向量。

sol.y

解。每列 sol.y(:,i) 包含时间 sol.x(i) 处的解。

sol.solver

求解器名称。

此外,如果指定了 Events 选项并且检测到事件,则 sol 还包括下列字段:

结构体字段说明

sol.xe

事件发生的点。sol.xe(end) 包含终止事件(如果有)的确切点。

sol.ye

sol.xe 中的事件相对应的解。

sol.ie

Events 选项中指定的函数所返回的向量的索引。这些值指示求解器检测到的事件。

算法

ode23tb 是 TR-BDF2 的实现,这是一个隐式 Runge-Kutta 公式,该公式第一阶段计算采用梯形法则步,第二阶段计算采用二阶后向差分公式。从构造上而言,两级计算使用相同的迭代矩阵。与 ode23sode23t 一样,此求解器在解算具有宽松容差的问题时,可能效率高于 ode15s [1], [2]

参考

[1] Bank, R. E., W. C. Coughran, Jr., W. Fichtner, E. Grosse, D. Rose, and R. Smith, “Transient Simulation of Silicon Devices and Circuits,” IEEE Trans. CAD, 4 (1985), pp. 436–451.

[2] Shampine, L. F. and M. E. Hosea, “Analysis and Implementation of TR-BDF2,” Applied Numerical Mathematics 20, 1996.

在 R2006a 之前推出