link管理

链接快照平台

输入网页链接，自动生成快照
标签化管理网页链接

ode23s

求解刚性微分方程 - 低阶方法

语法


            
             [t,y] = ode23s(odefun,tspan,y0)


            
             [t,y] = ode23s(odefun,tspan,y0,options)


            
             [t,y,te,ye,ie] = ode23s(odefun,tspan,y0,options)


            
             sol = ode23s(
             
              ___
             
             )

说明

示例

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

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

示例

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

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

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

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

示例

全部折叠

具有一个解分量的 ODE

打开实时脚本

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

解算 ODE

$y^{'} = - 10 t .$

指定时间区间 [0 2] 和初始条件 y0 = 1 。

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

对解绘图。

plot(t,y,'-o')

解算刚性 ODE

打开脚本

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

方程组为：

%VDP1000 Evaluate the van der Pol ODEs for mu = 1000. % See also ODE15S, ODE23S, ODE23T, ODE23TB. % Jacek Kierzenka and Lawrence F. Shampine dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

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

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

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

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

向 ODE 函数传递额外的参数

打开脚本

ode23s 仅适用于使用两个输入参量（ t 和 y ）的函数。但是，通过在函数外部定义参数并在指定函数句柄时传递这些参数，可以传入额外参数。

解算 ODE

y0 = [0 0.01]; [t,y] = ode23s(@(t,y) odefcn(t,y,A,B), tspan, y0);

绘制结果。

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

比较刚性 ODE 求解器

打开实时脚本

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

请考虑以下测试方程：

$y^{'} = - λ y .$

随着 $λ$ 的量级增加，方程的刚性逐渐增强。使用 $λ = 1 \times 1 0^{9}$ 、初始条件 $y (0) = 1$ 和时间区间 [0 0.5] 。这些值使该问题足够刚性，从而使 ode45 和 ode23 需要对该方程进行积分。此外，还使用 odeset 传入常量雅可比 $J = \frac{\partial f}{\partial y} = - λ$ 并打开求解器统计信息的显示。

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

使用 ode15s 、 ode23s 、 ode23t 和 ode23tb 对方程求解。生成子图进行比较。

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.117982 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.179014 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.233392 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.234603 seconds.

title('ode23tb')

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

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

输入参数

全部折叠

`odefun` — 要求解的函数
函数句柄

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

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

例如，要解算 $y' = 5 y - 3$ ，请使用此函数：

function dydt = odefun(t,y)
dydt = 5*y-3;
对于方程组，odefun 的输出为向量。向量中的每个元素均为一个方程右侧的计算值。例如，考虑以下包含两个方程的方程组
 $\begin{array}{l} y'_{1} = y_{1} + 2 y_{2} \\ y'_{2} = 3 y_{1} + 2 y_{2} \end{array}$ 
计算每个时间步下每个方程右侧值的函数为
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` — 积分区间
向量

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

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

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

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

求解器使用 tspan 的值计算 InitialStep 和 MaxStep 的合适值：

如果 tspan 包含多个中间点 [t0,t1,t2,...,tf] ，则指定的点表示了问题的规模，这可能影响求解器使用的 InitialStep 的值。因此，根据您是将 tspan 指定为二元素向量还是包含中间点的向量，求解器获得的解可能有所不同。
tspan 中的初始值和最终值用于计算最大步长 MaxStep 。因此，更改 tspan 中的初始值或最终值可能导致求解器使用不同步长序列，从而可能会更改解。

示例: [1 10]

示例: [1 3 5 7 9 10]

数据类型: single | double

`y0` — 初始条件
向量

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

数据类型: single | double

`options` — options 结构体
结构体数组

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

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

数据类型: struct

输出参量

全部折叠

`t` — 求值点
列向量

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

如果 tspan 包含两个元素 [t0 tf] ，则 t 包含用于执行积分的内部计算点。
如果 tspan 包含两个以上元素，则 t 与 tspan 相同。

`y` — 解
数组

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

`te` — 事件的时间
列向量

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

`ye` — 事件时间处的解
数组

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

`ie` — 触发的事件函数的索引
列向量

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

`sol` — 用于计算的结构体
结构体数组

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

结构体字段	描述
`sol.x`	求解器选择的步的行向量。
`sol.y`	解。每列 `sol.y(:,i)` 包含时间 `sol.x(i)` 处的解。
`sol.solver`	求解器名称。

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

结构体字段	描述
`sol.xe`	事件发生的点。 `sol.xe(end)` 包含终止事件（如果有）的确切点。
`sol.ye`	与 `sol.xe` 中的事件相对应的解。
`sol.ie`	`Events` 选项中指定的函数所返回的向量的索引。这些值指示求解器检测到的事件。

算法

ode23s 基于经过改进的二阶罗森布罗克公式。由于它是一个单步求解器，所以在解算允许宽松容差的问题或者解变化很快的问题时，效率可能高于 ode15s 。它可以解算一些刚性问题，而使用 ode15s 解算这些问题的效率非常低下。 ode23s 求解器在积分的每一步都计算雅可比，所以提供雅可比矩阵对它的可靠性和效率至关重要 [1] 。

参考

[1] Shampine, L. F. and M. W. Reichelt, “ The MATLAB ODE Suite ,” SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.

版本历史记录

在 R2006a 之前推出

另请参阅

ode15s | odeset | odeget | deval

ode23s

语法

说明

示例

具有一个解分量的 ODE

解算刚性 ODE

向 ODE 函数传递额外的参数

比较刚性 ODE 求解器

输入参数

odefun — 要求解的函数 函数句柄

tspan — 积分区间 向量

y0 — 初始条件 向量

options — options 结构体 结构体数组

输出参量

t — 求值点 列向量

y — 解 数组

te — 事件的时间 列向量

ye — 事件时间处的解 数组

ie — 触发的事件函数的索引 列向量

sol — 用于计算的结构体 结构体数组

算法

参考

版本历史记录

另请参阅

主题

`odefun` — 要求解的函数
函数句柄

`tspan` — 积分区间
向量

`y0` — 初始条件
向量

`options` — options 结构体
结构体数组

`t` — 求值点
列向量

`y` — 解
数组

`te` — 事件的时间
列向量

`ye` — 事件时间处的解
数组

`ie` — 触发的事件函数的索引
列向量

`sol` — 用于计算的结构体
结构体数组