WOLFRAM

Wolfram 语言与系统参考资料中心
NDSolve 的 “EventLocator” 方法

NDSolve 的 “EventLocator” 方法

引言局限性
事件定位方法选项总结
实例

引言

能够探测和精确定位在一个微分系统中的变化往往是很有用的. 例如，随着一个奇点或者状态变化的检测，可以采取适当的行动，比如重新启动积分.

一个微分系统

的一个事件，是一个沿着解，实值事件函数为零的点：

也可以考虑布尔值事件函数，在这种情况下，当函数从 True 变为 False 时事件发生，反之亦然.

内置于 NDSolve 的 "EventLocator" 方法实际上相当于一种控制器方法；它处理对事件的检查并且采取适当的动作，但是除此以外，微分系统的积分则完全留给一个下层的基本方法.

本节将举例说明 "EventLocator" 方法和选项的基本用途. 随后的章节给出了更多的事件定位的具体应用，如周期检测、庞加莱截面和不连续点处理.

这些初始化命令加载了一些有用的程序包，这些程序包具有一些要求解的微分方程，并且定义了一些功用函数：
Wolfram Language code: Needs["DifferentialEquations`NDSolveProblems`"]; Needs["DifferentialEquations`NDSolveUtilities`"]; Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]; Needs["GUIKit`"];

一个简单的例子是定位一个事件，比如从一个非平衡位置开始的钟摆，何时通过其最低点，并将在该点停止积分.

以下对钟摆方程进行积分，积分到解 y[t] 穿过轴的第一个点：
Wolfram Language code: sol = NDSolve[{y''[t] + Sin[y[t]] == 0, y'[0] == 0, y[0] == 1}, y, {t, 0, 10}, Method -> {"EventLocator", "Event" -> y[t]}]

从该解，用户可以看到钟摆在大约 时到达最低点 . 使用 InterpolatingFunctionAnatomy 程序包，就可能从 InterpolatingFunction 对象提取值.

以下提取事件发生的点，并且画出解（黑色）及其导数（蓝色）至该点的图形：
Wolfram Language code: end = InterpolatingFunctionDomain[First[y /. sol]][[1, -1]]; Plot[Evaluate[{y[t], y'[t]} /. First[sol]], {t, 0, end}, PlotStyle -> {{StandardBrown}, {StandardBlue}}]

当用户使用事件定位器方法的时候，要定位的事件和当找到一个事件时要采取的行动都是通过 "EventLocator" 方法的方法选项指定的.

检测事件的默认动作是前面演示的那样停止积分. 事件动作可以是任意表达式. 每当探测到一个事件时，便使用数值值取代问题变量对该表达式进行计算.

以下显示对一个阻尼摆，每当一个事件 被探测到时的时间和值：
Wolfram Language code: NDSolve[{y''[t] + .1 y'[t] + Sin[y[t]] == 0, y'[0] == 0, y[0] == 1}, y, {t, 0, 10}, Method -> {"EventLocator", "Event" -> y'[t] - y[t], "EventAction" :> Print["y'[", t, "] = y[", t, "] = ", y[t]]}]

请注意，在这个例子中，"EventAction" 选项使用 RuleDelayed（）给出，用以防止计算进行，除非事件被定位.

从显示的输出用户可以看到，当动作没有停止积分时，一个实件的多个事例可以被检测到. 事件是当事件表达式的符号改变时检测到的. 用户可以使用 "Direction" 选项来把事件限制到仅在一个特定的方向的符号变化.

这里是对一个阻尼摆，收集速度从负值变为正值的点. Reap 和 Sow 是编程组成结构. 当用户事先不知道会有多少数据时，它们对于收集数据是很有用的. Reap[expr] 给出 expr 的值，以及在计算过程中，Sow 所应用的所有表达式. 这里 Reap 封装了 NDSolve 的使用，而 Sow 是事件动作的一部分，这允许用户对于一个事件的每个实例收集数据：
Wolfram Language code: Reap[NDSolve[{y''[t] + .1 y'[t] + Sin[y[t]] == .1 Cos[t], y'[0] == 0, y[0] == 1}, y, {t, 0, 50}, Method -> {"EventLocator", "Event" -> y'[t], "Direction" -> 1, "EventAction" :> Sow[{t, y[t], y'[t]}]}]]

从前面例子的输出，用户可能会注意到，事件是在导数近似为零时被探测到. 当该方法在底层积分器的一个步骤中探测一个事件的存在时（通过事件表达式符号的改变），那么它采用一个数值方法来近似地找到根的位置. 由于定位过程是数值的，用户应该只能得到近似的结果. 定位方法选项 AccuracyGoal、PrecisionGoal 和 MaxIterations 可以用予那些利用 FindRoot 控制容差来求根的定位方法.

对于布尔值事件函数，当函数从 True 切换到 False 时，一个事件发生，反之亦然. "Direction" 选项可以用来限制仅从 True 改变为 False 的事件（"Direction"->-1）或者仅从 False 改变为 True 的事件 （"Direction"->1）.

以下将打开一个带有一个按钮的小窗口，点击按钮时，变量 stop 的值从初始值 False 变为 True：
Wolfram Language code: NDSolve`stop = False; GUIRun[Widget["Panel", {Widget["Button", { "label" -> "Stop", BindEvent["action", Script[NDSolve`stop = True]]}]}]];
以下对钟摆方程积分，直到按钮被点击（或者系统耗尽了内存）：
Wolfram Language code: NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 1, y'[0] == 0}, y, {t, 0, ∞}, Method -> {"EventLocator", "Event" :> NDSolve`stop}, MaxSteps -> ∞]

注意，在这个例子中，"Event" 选项使用 RuleDelayed（:>）指定，用以避免计算 stop 的立即值，并且将其设置为函数.

