NDSolveの"EventLocator"メソッド
はじめに
微分方程式系の変化を検出し，その位置を正確に知ることができると便利なことがよくある．例えば，特異性や状態変化の検出では，積分のやり直し等の適切な操作が行われる．
のようにゼロになる解の付近の点のことである．あるいは，ブール値のイベント関数を考えることもできる．この場合，イベントは関数がTrueからFalseへ，あるいはその逆に変化するときに起る．
NDSolveに組み込まれている"EventLocator"メソッドは，コントローラ法として効率的に動作する．これはイベントのチェックおよび適切なアクションの実行を司るが，そうでない場合は微分方程式系の積分は，内在するメソッドに完全に委ねられる．
このセクションには，"EventLocator"メソッドとそのオプションの基本的な使用方法を示す例題が含まれている．後続のセクションでは，期間の検出，ポアンカレ(Poincaré)断面，不連続の取扱い等の高度な応用例を示す．
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["GUIKit`"];簡単な例として，非平衡位置で開始した振子が最低点を通過する時等のイベントを見付けて，その時点で積分を中止するというものがある．
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オブジェクトからその値を抽出することができる．
end = InterpolatingFunctionDomain[First[y /. sol]][[1, -1]];
Plot[Evaluate[{y[t], y'[t]} /. First[sol]], {t, 0, end}, PlotStyle -> {{StandardBrown}, {StandardBlue}}]イベント位置特定メソッドを使うときは，"EventLocator"メソッドのメソッドオプションを使って，位置特定の対象となるイベントと，イベントを見付けたときのアクションを指定する．
イベントが検出されたときのデフォルトのアクションは，上のように積分を中止することである．イベントアクションはどのような式でもよい．これは，イベントが検出されるたびに問題の変数に数値を代入して評価される．
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[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)に制限することができる．
NDSolve`stop = False;
GUIRun[Widget["Panel", {Widget["Button", {
"label" -> "Stop",
BindEvent["action",
Script[NDSolve`stop = True]]}]}]];NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 1, y'[0] == 0}, y, {t, 0, ∞}, Method -> {"EventLocator", "Event" :> NDSolve`stop}, MaxSteps -> ∞]
上記の例ではstopの即時型の値が評価され，関数として設定されないようにするために，RuleDelayed (:>)を使って"Event"オプションを指定している．
複数のイベントを指定することもできる．イベント関数を数値的評価してリストになった場合，リストの各要素は別のイベントとみなされる．各イベントについて異なるアクション，方向等を指定したい場合は，そのオプションの値を，適切な長さのリストとして指定するとよい．
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"オプションを使うと，イベントのテストで満足しなければならない追加のブール条件を指定することができる．ブールイベントの代りにこれを使った方が根の探索がより効率的に行われるため，有利である．
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"]}]Plot[Evaluate[{y[t], y'[t], y'[t] ^ 2 / 2 - Cos[y[t]], -.9} /. First[sol]], {t, 0, end}, PlotStyle -> {{StandardBrown}, {StandardBlue}, {StandardGreen}, {StandardRed}}]"EventLocator"のMethodオプションを使うと，積分で使用する数値メソッドの指定が可能である．
イベントの位置特定メソッド
"EventLocator"メソッドは，内在するメソッドのステップを取り，イベント関数のいずれかの符号（パリティ）がステップの終点で異なるかどうかをチェックすることにより動作する．イベント関数には，実数値かブール値が想定されているため，符号の変化があれば，そのステップ区間でイベントがあったことになる．イベント関数で，1ステップ内にイベントを含むものについては，その区間内におけるイベントの位置を特定するために，絞込み操作が行われる．
位置の絞込みには，いくつかの異なるメソッドがある．このようなメソッドには，積分区間の最初と最後で単純に解を取るというもの，イベント値の線形補間，小区間内の根の探索メソッドの使用等がある．適切なメソッドは，実行速度か位置の正確さのどちらを犠牲にするかによる．
イベントアクションが積分の中止となっている場合は，その積分が中止される特定の値は，"EventLocator"のEventLocationMethodオプションで得られる値に依存する．
1つのイベントの位置は通常迅速に見付かるので，使用するメソッドが全体の計算時間に大きく影響を与えることはない．しかし，イベントが複数回検出されると，位置の絞込みメソッドにより影響が出る可能性がある．
"StepBegin"メソッドと"StepEnd"メソッド
イベントの正確な位置は大して問題ではない場合，あるいはイベントの正確な位置が分かっても，内在する計算に確度を持ったものが何も反映されない場合は，最も大雑把なメソッドが適している．前セクションのStopボタンの例がそれである．時間ステップは非常に迅速に計算されるので，一定の時間ステップ内でボタンをクリックするのにかかる時間を計測する方法はない．ましてや1つの時間ステップ内の特定の点における時間等測ることはできない．従って，継承されたイベントの精度に基づいている場合は，絞込みを行う意味がない．これは"StepBegin"あるいは"StepEnd"位置特定メソッドを使って指定することができる．イベントの定義がヒューリスティックであったり，幾分不明瞭であったりした場合は，このメソッドが適切である．
"LinearInterpolation"メソッド
グラフにプロットする点のためにイベントの結果が必要な場合は，グラフの解像度に合わせてイベントの位置を見付けるだけでよい．これにはステップの終点を使うだけでは通常大雑把すぎるが，イベント関数の値に基づく単独の線形補間で十分である．
数値積分の連続したメッシュ点におけるイベント関数を以下のように表す．
線形補間は，イベント時間における解の近似にも使うことができる．しかし，導関数値 
および 
がメッシュ点で利用可能なので，イベントにおける解のよりよい近似は，次のように三次エルミート補間を使って安価に計算することができる．
設定を"LinearInterpolation"とすると，単独の線形補間に基づく絞込みが指定できる．
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}}]周期全体に渡るプロットの解像度において，導関数が軸と交わる正確な点が終点にはならないことは分からない．しかし，十分ズームすると，誤差が見える．
Plot[Evaluate[y'[t] /. sol], {t, end * (1 - .001), end}, PlotStyle -> StandardBlue]線形補間法は，以下のセクションで示すポアンカレ断面の例のように，見ることを目的とするほとんどの場合には十分な方法である．しかし，ブール値を取るイベント関数では，線形補間は実質的には2分法の1ステップでしかないので，線形補間法はグラフィックスには不適切である場合がある．
ブレント法(Brent's Method)
デフォルトの位置特定メソッドは，イベント位置特定メソッドの"Brent"である．これはFindRootでブレント法を使ってイベントの位置を見付けるものである．ブレント法は小区間の根から開始して，補間法と2分法に基づくステップを組み合せることによって，少なくとも2分法と同じくらいの収束率を保証する．FindRootがイベント関数の根をどの程度の確度・精度で見付けるようにするかは，イベント位置特定メソッドのオプション"Brent"を使って制御することができる．デフォルトは，NDSolveが局所誤差制御で使っているのと同じ確度・精度で根を見付けることである．
密な出力や難解な出力をサポートするメソッドでは，イベント関数の引数は，連続出力公式を使うことにより非常に効率よく求めることができる．しかし，連続出力をサポートしないメソッドでは，内在するメソッドのステップを取ることで解を計算する必要があるため，比較的高価になる可能性がある．解の近似を求める別の方法で，メソッドの順序通りではないが，NDSolveから返されるInterpolatingFunctionオブジェクトにFindRootを使うことに見合ったものとして，三次エルミート(Hermite)補間がある．これは，ステップの終点における解の値と解の導関数値に基づく補間によって，ステップの中盤で解の近似値を得るものである．
比較
以下の例では，多くの異なるイベント位置特定メソッドで振り子方程式を積分し，イベントが見付かった時間を比較する．
eventmethods = {"StepBegin", "StepEnd", "LinearInterpolation", Automatic};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]を参照）．
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"]}]]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]ファン・デル・ポル（van der Pol）発振器の周期
ファン・デル・ポル発振器は，極度に硬い常微分方程式系の例である．イベントの位置特定メソッドは，常微分方程式系の積分を実際に行うのにどのようなメソッドを呼び出すこともできる．デフォルトメソッドのAutomaticは必要に応じて硬い系に適したメソッドに切り替えるので，硬さによる問題はない．
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の解の終点を選ぶことにより，周期を返す関数を 
の関数として書くことができる．
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]]];
vper[1000]もちろん，周期解を使ってこのメソッドをすべての系に一般化することは簡単である．
ポアンカレ断面
ポアンカレ断面は，高次元の微分方程式系の解を視覚化するのに便利である．
インタラクティブなグラフィカルインターフェースについては，EquationTrekkerパッケージをご参照いただきたい．
ヘノン–ハイレス(Hénon–Heiles)系
ここでは銀河系の恒星運動をモデル化するヘノン–ハイレス系を定義する．
system = GetNDSolveProblem["HenonHeiles"];
vars = system["DependentVariables"];eqns = {system["System"], system["InitialConditions"]}この場合，関心のあるポアンカレ断面は，軌道が 
を通過するときの 
平面の点の集合である．
数値積分の実際の結果は必要ではないので，出力変数リスト（NDSolveの第2引数）を空，つまり{}と指定することにより，InterpolatingFunctionにすべてのデータを保存することを避けることができる．こうすると，NDSolveが出力にInterpolatingFunctionを生成しないため，多量の不必要なデータを保存する必要がなくなるのである．NDSolveはNDSolve::nooutというメッセージで出力関数がないという警告を出すが，この場合は必要なデータがイベントアクションから集められているので，メッセージをオフにしても差し支えない．
ここでの計算目的は，比較的低解像度のグラフで結果を見ることにあるため，線形補間によるイベント位置特定メソッドが使われる．ポアンカレ地図の固定点等，詳細を見たり，特徴を調べたりするためにグラフをズームする必要がある例題を行っている場合は，デフォルトの位置特定メソッドを使った方が適切であろう．
Off[NDSolve::noout];
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 -> ∞];
];psdata = data[[-1, 1]];
ListPlot[psdata, Axes -> False, Frame -> True, AspectRatio -> 1]ヘノン-ハイレス系はハミルトニアン(Hamiltonian)なので，この例ではシンプレクティック法を使った方が良質な結果が得られる．
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 -> ∞];
];psdata = sdata[[-1, 1]];
ListPlot[psdata, Axes -> False, Frame -> True, AspectRatio -> 1]ABCフロー
system = GetNDSolveProblem["ArnoldBeltramiChildress"];
eqs = system["System"];
vars = system["DependentVariables"];
icvars = vars /. T -> 0;Y1 = eqs; Y1[[2, 2]] = 0;Y1Y2 = eqs;Y2[[{1, 3}, 2]] = 0; Y2ImplicitMidpoint = {"ImplicitRungeKutta", "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients", "DifferenceOrder" -> 2, "ImplicitSolver" -> {"FixedPoint", AccuracyGoal -> 10, PrecisionGoal -> 10, "IterationSafetyFactor" -> 1}};ABCSplitting = {"Splitting",
"DifferenceOrder" -> 2,
"Equations" -> {Y2, Y1, Y2},
"Method" -> {"LocallyExact", ImplicitMidpoint, "LocallyExact"}};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]]
];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を位置特定メソッドで使って，ある動作をモデル化するよい例である．
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}]
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}]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]}]]ramp[x_] := If[x < 1, 1 - x, 0];
BouncingBall[.7, ramp, {0, 1.25}]circle[x_] := If[x < 1, Sqrt[1 - x ^ 2], 0];
BouncingBall[.7, circle, {.1, 1.25}]wavyramp[x_] := If[x < 1, 1 - x + .05 Cos[11 Pi x], 0];
BouncingBall[.75, wavyramp, {0, 1.25}]イベントの向き
常微分方程式
この例では，4つの方程式の標準的な硬くないテストシステムである制限三体問題の解を示す．この例は，月の周囲を回って地球に戻る宇宙船の軌道を追跡する（[SG75]の246ページを参照）．複数イベントおよび零交差の方向が指定できるということが重要である．
μ = (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}};ddist = 2(Subscript[y, 3][t](Subscript[y, 1][t] - 1.2) + Subscript[y, 4][t]Subscript[y, 2][t]);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}}]]最初の2つの解の要素は，無限小量の物体の座標である．従って，ひとつの要素に対して別の要素をプロットすると，物体の軌道が得られる．
ParametricPlot[{Subscript[y, 1][t], Subscript[y, 2][t]} /. sol, {t, 0, tfar}]ParametricPlot[{Subscript[y, 1][t], Subscript[y, 2][t]} /. sol, {t, 0, tfinal}]遅延微分方程式
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]};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}
];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]も参照のこと）．
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 = Automatic;data = Reap[
sol = y[t] /. First[NDSolve[system, y, {t, t0, 1}, Method -> odemethod, MaxStepFraction -> 1, StepMonitor :> Sow[t]]]
][[2, 1]];
soldirsol = Plot[sol, {t, t0, 1}]解が実際に得られたが，効率的に得られたのかどうかは明らかではない．
以下の例は，不連続性の通過により数値ソルバの使用が困難になることを示している．
StepPlot[data_, opts___ ? OptionQ] :=
Module[{sdata},
sdata = Transpose[{data, Range[Length[data]]}];
ListPlot[sdata, opts, Axes -> False, Frame -> True, PlotRange -> All]
];StepPlot[data]不連続性を通過する最も有効なメソッドのひとつに，積分を中止して不連続点で再開するというものがある．
以下の例では，これを行うためにどのように"EventLocator"メソッドを使えばよいかを例示する．
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"メソッドで見付かった不連続性を新しい初期条件として使うと，連続積分できるようになる．
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]];
sol2evsol = Plot[If[t ≤ t1, sol1, sol2], {t, 0, 1}]StepPlot[Join[data1, data2]]不連続の値は[HNW93]では0.6234で与えてある．これは"EventLocator"メソッドで見付かった値と偶然にも一致している．
以下の例では，系を解析的に解き，数値法を使ってその値をチェックすることが可能である．
dsol = FullSimplify[First[DSolve[system1, y[t], t]]]FindRoot[event /. dsol, {t, 3 / 5}]偏微分方程式での折返しの回避
多くの発展方程式は，特化された離散化メソッドを使わずに全領域を離散化することができなくなるほどの，無限もしくは十分大きい空間領域上で動作する．実際には，興味のある解が局所化されている限り，小さい計算領域を使うことができる場合もある．
計算領域の境界が，研究中の実際のモデルではなく実践的であるかどうかによって決められている場合，境界条件を適切に選ぶことができる．擬スペクトル法を周期境界条件で使用すると，その優れた近似結果のため，計算領域の範囲を拡張することが可能である．周期境界条件の欠点は，境界を超えて伝播する信号が領域の反対側に残り，折返しによる解への影響が出る．この影響を最小限にするために，境界付近に吸収層を使うこともできるが，常に完全に排除できるというわけではない．
サインゴルドン(Sin-Gordon)方程式は，微分幾何学および相対論的場の理論で現れる．以下の例では，散らばっている局所化された初期条件から始めて，サインゴルドン方程式を積分する．積分には，周期的な擬スペクトル法が使われる．境界付近には吸収層を置いていないので，折返しが顕著になった時点で積分を中止するのが最も妥当である．この条件は"EventLocator"メソッドを使ったイベント位置特定で，簡単に検出できる．
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"}}]]]end = InterpolatingFunctionDomain[u /. sgsol][[1, -1]];
DensityPlot[u[t, x] /. sgsol, {x, -50, 50}, {t, 0, end}, Mesh -> False, PlotPoints -> 100]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"}}]]]パフォーマンスの比較
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 -> ∞]]]
];elmethods = {"StepBegin", "StepEnd", "LinearInterpolation", {"Brent", "SolutionApproximation" -> "CubicHermiteInterpolation"}, Automatic};
odemethods = {Automatic, "ExplicitRungeKutta"};
TableForm[Outer[EventLocatorTiming, elmethods, odemethods, 1], TableHeadings -> {elmethods, odemethods}, TableSpacing -> {1, .25}]単純なステップ開始／終了法と線形補間法は同様にコストが低いが，これよりもよい位置特定メソッドになると高価になる．デフォルトの位置特定メソッドは，連続出力公式をサポートしていないので，陽的ルンゲ・クッタ法では特に高価になる．従って，デフォルトメソッドは，局所的な最小化の過程で，異なる刻み幅で繰返しメソッドを起動される必要がある．
ここで特筆すべきことは，イベントの計算にかかる付加的な時間の大分部は，符号の変化の可能性をチェックするために，各時間ステップでイベント関数を評価する必要があるということにより生じる．
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オプションを使って，メソッドが取れる最大刻み幅を削減するというものが最も簡単な解決策である．高度な方法もあるが，最善策は何を計算しているかによる．以下の例では，この問題とそれに対する2種類の解決策を示す．
Block[{n = 0}, NDSolve[{y'[t] == y[t], y[-1] == E^-1}, y, {t, 5}, Method -> {"EventLocator", "Event" -> Sin[π y[t]], "EventAction" :> n++}];n]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]この例から，終点が大きくなるにつれ，最大刻み幅を小さくしていく必要があるという可能性が高いということは明らかである．いたるところで非常に小さい刻み幅を取るのは，極めて非効率的である．イベント関数と同じ時間スケールで変化する変数を設定することにより，適応型の時間刻み制限を導入することが可能である．
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 | 常微分方程式系を積分するのに使うメソッド |
"EventLocationMethod"オプション
| "Brent" | Method->"Brent"でFindRootを使い，イベントの位置を特定する．これは設定Automaticでのデフォルトである |
| "LinearInterpolation" | 線形補間を用いてイベント時間を特定する．この後，イベント時間での会を見付けるために，三次エルミート補間が使われる |
| "StepBegin" | イベントはステップの最初の解で与えられる |
| "StepEnd" | イベントはステップの最後の解で与えられる |
"Brent"オプション
オプション名
|
デフォルト値
| |
| "MaxIterations" | 100 | メソッドの1ステップ内でイベントを見付けるために使う最大反復回数 |
| "AccuracyGoal" | Automatic | FindRootに渡される目標確度設定．Automaticのとき，FindRootに渡される値は，NDSolveの局所誤差設定に基づく |
| "PrecisionGoal" | Automatic | FindRootに渡される目標精度設定．Automaticのとき，FindRootに渡される値は，NDSolveの局所誤差設定に基づく |
| "SolutionApproximation" | Automatic | 絞込みの過程でイベント関数を評価するための解を近似する方法．Automaticあるいは"CubicHermiteInterpolation"のいずれかを取る |