matlab使用教程(25)—常微分方程(ODE)选项

1.ODE 选项摘要

解算 ODE 经常要求微调参数、调整误差容限或向求解器传递附加信息。本主题说明如何指定选项以及每个选项与哪些微分方程求解器兼容。

1.1 选项语法

使用 odeset 函数创建 options 结构体，然后将其作为第四个输入参数传递给求解器。例如，要调整相对和绝对误差容限，请使用以下命令：

opts = odeset('RelTol',1e-2,'AbsTol',1e-5);
[t,y] = ode45(@odefun,tspan,y0,opts);

如果您使用不带输入参数的 odeset 命令，MATLAB 将显示每个选项的可能值列表，并以花括号 {} 表示默认值。

odeget 函数可以查询现有结构体中的选项值，您可以用该函数来根据条件动态更改选项值。例如，以下代码检测 Stats 是否设置为 'on' ，并根据需要更改其值：

if isempty(odeget(opts,'Stats'))
odeset(opts,'Stats','on')
end

1.2 选项与每个求解器的兼容性

在 odeset 中，有些选项是通用的，可与任何求解器兼容，还有一些选项是特定于求解器的。下表总结了每个选项与不同求解器之间的兼容性。

* 只有在解算不涉及质量矩阵的问题时，才对 ode15s 、 ode23t 和 ode23tb 求解器使用 NonNegative 参数。

** ode15i 的事件函数必须接受 yp 的第三个输入参数。

1.3 用法示例

MATLAB 包含几个示例文件，说明了如何使用各种选项。例如，键入 edit ballode 查看使用 'Events' 指定事件函数的示例，或键入 edit batonode 查看使用 'Mass' 指定质量矩阵的示例。有关示例文件及其使用的选项的完整总结，请参阅 “ODE 示例和文件摘要” 。

2.ODE 事件位置

2.1 什么是事件位置？

某些 ODE 方程组的解算难度部分在于确定停止求解的合适时间。积分区间中的最终时间可能由具体事件（而非数字）进行定义。从树上落下的苹果便是一个示例。ODE 求解器应当在苹果落地后立即停止，但您事前并不知道该事件会在何时发生。与之相似地，有些问题所涉及的事件则不会终止求解。沿行星轨道运行的卫星便是一个示例。这种情况下，您可能不希望提前停止积分，但仍希望检测到卫星每次围绕较大天体完成一个运行周期的时刻。

在 ODE 的求解过程中，使用事件函数来检测发生特定事件的时刻。事件函数采用您指定的表达式，并在该表达式等于零时检测事件。它们还能在检测到事件时提示 ODE 求解器停止积分。

2.2 编写事件函数

使用 odeset 函数的 'Events' 选项指定事件函数。事件函数必须具有一般形式

[value,isterminal,direction] = myEventsFcn(t,y)

在 ode15i 的情况下，事件函数还必须接受 yp 的第三个输入参数。

输出参数 value 、 isterminal 和 direction 均为向量，其第 i 个元素与第 i 个事件相对应：

• value(i) 是描述第 i 个事件的数学表达式。当 value(i) 等于零时发生事件。

• 如果当第 i 个事件发生时停止积分，则 isterminal(i) = 1 。否则为 0 。

• 如果需要查找全零值（默认值），则 direction(i) = 0 。值为 +1 时仅在事件函数递增的位置查找零，值为 -1 时仅在事件函数递减的位置查找零。指定 direction = [] 将对所有事件使用默认值 0。再次考虑苹果从树上落下的情形。表示落体的 ODE 为

初始条件为 y 0 = 1 和 y ′ 0 = 0 。您可以使用事件函数来确定 y t = 0 的时刻，即苹果落地的时刻。对于此问题，检测苹果何时落地的事件函数为

function [position,isterminal,direction] = appleEventsFcn(t,y)
position = y(1); % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end

2.3 事件信息

如果指定了事件函数，则使用三个额外的输出参数调用 ODE 求解器，如下所示

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

求解器返回的三个附加输出对应于检测到的事件：

• te 是发生事件的时刻的列向量。

• ye 包含 te 中的每个事件时刻对应的解值。

• ie 包含事件函数返回的向量的索引。这些值指示求解器检测到的事件。

您也可以使用单个输出来调用求解器，如下所示

sol = odeXY(odefun,tspan,y0,options)

这种情况下，事件信息以 sol.te 、 sol.ye 和 sol.ie 的形式存储在结构体中。

2.4 局限性

ODE 求解器采用的与事件函数配合使用的求根机制存在下列局限性：

• 如果在积分的第一步即发生终止事件，则求解器会将该事件记录为非终止事件并继续积分。

• 如果在第一步发生多个终止事件，则仅记录第一个事件，并且求解器会继续积分。

• 零值由每步之间的符号交叉确定。因此，对于两步间有偶数个交叉的函数而言，可能会错失其零值。如果求解器步长跨越了多个事件，请尝试减小 RelTol 和 AbsTol 以提高精度。也可以设置 MaxStep 以便为步长设置上限。调整 tspan 不会更改求解器所用的步长。

2.5 简单事件位置：弹球

此示例说明如何编写一个与 ODE 求解器配合使用的简单事件函数。示例文件 ballode 将模拟弹球的运动。事件函数在球每次弹起时停止积分，然后使用新的初始条件重新开始积分。在球的弹跳过程中，积分多次停止并重新开始。弹球的方程为

function [value,isterminal,direction] = bounceEvents(t,y)
value = y(1); % Detect height = 0
isterminal = 1; % Stop the integration
direction = -1; % Negative direction only

键入 ballode 以运行该示例并演示使用事件函数模拟球弹跳的过程。

ballode

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

此示例说明如何使用事件函数的定向分量。示例文件 orbitode 模拟限制性三体问题，其中一个主体环绕两个大得多的主体做轨道运行。事件函数将确定轨道中距离环绕主体最近和最远的点。由于事件函数在轨道最近点和最远点的值相同，因此将使用过零的方向来区分二者。限制性三体问题的方程为

前两个解分量是微小物体的坐标，因此针对一个分量绘制另一个分量可以得到物体的轨迹。orbitode.m 中嵌套的事件函数将搜索两个事件。一个事件查找距离起点最远的点，另一个事件查找宇宙飞船返回到起点的点。即使积分器使用的步长并非通过事件位置确定，也会准确定位事件。在此示例中，指定过零方向的功能非常重要。返回到起点的点和距离起点最远的点具有相同的事件值，并由交叉方向来区分这两个点。为此行为编码的事件函数为

function [value,isterminal,direction] = orbitEvents(t,y)
% dDSQdt is the derivative of the equation for current distance. Local
% minimum/maximum occurs when this value is zero.
dDSQdt = 2 * ((y(1:2)-y0(1:2))' * y(3:4));
value = [dDSQdt; dDSQdt];
isterminal = [1; 0]; % stop at local minimum
direction = [1; -1]; % [local minimum, local maximum]
end

键入 orbitode 以运行该示例。