用户可以指定多个事件. 如果事件函数将数值计算结果存入一个列表，那么该列表的每个组成部分被认为是一个单独的事. 通过把这些选项的值指定为具有适当长度的列表，用户可以对每个事件指定不同的动作、方向等等.

以下对钟摆方程积分，直到按钮被点击. 在积分过程中记录钟摆的完整摆动的数目：
Wolfram Language code: NDSolve`stop = False; swings = 0;{ NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 0, y'[0] == 1}, y, {t, 0, 1000000}, Method -> {"EventLocator", "Event" :> {y[t], NDSolve`stop}, "EventAction" :> {swings++, Throw[Null, "StopIntegration"]}, "Direction" -> {1, All}}, MaxSteps -> Infinity], swings}

正如从前面的例子中可以看到的，我们有可能混合实值和布尔值事件函数. 组成部分的期望数目和每个组成部分的类型都是基于初始条件的值，并且需要在整个积分过程中保持一致.

"EventLocator" 的 "EventCondition" 选项允许用户指定事件检测需要满足的额外的布尔条件. 在可能的情况下，用此而不是一个布尔事件是有利的，因为求根过程可以更有效地完成.

以下是一旦衰减使能量积分减少为 后，第一次满足 时，终止一个阻尼摆的积分：
Wolfram Language code: sol = NDSolve[{y''[t] + .1 y'[t] + Sin[y[t]] == 0, y'[0] == 1, y[0] == 0}, y, {t, 0, 100}, Method -> {"EventLocator", "Event" -> y[t], "EventCondition" -> (y'[t] ^ 2 / 2 - Cos[y[t]] < -0.9), "EventAction" :> Throw[end = t, "StopIntegration"]}]
以下生成解（棕色）、导数（蓝色）和能量积分（绿色）的图形. 能量临界值以红色显示：
Wolfram Language code: Plot[Evaluate[{y[t], y'[t], y'[t] ^ 2 / 2 - Cos[y[t]], -.9} /. First[sol]], {t, 0, end}, PlotStyle -> {{Black}, {Blue}, {Green}, {Red}}]

"EventLocator" 的 Method 选项允许指定在积分中使用的具体数值方法.

事件定位方法

"EventLocator" 方法的工作原理是通过采取一步低层方法、并且检查某一事件函数的符号（或者奇偶性）在步骤端点处是否不同来实现的. 事件函数应该是实值函数或者布尔值函数，所以，如果有改变，在步骤间隔中必须有一个事件. 对于每个在一个步进中有事件发生的事件函数，将执行一个细化的程序来定位事件在该区间中的位置.

有几种不同的方法可用于细化事件位置. 其中包括简单地采用积分区间开头和末尾的解、事件值的线性插值、以及使用括号求根方法. 应采用的适当方法取决于执行速度和定位精度之间的取舍.

如果该事件动作是用来停止积分，那么积分停止的特定值取决于从 "EventLocator" 的 "EventLocationMethod" 选项获得的值.

单一事件的定位通常足够快，使得所采用的方法不会显著影响整体计算时间. 然而，当一个事件被多次检测，定位的细化方法可以产生很大的影响.

“StepBegin” 和 “StepEnd” 方法

当事件定位的精确位置并不重要、或者不反映底层计算的任何精度信息时，最粗略的方法是适用的. 上一节的 stop 按钮的例子就是这样一个例子：时间步进的计算如此之快，以至于无法在一个特定的时间步进内对按钮的点击计时，更不用说在某一时间步进内的一个特定点. 因此，基于事件本身的准确程度，细化没有任何意义. 用户可以通过使用 "StepBegin" 或者 "StepEnd" 定位方法对此进行指定. 在事件的定义是试探性的或者不够确切的任何例子中，这会是一个合适的选择.

“LinearInterpolation” 方法

当事件结果用于图形绘制中的点的时候，用户只需要把该事件定位到图形的分辨率上. 虽然仅使用步进端点过于草率，基于事件函数值的单一线性插值就足够了.

如下表示数值积分的连续网点上的事件函数值：

线性插值给出：

那么，事件时间的线性近似是：

线性插值也可以用来近似事件时间上的解. 但是，由于导数值 和 在网点是可用的，事件上的解的更好的近似可以廉价地使用三次厄米（Hermite）插值计算：

其中，下面是适当定义的插值权重：

用户可以指定基于具有设置 "LinearInterpolation" 的单一线性插值的细化方法.

以下计算钟摆方程单个周期的解，并且绘制出该周期的解：
Wolfram Language code: sol = First[NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0}, y, {t, 0, ∞}, Method -> {"EventLocator", "Event" -> y'[t], "EventAction" :> Throw[end = t, "StopIntegration"], "Direction" -> -1, "EventLocationMethod" -> "LinearInterpolation", Method -> "ExplicitRungeKutta"}]]; Plot[Evaluate[{y[t], y'[t]} /. sol], {t, 0, end}, PlotStyle -> {{StandardBrown}, {StandardBlue}}]

在整个周期内的图像的分辨率上，用户看不出端点可能不是确切位于导数与轴相交的位置. 但是，如果充分放大的话，用户就可以看到误差.

以下显示接近末端点的图形：
Wolfram Language code: Plot[Evaluate[y'[t] /. sol], {t, end * (1 - .001), end}, PlotStyle -> StandardBlue]

线性插值方法对于大部分图示目的是足够的，例如接下来的章节中所示的庞加莱截面的例子. 请注意，对于布尔值事件函数，线性插值实际上只是一个二分步骤，所以线性插值方法可能对于图形是不足的.

布伦特 （Brent） 方法

默认的定位方法是事件定位方法 "Brent"，即 使用带有布伦特方法的FindRoot 寻找事件的位置. 布伦特方法始于一个带括号的根，并且把基于插值法和二分法的步骤结合起来，确保收敛速度至少与二分法一样好. 用户可以使用 "Brent" 事件定位方法的方法选项来控制 FindRoot 试图获取的事件函数的根的准确度和精度. 默认情况下是采用与 NDSolve 用于局部误差控制相同的准确度和精度来求得根.

对于支持连续或者稠密输出的方法，事件函数的变量可以通过使用连续输出公式很有效很简单地找到. 不过，对于不支持连续输出的方法，解需要执行一步底层方法来计算，这样做的代价可能相对比较昂贵. 另外一种获取准确度不及方法的阶数、但与将 FindRoot 用于从 NDSolve 返回的 InterpolatingFunction 对象一致的近似解的方法是，使用三次厄米（Hermite）插值，通过基于在步骤结束时的解的值和解的导数值的插值方法获取步骤中间的解的近似值.

比较

这个例子使用数个不同的事件定位方法对钟摆方程进行积分，并且比较事件被发现的时间.

以下定义将要采用的事件定位方法：
Wolfram Language code: eventmethods = {"StepBegin", "StepEnd", "LinearInterpolation", Automatic};
以下对该系统进行积分，并且显示使用的方法，和当积分终止时的自变量的值：
Wolfram Language code: Map[ NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0}, y, {t, 0, ∞}, Method -> {"EventLocator", "Event" -> y'[t], "EventAction" :> Throw[Print[#, ": t = ", t, " y'[t] = ", y'[t]], "StopIntegration"], "Direction" -> -1, Method -> "ExplicitRungeKutta", "EventLocationMethod" -> #}]&, eventmethods ];

实例

自由落体

该系统对在遇到空气阻力的情况下、由于引力作用的自由落体进行建模（见 [M04]）.

事件行动存储落体撞击地面的时间，并且终止积分过程：
Wolfram Language code: sol = y[t] /. First[NDSolve[{y''[t] == -1 + y'[t] ^ 2, y[0] == 1, y'[0] == 0}, y, {t, 0, Infinity}, Method -> {"EventLocator", "Event" :> y[t], "EventAction" :> Throw[tend = t, "StopIntegration"]}]]
以下画出解的图形，并且通过画圈突出显示初始点和最终点（绿色和红色）：
Wolfram Language code: plt = Plot[sol, {t, 0, tend}, Frame -> True, Axes -> False, PlotStyle -> StandardBlue, DisplayFunction -> Identity]; grp = Graphics[{{Green, Circle[{0, 1}, 0.025]}, {StandardRed, Circle[{tend, sol /. t -> tend}, 0.025]}}]; Show[plt, grp, DisplayFunction -> $DisplayFunction]

范德波尔振荡器的周期

范德波尔振荡器是常微分方程的一个极端刚性系统的例子. 事件定位方法可以调用任何方法来对常微分方程系统进行实际积分. 默认方法，Automatic，在必要时自动切换至适用于刚性系统的方法，以便刚性不至引出什么问题.

以下对参数的一个特定值 ，对范德波尔振荡器进行积分，直至变量 到达初始值和方向为止：
Wolfram Language code: vsol = NDSolve[(| | | | ------------------------------------------- | --------- | | y1'[t] == y2[t] | y1[0] == 2 | | y2'[t] == 1000 (1 - y1[t] ^ 2) y2[t] - y1[t] | y2[0] == 0 |), {Subscript[y, 1], Subscript[y, 2]}, {t, 3000}, Method -> {"EventLocator", "Event" -> Subscript[y, 2][t], "Direction" -> -1}]

请注意，不考虑位于初始条件的事件.

通过选择 NDSolve 解的端点，就可以编写一个函数，以返回作为 的一个函数的周期.

以下定义一个函数，它返回作为 的一个函数的周期：
Wolfram Language code: vper[μ_] := Module[{vsol}, vsol = First[Subscript[y, 2] /. NDSolve[(| | | | ---------------------------------------- | --------- | | y1'[t] == y2[t] | y1[0] == 2 | | y2'[t] == μ (1 - y1[t] ^ 2) y2[t] - y1[t] | y2[0] == 0 |), {Subscript[y, 1], Subscript[y, 2]}, {t, Max[100, 3 μ]}, Method -> {"EventLocator", "Event" -> Subscript[y, 2][t], "Direction" -> -1}]]; InterpolatingFunctionDomain[vsol][[1, -1]]];
这里使用该函数来计算 时的周期：
Wolfram Language code: vper[1000]

当然，我们很容易把这个方法推广到任意具有周期性解的系统中.

庞加莱截面

使用庞加莱截面是可视化高维微分系统的解的一个有用的技术.

对于一个交互式图形界面，参阅程序包 EquationTrekker.

Hénon–Heiles系统

定义 Hénon–Heiles系统，该系统对星系中的恒星运动进行建模.

以下从 NDSolveProblems 程序包获取 Hénon–Heiles 系统：
Wolfram Language code: system = GetNDSolveProblem["HenonHeiles"]; vars = system["DependentVariables"];eqns = {system["System"], system["InitialConditions"]}

在这种情况下，我们感兴趣的庞加莱截面是当轨道通过 时 平面中点的集合.

由于数值积分的实际结果并不是必须的，有可能可以通过把输出变量列表（NDSolve 的第二个变量中）指定为空，或者 {} 来避免把所有数据存储在 InterpolatingFunction 中. 这意味着 NDSolve 不会产生InterpolatingFunction 作为输出，避免了存储大量不必要的数据. NDSolve 的确发出一条消息 NDSolve::noout 警告将不会有输出函数，但是在这种情况下，它可以安全地被关闭，因为我们所感兴趣的数据是从事件动作中收集到的.

这里使用了线性插值事件定位方法，因为这里计算的目的是在分辨率相对较低的图像中查看结果. 如果用户正在做一个需要放大图像以查看详细信息或者找到某个特征的例子，比如庞加莱 map 中的一个固定点，那么使用默认的定位方法就更加适合了.

这里将关闭关于没有输出的警告消息：
Wolfram Language code: Off[NDSolve::noout];
这里使用具有固定步长0.25的四阶显式朗格-库塔方法对 Hénon–Heiles系统进行积分.事件动作是在 和 的值上使用 Sow：
Wolfram Language code: data = Reap[ NDSolve[eqns, {}, {T, 10000}, Method -> {"EventLocator", "Event" -> Subscript[Y, 1][T], "EventAction" :> Sow[{Subscript[Y, 2][T], Subscript[Y, 4][T]}], "EventLocationMethod" -> "LinearInterpolation", "Method" -> {"FixedStep", "Method" -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4}}}, StartingStepSize -> 0.25, MaxSteps -> ∞]; ];
以下绘制庞加莱截面的图形. 所收集的数据在 Reap 的结果的最后一部分可以找到，而点列表在第一部分可以找到：
Wolfram Language code: psdata = data[[-1, 1]]; ListPlot[psdata, Axes -> False, Frame -> True, AspectRatio -> 1]

由于 Hénon–Heiles系统是哈密顿的，一个辛方法提供了这个例子的更好的定性结果.

以下使用具有固定步长0.25的四阶辛分区朗格-库塔方法对 Hénon–Heiles系统进行积分. 事件行动是对 和 上的值采用 Sow：
Wolfram Language code: sdata = Reap[ NDSolve[eqns, {}, {T, 10000}, Method -> {"EventLocator", "Event" -> Subscript[Y, 1][T], "EventAction" :> Sow[{Subscript[Y, 2][T], Subscript[Y, 4][T]}], "EventLocationMethod" -> "LinearInterpolation", "Method" -> {"SymplecticPartitionedRungeKutta", "DifferenceOrder" -> 4, "PositionVariables" -> {Subscript[Y, 1][T], Subscript[Y, 2][T]}}}, StartingStepSize -> 0.25, MaxSteps -> ∞]; ];
以下绘制庞加莱截面的图形. 收集到的数据可以在 Reap 结果的最后一部分找到，而点列表可以在第一部分找到：
Wolfram Language code: psdata = sdata[[-1, 1]]; ListPlot[psdata, Axes -> False, Frame -> True, AspectRatio -> 1]

ABC 流

以下加载 Arnold–Beltrami–Childress (ABC) 流的一个实例问题，用于对三维欧拉方程的层流中的混沌（chaos）进行建模：
Wolfram Language code: system = GetNDSolveProblem["ArnoldBeltramiChildress"]; eqs = system["System"]; vars = system["DependentVariables"]; icvars = vars /. T -> 0;
以下通过把某些右侧分量设置为零，定义系统的一个 Y1、Y2 分裂：
Wolfram Language code: Y1 = eqs; Y1[[2, 2]] = 0;Y1
Wolfram Language code: Y2 = eqs;Y2[[{1, 3}, 2]] = 0; Y2
以下定义隐式中点法：
Wolfram Language code: ImplicitMidpoint = {"ImplicitRungeKutta", "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients", "DifferenceOrder" -> 2, "ImplicitSolver" -> {"FixedPoint", AccuracyGoal -> 10, PrecisionGoal -> 10, "IterationSafetyFactor" -> 1}};
以下构造了二阶分裂法，它保持体积（volume）并且反转对称性：
Wolfram Language code: ABCSplitting = {"Splitting", "DifferenceOrder" -> 2, "Equations" -> {Y2, Y1, Y2}, "Method" -> {"LocallyExact", ImplicitMidpoint, "LocallyExact"}};
以下定义给出特定初始条件的庞加莱截面的一个函数：
Wolfram Language code: psect[ics_] := Module[{reapdata}, reapdata = Reap[ NDSolve[{eqs, Thread[icvars == ics]}, {}, {T, 1000}, Method -> {"EventLocator", "Event" -> Subscript[Y, 2][T], "EventAction" :> Sow[{Subscript[Y, 1][T], Subscript[Y, 3][T]}], "EventLocationMethod" -> "LinearInterpolation", Method -> ABCSplitting}, StartingStepSize -> 1 / 4, MaxSteps -> ∞] ]; reapdata[[-1, 1]] ];
以下对于一些不同的初始条件，寻找庞加莱截面，并且把它们简化（flatten）为一个一些点组成的单一列表：
Wolfram Language code: data = Mod[Map[psect, {{4.267682454609692, 0, 0.9952906114885919}, {1.6790790859443243, 0, 2.1257099470901704}, {2.9189523719753327, 0, 4.939152797323216}, {3.1528896559036776, 0, 4.926744120488727}, {0.9829282640373566, 0, 1.7074633238173198}, {0.4090394012299985, 0, 4.170087631574883}, {6.090600411133905, 0, 2.3736566160602277}, {6.261716134007686, 0, 1.4987884558838156}, {1.005126683795467, 0, 1.3745418575363608}, {1.5880780704325377, 0, 1.3039536044289253}, {3.622408133554125, 0, 2.289597511313432}, {0.030948690635763183, 0, 4.306922133429981}, {5.906038850342371, 0, 5.000045498132029}}], 2π]; ListPlot[data, ImageSize -> Medium]

弹跳球

这个例子是 [SGT03] 中一个例子的推广. 它对一个沿着具有给定形状的坡道弹跳的球进行建模. 这个例子很好地演示了用户可以如何使用具有事件定位的 NDSolve 的多次调用来对某些行为进行建模.

以下定义一个函数，计算从一次弹跳到下一次弹跳的解. 这里计算该解，直到下一次路径与坡道相交：
Wolfram Language code: OneBounce[k_, ramp_][{t0_, x0_, xp0_, y0_, yp0_}] := Module[{sol, t1, x1, xp1, y1, yp1}, sol = First[NDSolve[ {x''[t] == 0, x'[t0] == xp0, x[t0] == x0, y''[t] == -9.8, y'[t0] == yp0, y[t0] == y0}, {x, y}, {t, t0, ∞}, Method -> {"EventLocator", "Event" :> y[t] - ramp[x[t]]}, MaxStepSize -> 0.01]]; t1 = InterpolatingFunctionDomain[x /. sol][[1, -1]]; {x1, xp1, y1, yp1} = Reflection[k, ramp][{x[t1], x'[t1], y[t1], y'[t1]} /. sol]; Sow[{x[t] /. sol, t0 ≤ t ≤ t1}, "X"]; Sow[{ y[t] /. sol, t0 ≤ t ≤ t1}, "Y"]; Sow[{x1, y1}, "Bounces"]; {t1, x1, xp1, y1, yp1}]
这里定义当球击中坡道时的弹跳函数. 公式的建立是基于关于坡道法线的反射，假设弹跳后，只剩比例为 的能量：
Wolfram Language code: Reflection[k_, ramp_][{x_, xp_, y_, yp_}] := Module[{gramp, gp, xpnew, ypnew}, gramp = -ramp'[x]; If[Not[NumberQ[gramp]], Print["Could not compute derivative "]; Throw[$Failed]]; gramp = {-ramp'[x], 1}; If[ gramp.{xp, yp} == 0, Print["No reflection"]; Throw[$Failed]]; gp = {1, -1} Reverse[gramp]; {xpnew, ypnew} = (k / (gramp.gramp)) (gp gp.{xp, yp} - gramp gramp.{xp, yp}); {x, xpnew, y, ypnew}]
这里定义函数，它对于一个给定的反射率、坡道和起始位置，运行弹跳球的模拟：
Wolfram Language code: BouncingBall[k_, ramp_, {x0_, y0_}] := Module[{data, end, bounces, xmin, xmax, ymin, ymax}, If[y0 < ramp[x0], Print["Start above the ramp"]; Return[$Failed]]; data = Reap[ Catch[Sow[{x0, y0}, "Bounces"]; NestWhile[OneBounce[k, ramp], {0, x0, 0, y0, 0}, Function[1 - #1[[1]] / #2[[1]] > 0.01], 2, 25]], _, Rule]; end = data[[1, 1]]; data = Last[data]; bounces = ("Bounces" /. data); xmax = Max[bounces[[All, 1]]]; xmin = Min[bounces[[All, 1]]]; ymax = Max[bounces[[All, 2]]]; ymin = Min[bounces[[All, 2]]]; Show[{Plot[ramp[x], {x, xmin, xmax}, PlotRange -> {{xmin, xmax}, {ymin, ymax}}, Epilog -> {PointSize[.025], Map[Point, bounces]}, AspectRatio -> (ymax - ymin) / (xmax - xmin)], ParametricPlot[Evaluate[{Piecewise["X" /. data], Piecewise["Y" /. data]}], {t, 0, end}, PlotStyle -> StandardRed]}]]
这是 [SGT03] 中实现的例子：
Wolfram Language code: ramp[x_] := If[x < 1, 1 - x, 0]; BouncingBall[.7, ramp, {0, 1.25}]
现在，定义坡道为一个四分之一圆：
Wolfram Language code: circle[x_] := If[x < 1, Sqrt[1 - x ^ 2], 0]; BouncingBall[.7, circle, {.1, 1.25}]
这里对坡道增加一个轻微的波动：
Wolfram Language code: wavyramp[x_] := If[x < 1, 1 - x + .05 Cos[11 Pi x], 0]; BouncingBall[.75, wavyramp, {0, 1.25}]

事件方向

常微分方程

这个例子说明受限的三体问题（由四个方程组成的标准非刚性测试系统）的解. 这个例子追踪了绕着月球旅行、并且返回地球的飞船的路径（见 [SG75]的第246页）. 指明多个事件和零点交叉方向的能力是重要的.

初始条件的选择是使得轨道具有周期性. 的值对应于绕月球和地球旅行的飞船：
Wolfram Language code: μ = (1/82.45); μ^ * = 1 - μ; Subscript[r, 1] = Sqrt[(Subscript[y, 1][t] + μ)^2 + Subscript[y, 2][t]^2]; Subscript[r, 2] = Sqrt[(Subscript[y, 1][t] - μ^ * )^2 + Subscript[y, 2][t]^2]; eqns = {{Derivative[1][Subscript[y, 1]][t] == Subscript[y, 3][t], Subscript[y, 1][0] == 1.2}, {Derivative[1][Subscript[y, 2]][t] == Subscript[y, 4][t], Subscript[y, 2][0] == 0}, {Derivative[1][Subscript[y, 3]][t] == 2 Subscript[y, 4][t] + Subscript[y, 1][t] - (μ^ * (Subscript[y, 1][t] + μ)/Subsuperscript[r, 1, 3]) - (μ (Subscript[y, 1][t] - μ^ * )/Subsuperscript[r, 2, 3]), Subscript[y, 3][0] == 0}, {Derivative[1][Subscript[y, 4]][t] == -2 Subscript[y, 3][t] + Subscript[y, 2][t] - (μ^ * Subscript[y, 2][t]/Subsuperscript[r, 1, 3]) - (μ Subscript[y, 2][t]/Subsuperscript[r, 2, 3]), Subscript[y, 4][0] == -1.0493575098303199073}};
事件函数是离开初始条件的距离的导数. 当该值穿过零点时，出现一个局部最大值或者最小值：
Wolfram Language code: ddist = 2(Subscript[y, 3][t](Subscript[y, 1][t] - 1.2) + Subscript[y, 4][t]Subscript[y, 2][t]);
有两个事件，在这个例子中是相同的. 第一个事件（具有 Direction 1）对应于从初始点开始的距离是局部最小值的点，所以飞船返回初始位置. 该事件的动作是存储事件的时间到变量 tfinal 中，并且停止积分. 第二个事件对应于一个局部最大值. 事件动作是存储飞船离初始位置最远处的时间到变量 tfar 中：
Wolfram Language code: sol = First[NDSolve[eqns, {Subscript[y, 1], Subscript[y, 2], Subscript[y, 3], Subscript[y, 4]}, {t, ∞}, Method -> {"EventLocator", "Event" -> {ddist, ddist}, "Direction" -> {1, -1}, "EventAction" :> {Throw[tfinal = t, "StopIntegration"], tfar = t}}]]

前两个解的组成部分是无穷小质量的物体的坐标，所以画出一个对另一个的图形就能给出物体的轨道.

以下显示当飞船位于距离初始位置最远的点时的半个轨道：
Wolfram Language code: ParametricPlot[{Subscript[y, 1][t], Subscript[y, 2][t]} /. sol, {t, 0, tfar}]
以下显示当飞船返回初始位置时的一个完整的轨道：
Wolfram Language code: ParametricPlot[{Subscript[y, 1][t], Subscript[y, 2][t]} /. sol, {t, 0, tfinal}]

延迟微分方程

下面系统对传染病进行建模（参阅 [HNW93]、[ST00] 和 [ST01]）：
Wolfram Language code: system = {y1'[t] == -y1[t] y2[t - 1] + y2[t - 10], y1[t /; t ≤ 0] == 5, y2'[t] == y1[t] y2[t - 1] - y2[t], y2[t /; t ≤ 0] == 1 / 10, y3'[t] == y2[t] - y2[t - 10], y3[t /; t ≤ 0] == 1}; vars = {y1[t], y2[t], y3[t]};
当积分进行时，对每个组成部分的一个局部最大值的数据进行收集. 用 Sow 和 Reap 的一个不同的标签来区分组成部分：
Wolfram Language code: data = Reap[ sol = First[NDSolve[system, vars, {t, 0, 40}, Method -> {"EventLocator", "Event" :> {y1'[t], y2'[t], y3'[t]}, "EventAction" :> {Sow[{t, y1[t]}, 1], Sow[{t, y2[t]}, 2], Sow[{t, y3[t]}, 3]}, "Direction" -> {-1, -1, -1}}]], {1, 2, 3} ];
同时显示局部最大值，以及解的组成部分：
Wolfram Language code: colors = {{StandardRed}, {StandardBlue}, {StandardGreen}}; plots = Plot[Evaluate[vars /. sol], {t, 0, 40}, PlotStyle -> colors]; max = ListPlot[Part[data, -1, All, 1], PlotStyle -> colors]; Show[plots, max]

不连续方程和切换函数

在许多应用中，微分系统中的函数可能并不是处处解析的或者连续的.

实际应用中常见的不连续性问题包含一个切换函数 ：

为了说明穿过一个不连续点时的困难，考虑以下例子 [GØ84]（也请参阅 [HNW93]）：

下面是整个系统的输入. 切换函数赋给了符号事件 event，而定义该系统的函数取决于切换函数的符号：
Wolfram Language code: t0 = 0; ics0 = (3/10); event = (t + (1/20))^2 + (y[t] + (3/20))^2 - 1; system = {y'[t] == If[event <= 0, t^2 + 2y[t]^2, 2t^2 + 3y[t]^2 - 2], y[t0] == ics0};
符号 odemethod 用于表明应该用于积分的数值方法. 作为比较，用户可能想要定义一个不同的方法，比如 "ExplicitRungeKutta"，并且重新运行本节中的计算，看看其他方法如何工作：
Wolfram Language code: odemethod = Automatic;
以下求解区间 [0, 1] 上的系统，并且使用 Reap 和 Sow 收集积分的网格点的数据：
Wolfram Language code: data = Reap[ sol = y[t] /. First[NDSolve[system, y, {t, t0, 1}, Method -> odemethod, MaxStepFraction -> 1, StepMonitor :> Sow[t]]] ][[2, 1]]; sol
以下是解的图形：
Wolfram Language code: dirsol = Plot[sol, {t, t0, 1}]

尽管事实是已经得到一个解，但是我们还不清楚这个解是否高效地得到.

下面的例子表明了对于数值求解器，穿过不连续点的困难之处.

以下定义了一个函数，它显示了积分的网格点以及采用的积分步的数目：
Wolfram Language code: StepPlot[data_, opts___ ? OptionQ] := Module[{sdata}, sdata = Transpose[{data, Range[Length[data]]}]; ListPlot[sdata, opts, Axes -> False, Frame -> True, PlotRange -> All] ];
当积分穿过不连续点（值接近0.6），积分方法将遇到困难，并且采用大量的小步长的步进——有时候也可以观察到骤一些的被拒步骤：
Wolfram Language code: StepPlot[data]

穿过不连续点的最有效的方法之一是通过在不连续点中断积分，重新开始.

以下例子表明如何使用 "EventLocator" 方法实现这个目的.

以下对该系统的第一部分进行数值积分，直到不连续点. 切换函数作为事件给出. 该事件的方向限制为从负值到正值的改变. 当找到该事件时，解和事件时间由事件动作存储：
Wolfram Language code: system1 = {y'[t] == t^2 + 2y[t]^2, y[t0] == ics0}; data1 = Reap[sol1 = y[t] /. First[NDSolve[system1, y, {t, t0, 1}, Method -> {"EventLocator", "Event" -> event, Direction -> 1, EventAction :> Throw[t1 = t; ics1 = y[t]; , "StopIntegration"], Method -> odemethod}, MaxStepFraction -> 1, StepMonitor :> Sow[t]]] ][[2, 1]]; sol1

使用 "EventLocator" 方法找到的不连续点作为一个新的初始条件，现在积分可以继续.

以下定义一个系统和初始条件，数值求解该系统，并且收集用于网格点的数据：
Wolfram Language code: system2 = {y'[t] == 2t^2 + 3y[t]^2 - 2, y[t1] == ics1}; data2 = Reap[ sol2 = y[t] /. First[NDSolve[system2, y, {t, t1, 1}, Method -> odemethod, MaxStepFraction -> 1, StepMonitor :> Sow[t]]] ][[2, 1]]; sol2
两个解的图形与通过求解整个系统所得的非常相似：
Wolfram Language code: evsol = Plot[If[t ≤ t1, sol1, sol2], {t, 0, 1}]
检查网格点，显然该方法采用了少得多的步进，并且已经消除了接近不连续点处遇到的有问题的行为：
Wolfram Language code: StepPlot[Join[data1, data2]]

不连续点的值在 [HNW93] 中给出为 0.6234，这与 "EventLocator" 方法找到的值一致.

在本例中，有可能解析求解该系统，并且使用数值方法来检查值.

该系统至不连续点的解可以使用贝塞尔和伽马函数表示：
Wolfram Language code: dsol = FullSimplify[First[DSolve[system1, y[t], t]]]
把解代入切换函数，一个局部极小化确认了不连续点的值：
Wolfram Language code: FindRoot[event /. dsol, {t, 3 / 5}]

避免偏微分方程中的 Wraparound

许多演化方程对这样一个空间域内的行为建模，该空间域是无限的或者足够大以至于无法不使用专门的离散方法对整个域进行离散化. 在实践中，通常的情况是，只要我们感兴趣的解依然是局域化的，我们就有可能一直使用较小的计算域.

当计算域的边界由实际考虑制定，而不是根据研究的实际模型制定时，我们就可能适当地选择边界条件. 因为周期性拟谱逼近的良好的分辨率，使用具有周期性边界条件的拟谱方法可以增加计算域的范围. 周期性边界条件的缺点是传播过边界的信号在域的另一边继续，从而通过wraparound影响了解. 我们有可能在接近边界处使用一个吸收层以最小化这些效果，但是并不总是可能完全消除它们.

sine-Gordon 方程出现在微分几何和相对场论中. 这个例子从一个局部化的初始条件开始，对方程进行积分.周期性的拟谱方法用于积分. 由于没有在边界附近设置吸收层，最适当的办法是，一旦wraparound变得显著，即停止积分. 这种条件很容易通过使用 "EventLocator" 方法，用事件定位检测到.

当周期性wraparound点处的解的大小超过临界值0.01时，积分停止，超过临界值时，波的形式将受到周期性的影响：
Wolfram Language code: Timing[sgsol = First[NDSolve[{Subscript[∂, t, t]u[t, x] == Subscript[∂, x, x]u[t, x] - Sin[u[t, x]], u[0, x] == E^-(x - 5)^2 + E^-(x + 5)^2 / 2, u^(1, 0)[0, x] == 0, u[t, -50] == u[t, 50]}, u, {t, 0, 1000}, {x, -50, 50}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"}, Method -> {"EventLocator", "Event" :> Abs[u[t, -50]] - 0.01, "EventLocationMethod" -> "StepBegin"}}]]]
以下从 InterpolatingFunction 对象提取结束时间，并且对计算所得的解生成图形. 用户可以看到第一组波开始到达边界时，积分已经停止：
Wolfram Language code: end = InterpolatingFunctionDomain[u /. sgsol][[1, -1]]; DensityPlot[u[t, x] /. sgsol, {x, -50, 50}, {t, 0, end}, Mesh -> False, PlotPoints -> 100]
"DiscretizedMonitorVariables" 选项影响对偏微分方程事件解释的方式；使用设置 True，u[t,x] 由离散值组成的一个向量代替. 这样更加有效，因为它避免了明确构建InterpolatingFunction 来对事件进行计算：
Wolfram Language code: Timing[sgsol = First[NDSolve[{Subscript[∂, t, t]u[t, x] == Subscript[∂, x, x]u[t, x] - Sin[u[t, x]], u[0, x] == E^-(x - 5)^2 + E^-(x + 5)^2 / 2, u^(1, 0)[0, x] == 0, u[t, -50] == u[t, 50]}, u, {t, 0, 1000}, {x, -50, 50}, Method -> {"MethodOfLines", "DiscretizedMonitorVariables" -> True, "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"}, Method -> {"EventLocator", "Event" :> Abs[First[u[t, x]]] - 0.01, "EventLocationMethod" -> "StepBegin"}}]]]

性能比较

下面的例子构建了一个表格，对两种不同的积分方法进行比较.

这里定义一个函数，它返回计算一个轻度阻尼摆方程的解所需的时间，直至摆动经历了1000次的瞬时静止后的时间点：
Wolfram Language code: EventLocatorTiming[locmethod_, odemethod_] := Block[{Second = 1, y, t, p = 0}, First[Timing[NDSolve[{y''[t] + (1/1000) y'[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0}, y, {t, ∞}, Method -> {"EventLocator", "Event" -> y'[t], "EventAction" :> If[p++ ≥ 1000, Throw[end = t, "StopIntegration"]], "EventLocationMethod" -> locmethod, "Method" -> odemethod}, MaxSteps -> ∞]]] ];
以下使用函数来生成一个表格，对两种不同的常微分方程积分方法进行不同的定位方法的比较：
Wolfram Language code: elmethods = {"StepBegin", "StepEnd", "LinearInterpolation", {"Brent", "SolutionApproximation" -> "CubicHermiteInterpolation"}, Automatic}; odemethods = {Automatic, "ExplicitRungeKutta"}; TableForm[Outer[EventLocatorTiming, elmethods, odemethods, 1], TableHeadings -> {elmethods, odemethods}, TableSpacing -> {1, .25}]

虽然简单的步进开始/结束以及线性插值定位的代价基本一样低，更好的定位方法则更加昂贵. 默认的定位方法对于显式朗格-库塔方法尤其昂贵，因为它还不支持一个连续的输出公式——因此，它需要在局域极小化过程中反复调用具有不同步长的方法.

值得指出的是，往往计算事件的额外时间的重要组成部分来自在每个时间步进对事件函数进行计算的需求，以便检查可能的符号变化.

Wolfram Language code: TableForm[ {Map[ Block[{Second = 1, y, t, p = 0}, First[Timing[NDSolve[{y''[t] + (1/1000) y'[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0}, y, {t, end}, Method -> #, MaxSteps -> ∞]]]]&, odemethods ]}, TableHeadings -> {None, odemethods}]

对于仅包含自变量的事件函数，有一个最优化操作. 这样的事件在初始化时被自动检测. 例如，这样做的优势是应变量的解的插值在局部优化搜索的每个步进中并不执行——而是延迟直到找到自变量的值以后.

局限性

事件定位方法的一个局限性是：由于事件函数在一个步进区间内仅检查符号改变，如果在一个步进区间内该事件函数具有多个根，那么所有的或者一些事件可能被遗漏. 这种情况通常仅当常微分方程的解比事件函数变化得更慢时才发生. 当用户怀疑这种情况可能发生时，最简单的解决方案是通过使用 NDSolve 的 MaxStepSize 选项，降低方法所能采用的最大步长. 我们也可以采用成熟的方法，但是最好的方法取决于计算的内容. 接下来的一个例子表明了这个问题，并且给出两种方法来解决这个问题.

这里应该计算小于 的正整数数目（有 148个）. 但是，大部分都被漏掉了，因为解 x[t] 很简单，所以方法采用较大的步长：
Wolfram Language code: Block[{n = 0}, NDSolve[{y'[t] == y[t], y[-1] == E^-1}, y, {t, 5}, Method -> {"EventLocator", "Event" -> Sin[π y[t]], "EventAction" :> n++}];n]
这里限制最大步长，所以所有的事件都可以找到：
Wolfram Language code: Block[{n = 0}, NDSolve[{y'[t] == y[t], y[-1] == E^-1}, y, {t, 5}, Method -> {"EventLocator", "Event" -> Sin[π y[t]], "EventAction" :> n++}, MaxStepSize -> 0.001];n]

从实例问题的本质可以很明显地看出，如果端点增加，则很有可能需要一个较小的最大步长. 在各个位置采用很小的步长，效率将很低. 我们可能通过建立一个按与事件函数相同的时间变化尺度变化的变量，引进一个自适应时间步进的限制.

这里引进另一个函数：事件函数， 进行积分. 通过这种修改，并且允许方法采用尽可能多的步进，我们有可能在一个至多为 的合理的时间内，找到正确的值：
Wolfram Language code: Block[{n = 0}, NDSolve[{y'[t] == y[t], y[-1] == E^-1, z'[t] == D[Sin[π y[t]], t], z[-1] == Sin[π E^-1]}, {y, z}, {t, 10}, Method -> {"EventLocator", "Event" -> z[t], "EventAction" :> n++}, MaxSteps -> ∞];n]

选项总结

“EventLocator” 选项

选项名
默认值
"Direction"All对该事件所允许的零点交叉的方向；1 意味着从负值到正值；-1 意味着从正值到负值；而 All 包括两个方向
"Event"None定义事件的一个表达式；一个事件发生在这样的点：即当替换问题变量的数值值时，表达式等于零
"EventAction"Throw[Null, "StopIntegration"]当一个事件发生时要采取的动作：在该事件中，问题变量使用它们的数值值代替；一般说来，用户需要使用RuleDelayed（）来防止选项被计算，除非选项具有数值值
"EventLocationMethod"Automatic细化给定事件的定位所采用的方法
"Method"Automatic对常微分方程组积分所采用的方法

"EventLocator" 方法选项.

“EventLocationMethod” 选项

"Brent"使用选择 Method->"Brent" 的 FindRoot 对事件进行定位；这是具有设置 Automatic 的默认值
"LinearInterpolation"使用线性插值法对事件事件时间进行定位；接着使用三次厄米（ Hermite） 插值法寻找该事件时间上的解
"StepBegin"事件由步骤开始处的解给出
"StepEnd"事件由步骤末尾处的解给出

"EventLocationMethod" 选项的设置.

“Brent” 选项

选项名
默认值
"MaxIterations"100在方法的一个步进中，用来定位一个事件所用的最大迭代次数
"AccuracyGoal"Automatic传递给 FindRoot 的准确度目标设置；如果是 Automatic，那么传递给 FindRoot 的值是基于 NDSolve 的局部误差设置的
"PrecisionGoal"Automatic传递给 FindRoot 的精度目标设置；如果是 Automatic，那么传递给 FindRoot 的值是基于 NDSolve 的局部误差设置的
"SolutionApproximation"Automatic如何在细化过程中对解求近似值以实现事件函数的计算；可以是 Automatic 或者 "CubicHermiteInterpolation"

事件定位方法 "Brent" 的选项.