主要内容

ODE事件位置

什么是事件位置?

求解某些ode系统的部分困难在于确定适当的停止解决方案的时间。积分区间中的最后时间可以由特定事件定义，而不是由数字定义。苹果从树上掉下来就是一个例子。ODE求解器应该在苹果落地后停止，但您可能不知道该事件将在何时发生。类似地，有些问题涉及到不终止解决方案的事件。一个例子是绕行星运行的月球。在这种情况下，你可能不想过早地停止整合，但你仍然想在月球围绕较大的天体完成一个周期的时候进行探测。

在ODE的解决过程中，使用事件函数检测某些事件何时发生。事件函数接受指定的表达式，并在该表达式等于零时检测事件。它们还可以在检测到事件时通知ODE求解器停止集成。

编写事件函数

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

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

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

输出参数价值，isterminal,方向是向量我Th元素对应于我th事件:

  • 值(我)是一个数学表达式描述我事件。当发生事件时值(我)等于0。

  • 终点(i) = 1如果集成要终止时我事件发生。否则，就是0．

  • 方向i = 0如果要定位所有的零(默认值)。值为＋1仅定位事件函数增加的零点，和－1仅定位事件函数递减的零点。指定方向= []的默认值0对于所有事件。

再考虑一个苹果从树上掉下来的例子。表示落体的ODE是

y ＇ ＇ ＝ − 1 + y ＇ 2 ，

在初始条件下 y （ 0 ） ＝ 1 而且 y ＇ （ 0 ） ＝ 0 ．您可以使用事件函数来确定时间 y （ t ） ＝ 0 ，也就是苹果落地的时候。对于这个问题，一个事件函数可以检测苹果何时触地

函数[位置，终端，方向]= appleEventsFcn(t,y)位置= y(1);我们希望为零的值终端= 1;%停止集成方向= 0;0可以从任意一个方向接近结束

事件信息

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

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

求解器返回的另外三个输出对应于检测到的事件:

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

  • 叶中每个事件时间的解决方案值te．

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

或者，您可以调用具有单个输出的求解器，如

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

在这种情况下，事件信息存储在结构as中sol.te，sol.ye,sol.ie．

限制

ODE求解器与事件函数一起使用的根查找机制有以下限制:

  • 如果在集成的第一步中发生了终结事件，那么求解器将该事件注册为非终结事件并继续集成。

  • 如果在第一步中发生了多个终端事件，那么只有第一个事件注册，求解器继续进行集成。

  • 零是通过步骤之间的符号交叉来确定的。因此，步骤之间交叉次数为偶数的函数的零位可能会被遗漏。

如果解算器步骤过去的事件，尝试减少RelTol而且AbsTol提高准确性。另外,组MaxStep在步长上设置一个上界。调整tspan不改变求解器所采取的步骤。

简单的活动地点:一个弹跳的球

打开脚本

这个例子展示了如何编写一个简单的事件函数来使用ODE求解器。示例文件ballode模拟弹跳球的运动。每次球反弹时，事件函数都会停止集成，然后以新的初始条件重新启动集成。随着球的反弹，集成会停止和重新启动几次。

弹跳球的方程是

数组$ $ \开始{}{cl} y ' _1 & # 38; = y_2 \ \ y ' _2 & # 38; = -9.8。\{数组}$ $

当球的高度发生反弹时，球就会反弹美元美元y_1 (t)减少后等于零。编码此行为的事件函数是

函数[value, terminal,direction] = bounceEvents(t,y) value = y(1);%检测高度= 0终端= 1;%停止集成方向= -1;%仅负方向

类型ballode运行该示例并说明如何使用事件函数来模拟球的弹跳。

ballode

高级事件地点:限制性三体问题

打开脚本

这个例子展示了如何使用事件函数的定向组件。示例文件orbitode模拟一个受限制的三体问题，其中一个物体围绕两个大得多的物体旋转。事件函数决定轨道上轨道体距离最近和最远的点。由于事件函数的值在轨道的最近点和最远点是相同的，所以零点交叉的方向是区分它们的地方。

限制性三体问题的方程是

数组$ $ \开始{}{cl} y ' _1 & # 38; = y_3 \ \ y ' _2 & # 38; = y y_4 \ \ ' _3 & # 38; = 2 y_4 + y_1 & # xA; \压裂{\μ^ * (y_1 + \μ)}{r_1 ^ 3} - \压裂{\μ(y_1 - \μ^ *}{r_2 ^ 3} \ \ y ' _4 & # 38; = & # xA; 2 y_3 + y_2 - \压裂{\μ^ * y_2} {r_1 ^ 3} - \压裂{\μy_2} {r_2 ^ 3}, & # xA; \{数组}$ $

在哪里

$ ${数组}{cl} \ \开始μ& # 38;= 1/82.45 \ \ \μ^ * & # 38;μ= 1 - \ \ \ r_1 & # 38; = & # xA; \√6 {(y_1 + \μ)^ 2 + y_2 ^ 2} \ \ r_2 & # 38; = \√6 {(y_1 \μ^ *)^ 2 + & # xA; y_2 ^ 2} \{数组}$ $结束。

解的前两个分量是质量无穷小的物体的坐标，所以把其中一个和另一个画出来就得到了物体的轨道。

事件功能嵌套在orbitode.m搜索两个事件。一个事件定位到距离起点最远的点，另一个定位到飞船返回起点的点。尽管积分器使用的步长不是由事件的位置决定的，但事件被准确地定位。在本例中，指定零点交叉方向的能力至关重要。返回起点的点和距离起点最大的点都具有相同的事件值，并使用交叉方向来区分它们。编码此行为的事件函数是

函数[value, terminal,direction] = orbitEvents(t,y)% dDSQdt是当前距离方程的导数。当地的当此值为零时，% minimum/maximum发生。dDSQdt = 2 * ((y(1:2)-y0(1:2))' * y(3:4));value = [dDSQdt;dDSQdt];终端= [1;0);%停止在本地最小值方向= [1;1);%[局部最小值，局部最大值]结束

类型orbitode运行示例。

orbitode
这是一个事件定位的例子，在这个例子中，指定过零点方向的能力是至关重要的。返回起始点的点和距离最大的点都有相同的事件函数值，用交叉方向来区分。调用活动事件函数的ODE45…请注意，积分器使用的步长不是由事件的位置决定的，并且事件仍然被准确地定位。

另请参阅

|

相关的话题