ODE 事件位置

本主题介绍如何在使用求解器函数（ode45、ode15s 等）求解 ODE 时检测事件。

什么是事件位置？

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

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

编写事件函数

使用 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'' = - 1 + y'^{2},$

初始条件为 $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

事件信息

如果指定了事件函数，则使用三个额外的输出参量调用 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 的形式存储在结构体中。

局限性

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

如果在积分的第一步即发生终止事件，则求解器会将该事件记录为非终止事件并继续积分。
如果在第一步发生多个终止事件，则仅记录第一个事件，并且求解器会继续积分。
零值由每步之间的符号交叉确定。因此，对于两步间有偶数个交叉的函数而言，可能会错失其零值。

如果求解器步长跨越了多个事件，请尝试减小 RelTol 和 AbsTol 以提高精度。也可以设置 MaxStep 以便为步长设置上限。调整 tspan 不会更改求解器所用的步长。

简单事件位置：弹球

打开脚本

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

弹球的方程为

$\begin{array}{cl} y'_1 &= y_2\\ y'_2 &= -9.8 . \end{array}$

当球的高度在逐渐减小后等于零时，发生球弹起。为此行为编码的事件函数为

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

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

打开脚本

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

限制性三体问题的方程为

$\begin{array}{cl} y'_1 &= y_3\\ y'_2 &= y_4 \\ y'_3 &= 2y_4 + y_1 -
\frac{\mu^* (y_1 + \mu)}{r_1^3} - \frac{\mu(y_1 - \mu^*}{r_2^3}\\ y'_4 &=
-2y_3 + y_2 - \frac{\mu^* y_2}{r_1^3} - \frac{\mu y_2}{r_2^3},
\end{array}$

其中

$\begin{array}{cl} \mu &= 1/82.45\\ \mu^* &= 1-\mu\\ r_1 &=
\sqrt{(y_1+\mu)^2+y_2^2}\\ r_2 &= \sqrt{(y_1-\mu^*)^2 +
y_2^2}.\end{array}$

前两个解分量是微小物体的坐标，因此针对一个分量绘制另一个分量可以得到物体的轨迹。

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 以运行该示例。

orbitode

This is an example of event location where the ability to
specify the direction of the zero crossing is critical.  Both
the point of return to the initial point and the point of
maximum distance have the same event function value, and the
direction of the crossing is used to distinguish them.

Calling ODE45 with event functions active...

Note that the step sizes used by the integrator are NOT
determined by the location of the events, and the events are
still located accurately.