NDSolveの"EventLocator"メソッド

はじめに

微分方程式系の変化を検出し,その位置を正確に知ることができると便利なことがよくある.例えば,特異性や状態変化の検出では,積分のやり直し等の適切な操作が行われる.

微分方程式系のイベント

とは,実数値のイベント関数が

のようにゼロになる解の付近の点のことである.あるいは,ブール値のイベント関数を考えることもできる.この場合,イベントは関数がTrueからFalseへ,あるいはその逆に変化するときに起る.

NDSolveに組み込まれている"EventLocator"メソッドは,コントローラ法として効率的に動作する.これはイベントのチェックおよび適切なアクションの実行を司るが,そうでない場合は微分方程式系の積分は,内在するメソッドに完全に委ねられる.

このセクションには,"EventLocator"メソッドとそのオプションの基本的な使用方法を示す例題が含まれている.後続のセクションでは,期間の検出,ポアンカレ(Poincaré)断面,不連続の取扱い等の高度な応用例を示す.

以下の初期化コマンドにより,解決する微分方程式を含んでおり,効用関数を定義する便利なパッケージをロードする:

簡単な例として,非平衡位置で開始した振子が最低点を通過する時等のイベントを見付けて,その時点で積分を中止するというものがある.

y[t]が軸と交差する最初の点まで,振子方程式を積分する:

この結果から,振子はおよそ で最低点 に達することが分かる.InterpolatingFunctionAnatomyパッケージを使うと,InterpolatingFunctionオブジェクトからその値を抽出することができる.

以下では,イベントが起る点を抽出し,その点までの解(黒)とその導関数(青)をプロットしている:

イベント位置特定メソッドを使うときは,"EventLocator"メソッドのメソッドオプションを使って,位置特定の対象となるイベントと,イベントを見付けたときのアクションを指定する.

イベントが検出されたときのデフォルトのアクションは,上のように積分を中止することである.イベントアクションはどのような式でもよい.これは,イベントが検出されるたびに問題の変数に数値を代入して評価される.

減衰振子でイベント が検出されるたびに,時間と値が出力される:

この例では,イベントが見付かったとき以外は"EventAction"オプションが評価しないようにするために,オプションの指定にRuleDelayed ()が使ってある点に注目のこと.

出力結果から分かるように,アクションにより積分が中止されない場合は,複数のイベントが見付かる.イベント式の符号が変化するときにイベントが検出される."Direction"オプションを使うと,符号が特定の向きに変化するときだけに,イベントを制限することができる.

以下は,減衰振子の速度が負から正に変化する点を集める.ReapSowはプログラミングコントラクトであり,初めにどれだけのデータがあるかが分かっていないときにデータを収集するのに便利である.Reap[expr]は,評価中にSowが適用された式すべてと,expr の値を返す.ここではReapNDSolveを囲んでおり,Sowがイベントアクションの一部になっているため,各イベントのデータを収集することができる:

前の例の出力からお気付きかもしれないが,導関数がほぼ零になるときだけ,イベントが検出される.メソッドが,根底にある積分器のステップにおいて(イベント式の符号の変化によって)イベントの存在を検出したら,メソッドは数値法を使ってイベントの位置を近似的に見付ける.位置特定のプロセスは数値的なので,近似的結果しか期待できない.位置特定メソッドのオプションであるAccuracyGoalPrecisionGoalMaxIterationsは,根を見付けるときに許容誤差を制御するためにFindRootを使用するこれらの位置特定メソッドに渡すことができる.

ブール値のイベント関数では,関数がTrueからFalseへ,あるいはその逆に変化するときにイベントが起る."Direction"オプションを使うと,イベントがTrueからFalseへ変わるときのみ("Direction"->-1),あるいはFalseからTrueへ変わるときのみ("Direction"->1)に制限することができる.

以下は,クリックすることにより,変数stopの値が初期化された値のFalseからTrueに変更されるボタンを持つ小さいウィンドウを開く:
これは,ボタンがクリックされる時点まで(あるいはシステムのメモリがなくなるまで)振子方程式を積分する:

上記の例ではstopの即時型の値が評価され,関数として設定されないようにするために,RuleDelayed (:>)を使って"Event"オプションを指定している.

複数のイベントを指定することもできる.イベント関数を数値的評価してリストになった場合,リストの各要素は別のイベントとみなされる.各イベントについて異なるアクション,方向等を指定したい場合は,そのオプションの値を,適切な長さのリストとして指定するとよい.

これは,ボタンがクリックされる時点まで振子方程式を積分する.振子の完全なスイング数は,積分の間記録されている:

上記の例から分かるように,実数値とブール値のイベント関数を共用することは可能である.想定される要素数と各要素のタイプは,初期条件での値に基づいており,積分の間,一貫している必要がある.

"EventLocator""EventCondition"オプションを使うと,イベントのテストで満足しなければならない追加のブール条件を指定することができる.ブールイベントの代りにこれを使った方が根の探索がより効率的に行われるため,有利である.

以下は,減衰によりエネルギー積分がに減少したとき最初に となるところで減衰振り子の積分を中止させる:
以下で解(黒),導関数(青),エネルギー積分(緑)のプロットを作成する.エネルギーの閾値は赤で示してある:

"EventLocator"Methodオプションを使うと,積分で使用する数値メソッドの指定が可能である.

イベントの位置特定メソッド

"EventLocator"メソッドは,内在するメソッドのステップを取り,イベント関数のいずれかの符号(パリティ)がステップの終点で異なるかどうかをチェックすることにより動作する.イベント関数には,実数値かブール値が想定されているため,符号の変化があれば,そのステップ区間でイベントがあったことになる.イベント関数で,1ステップ内にイベントを含むものについては,その区間内におけるイベントの位置を特定するために,絞込み操作が行われる.

位置の絞込みには,いくつかの異なるメソッドがある.このようなメソッドには,積分区間の最初と最後で単純に解を取るというもの,イベント値の線形補間,小区間内の根の探索メソッドの使用等がある.適切なメソッドは,実行速度か位置の正確さのどちらを犠牲にするかによる.

イベントアクションが積分の中止となっている場合は,その積分が中止される特定の値は,"EventLocator"EventLocationMethodオプションで得られる値に依存する.

1つのイベントの位置は通常迅速に見付かるので,使用するメソッドが全体の計算時間に大きく影響を与えることはない.しかし,イベントが複数回検出されると,位置の絞込みメソッドにより影響が出る可能性がある.

"StepBegin"メソッドと"StepEnd"メソッド

イベントの正確な位置は大して問題ではない場合,あるいはイベントの正確な位置が分かっても,内在する計算に確度を持ったものが何も反映されない場合は,最も大雑把なメソッドが適している.前セクションのStopボタンの例がそれである.時間ステップは非常に迅速に計算されるので,一定の時間ステップ内でボタンをクリックするのにかかる時間を計測する方法はない.ましてや1つの時間ステップ内の特定の点における時間等測ることはできない.従って,継承されたイベントの精度に基づいている場合は,絞込みを行う意味がない.これは"StepBegin"あるいは"StepEnd"位置特定メソッドを使って指定することができる.イベントの定義がヒューリスティックであったり,幾分不明瞭であったりした場合は,このメソッドが適切である.

"LinearInterpolation"メソッド

グラフにプロットする点のためにイベントの結果が必要な場合は,グラフの解像度に合わせてイベントの位置を見付けるだけでよい.これにはステップの終点を使うだけでは通常大雑把すぎるが,イベント関数の値に基づく単独の線形補間で十分である.

数値積分の連続したメッシュ点におけるイベント関数を以下のように表す.

線形補間で次が得られる.

よって,イベント時間の線形近似は次のようになる.

線形補間は,イベント時間における解の近似にも使うことができる.しかし,導関数値 および がメッシュ点で利用可能なので,イベントにおける解のよりよい近似は,次のように三次エルミート補間を使って安価に計算することができる.

ここで適切に定義された補間の重みは以下の通りである.

設定を"LinearInterpolation"とすると,単独の線形補間に基づく絞込みが指定できる.

これは,振子方程式の1周期の解を計算し,その周期の解をプロットする:

周期全体に渡るプロットの解像度において,導関数が軸と交わる正確な点が終点にはならないことは分からない.しかし,十分ズームすると,誤差が見える.

以下は,ちょうど終点付近のプロットである:

線形補間法は,以下のセクションで示すポアンカレ断面の例のように,見ることを目的とするほとんどの場合には十分な方法である.しかし,ブール値を取るイベント関数では,線形補間は実質的には2分法の1ステップでしかないので,線形補間法はグラフィックスには不適切である場合がある.

ブレント法(Brent's Method)

デフォルトの位置特定メソッドは,イベント位置特定メソッドの"Brent"である.これはFindRootブレント法を使ってイベントの位置を見付けるものである.ブレント法は小区間の根から開始して,補間法と2分法に基づくステップを組み合せることによって,少なくとも2分法と同じくらいの収束率を保証する.FindRootがイベント関数の根をどの程度の確度・精度で見付けるようにするかは,イベント位置特定メソッドのオプション"Brent"を使って制御することができる.デフォルトは,NDSolveが局所誤差制御で使っているのと同じ確度・精度で根を見付けることである.

密な出力や難解な出力をサポートするメソッドでは,イベント関数の引数は,連続出力公式を使うことにより非常に効率よく求めることができる.しかし,連続出力をサポートしないメソッドでは,内在するメソッドのステップを取ることで解を計算する必要があるため,比較的高価になる可能性がある.解の近似を求める別の方法で,メソッドの順序通りではないが,NDSolveから返されるInterpolatingFunctionオブジェクトにFindRootを使うことに見合ったものとして,三次エルミート(Hermite)補間がある.これは,ステップの終点における解の値と解の導関数値に基づく補間によって,ステップの中盤で解の近似値を得るものである.

比較

以下の例では,多くの異なるイベント位置特定メソッドで振り子方程式を積分し,イベントが見付かった時間を比較する.

以下で,使用するイベント位置特定メソッドを定義する:
系を積分する.積分が終了すると,使用したメソッドに加え,独立変数の値およびイベントの値を出力する:

例題

落下する物体

以下の系は空気抵抗を受けた重力のもとで落下する物体をモデル化する([M04]を参照).

イベントアクションは,落体が地面に接した時間を保存し,積分を中止することである:
解をプロットして,最初と最後の点をそれぞれ緑,赤で囲んでハイライトする

ファン・デル・ポル(van der Pol)発振器の周期

ファン・デル・ポル発振器は,極度に硬い常微分方程式系の例である.イベントの位置特定メソッドは,常微分方程式系の積分を実際に行うのにどのようなメソッドを呼び出すこともできる.デフォルトメソッドのAutomaticは必要に応じて硬い系に適したメソッドに切り替えるので,硬さによる問題はない.

以下は,変数 がその初期値と方向に到達する点まで,パラメータ の特定値に対するファン・デル・ポル系を積分する:

初期条件でのイベントは考慮されていない点に注意する.

NDSolveの解の終点を選ぶことにより,周期を返す関数を の関数として書くことができる.

以下は,周期を返す関数を の関数として定義する:
における周期を計算する関数を使用する:

もちろん,周期解を使ってこのメソッドをすべての系に一般化することは簡単である.

ポアンカレ断面

ポアンカレ断面は,高次元の微分方程式系の解を視覚化するのに便利である.

インタラクティブなグラフィカルインターフェースについては,EquationTrekkerパッケージをご参照いただきたい.

ヘノンハイレス(HénonHeiles)系

ここでは銀河系の恒星運動をモデル化するヘノンハイレス系を定義する.

NDSolveProblemsパッケージからヘノンハイレス系を呼び出す:

この場合,関心のあるポアンカレ断面は,軌道が を通過するときの 平面の点の集合である.

数値積分の実際の結果は必要ではないので,出力変数リスト(NDSolveの第2引数)を空,つまり{}と指定することにより,InterpolatingFunctionにすべてのデータを保存することを避けることができる.こうすると,NDSolveが出力にInterpolatingFunctionを生成しないため,多量の不必要なデータを保存する必要がなくなるのである.NDSolveNDSolve::nooutというメッセージで出力関数がないという警告を出すが,この場合は必要なデータがイベントアクションから集められているので,メッセージをオフにしても差し支えない.

ここでの計算目的は,比較的低解像度のグラフで結果を見ることにあるため,線形補間によるイベント位置特定メソッドが使われる.ポアンカレ地図の固定点等,詳細を見たり,特徴を調べたりするためにグラフをズームする必要がある例題を行っている場合は,デフォルトの位置特定メソッドを使った方が適切であろう.

出力がないという警告メッセージをオフにする:
固定刻み幅0.25で四次の陽的ルンゲ・クッタ法を使い,ヘノンハイレス系を積分する.イベントアクションは の値にSowを使うことである:
ポアンカレ断面をプロットする.集まったデータは,Reapの結果の最後の部分であり,点のリストはその最初の部分である:

ヘノン-ハイレス系はハミルトニアン(Hamiltonian)なので,この例ではシンプレクティック法を使った方が良質な結果が得られる.

刻み幅0.25で四次のシンプレクティック分割ルンゲ・クッタ法を使い,ヘノン-ハイレス系を積分する.イベントアクションは の値をSowすることである:
ポアンカレ断面をプロットする.集まったデータは,Reapの結果の最後の部分であり,点のリストはその最初の部分である:

ABCフロー

三次元オイラー方程式の層流のカオスをモデル化するために使われるABC(ArnoldBeltramiChildress)フローの例題をロードする:
右辺の要素のいくつかを零に設定して,系の分割したY1Y2を定義する:
陰的中点法を定義する:
体積と逆対称を維持する二次の分割法を構築する:
特定の初期条件に対するポアンカレ断面を与える関数を定義する:
以下で,いくつかの異なる初期条件に対するポアンカレ断面を見付け,それらすべてを1つの点のリストに平坦化する:

バウンドするボール

この例は,[SGT03]の例題の一般化である.これは,指定された特徴で斜面をバンドしながら降りてゆくボールをモデル化する.この例題は,複数のNDSolveを位置特定メソッドで使って,ある動作をモデル化するよい例である.

以下は,1つのバウンドから次のバウンドへの解を計算する関数を定義する.解は,軌道が次に斜面に接するまで計算される:
ボールが斜面にぶつかったときのバウンドの関数を定義する.この公式は,バウンドした後,エネルギーはほんの しか残っていないと仮定し,斜面に対する垂線の投影に基づいている:
指定された投影率,斜面,初期位置でのバウンドするボールのシミュレーションを実行する関数を定義する:
[SGT03]で実行されている例題である:
以下で斜面を四半円で定義する:
斜面にわずかな波形を加える:

イベントの向き

常微分方程式

この例では,4つの方程式の標準的な硬くないテストシステムである制限三体問題の解を示す.この例は,月の周囲を回って地球に戻る宇宙船の軌道を追跡する([SG75]の246ページを参照).複数イベントおよび零交差の方向が指定できるということが重要である.

軌道を周期的にするために,初期条件を選ぶ. の値は月と地球の周囲を回る宇宙船に対応している:
イベント関数は,初期条件からの距離の導関数である.値が零と交差するときに極大値・極小値となる:
2つのイベントがあり,ここでは同じものである.Direction 1の最初のイベントは,初期値からの距離が極小値となる点に当たるので,宇宙船はもとの位置に戻る.イベントアクションは,イベントの時間を変数tfinalで保存し,積分を中止することである.2つ目のイベントは,極大値に当たる.イベントアクションは宇宙船が最初の位置から一番遠いところにある時の時間を変数tfarで保存することである:

最初の2つの解の要素は,無限小量の物体の座標である.従って,ひとつの要素に対して別の要素をプロットすると,物体の軌道が得られる.

以下は,宇宙船が最初の位置から最も遠い地点にあるときの軌道の半分を示している:
宇宙船がもとの位置に戻ったときの完全な軌道を示したものである:

遅延微分方程式

以下の系は感染病をモデル化したものである([HNW93],[ST00],[ST01]を参照のこと):
積分が進むにつれての各要素の極大値のデータを収集する.要素を分離するためにSowおよびReapに対して別々のタグが使われる:
極大値を解の要素とともに表示する:

連続方程式とスイッチ関数

多くのアプリケーションでは,微分方程式系の関数がどこでも解析関数,あるいは連続関数とは限らない場合がある

実践で生じる一般的な不連続問題には,次のようにスイッチ関数の が関与している.

不連続性を渡る難しさを示すために,以下の例[GØ84]を考えてみる([HNW93]も参照のこと).

これは,系全体の入力である.シンボルeventにスイッチ関数が割り当てられており,系を定義する関数は,スイッチ関数の符号に依存している:
積分に使用する数値メソッドを示すために,シンボルodemethodが使われる.これと比較する目的で,"ExplicitRungeKutta"等の異なるメソッドを定義し,このセクションの計算を評価して,他のメソッドがどのように動作するかを見てみるとよい:
以下は区間[0, 1]上で系を解き,ReapSowを使って積分のメッシュ点のデータを集める:
解のプロットである:

解が実際に得られたが,効率的に得られたのかどうかは明らかではない.

以下の例は,不連続性の通過により数値ソルバの使用が困難になることを示している.

これは,積分のメッシュ点と積分ステップの数を表示する関数を定義する:
積分が不連続性を通過するとき(値で0.6付近),積分メソッドは困難に陥り,多くの小さいステップが取られる.却下されたステップが多数見られることもある:

不連続性を通過する最も有効なメソッドのひとつに,積分を中止して不連続点で再開するというものがある.

以下の例では,これを行うためにどのように"EventLocator"メソッドを使えばよいかを例示する.

これは,不連続点まで系の最初の部分を数値積分する.スイッチ関数はイベントとして与えられる.イベントの方向は,負から正への方向に限られている.イベントが見付かると,解とイベントの時間がイベントアクションによって保存される:

"EventLocator"メソッドで見付かった不連続性を新しい初期条件として使うと,連続積分できるようになる.

これは系と初期条件を定義し,その系を数値的に解き,メッシュ点に使われるデータを集める:
2つの解のプロットは,系全体を1度に解いて得られたものとよく似ている:
メッシュ点を調べてみると,このメソッドで取られたステップ数の方が格段に少なく,不連続付近で見られた問題動作がなくなっていることが分かる:

不連続の値は[HNW93]では0.6234で与えてある.これは"EventLocator"メソッドで見付かった値と偶然にも一致している.

以下の例では,系を解析的に解き,数値法を使ってその値をチェックすることが可能である.

不連続までの系の解はベッセル関数とガンマ関数で表すことができる:
スイッチ関数に解を代入すると,局所的最小化により不連続性の値が確認される:

偏微分方程式での折返しの回避

多くの発展方程式は,特化された離散化メソッドを使わずに全領域を離散化することができなくなるほどの,無限もしくは十分大きい空間領域上で動作する.実際には,興味のある解が局所化されている限り,小さい計算領域を使うことができる場合もある.

計算領域の境界が,研究中の実際のモデルではなく実践的であるかどうかによって決められている場合,境界条件を適切に選ぶことができる.擬スペクトル法を周期境界条件で使用すると,その優れた近似結果のため,計算領域の範囲を拡張することが可能である.周期境界条件の欠点は,境界を超えて伝播する信号が領域の反対側に残り,折返しによる解への影響が出る.この影響を最小限にするために,境界付近に吸収層を使うこともできるが,常に完全に排除できるというわけではない.

サインゴルドン(Sin-Gordon)方程式は,微分幾何学および相対論的場の理論で現れる.以下の例では,散らばっている局所化された初期条件から始めて,サインゴルドン方程式を積分する.積分には,周期的な擬スペクトル法が使われる.境界付近には吸収層を置いていないので,折返しが顕著になった時点で積分を中止するのが最も妥当である.この条件は"EventLocator"メソッドを使ったイベント位置特定で,簡単に検出できる.

積分は,周期的な折返し点における解の大きさが,閾値0.01を超えると終了する.そこを超えると,周期性による影響が波形に及ぶ:
以下は,InterpolatingFunctionオブジェクトの終了時間を抽出し,計算された解のプロットを作成する.最初の波形が境界にちょうど達したところで,積分が終了していることが分かる:
"DiscretizedMonitorVariables"オプションは,イベントが偏微分方程式でどのように解釈されるかを指定する.Trueに設定すると,u[t,x]は離散化された値のベクトルで置き換えられる.これは,イベントを評価するために明示的にInterpolatingFunctionを構築することを避けるため,ずっと効率的である:

パフォーマンスの比較

以下の例では,2つの異なる積分法を比較する表を作成する.

以下で,やや減衰した振子方程式が一時的に1000回止まった点までの解を計算するのにかかった時間を返す関数を定義する:
この関数を使って,2つの異なる常微分方程式積分法に対する異なる位置特定メソッドを比較する表を作成する:

単純なステップ開始/終了法と線形補間法は同様にコストが低いが,これよりもよい位置特定メソッドになると高価になる.デフォルトの位置特定メソッドは,連続出力公式をサポートしていないので,陽的ルンゲ・クッタ法では特に高価になる.従って,デフォルトメソッドは,局所的な最小化の過程で,異なる刻み幅で繰返しメソッドを起動される必要がある.

ここで特筆すべきことは,イベントの計算にかかる付加的な時間の大分部は,符号の変化の可能性をチェックするために,各時間ステップでイベント関数を評価する必要があるということにより生じる.

最適化は,独立変数のみを含むイベント関数について実行される.そのようなイベントは,初期化時に自動的に検出される.このことには,例えば,従属変数の解の補間が各ステップで実行されないという利点がある.局所的最適化検索の場合は,独立変数の値が見付かるまで遅延されるのである.

限界

イベント関数は,ステップ区間での符号の変化しかチェックされないため,あるステップ区間に複数の根がある場合,すべてあるいはいくつかのイベントが見失われる可能性がある.これがイベント位置特定メソッドのひとつの限界である.これは一般に,常微分方程式の解がイベント関数より格段に遅く変化するときにのみ起る.これが起ったと思われる場合,NDSolveMaxStepSizeオプションを使って,メソッドが取れる最大刻み幅を削減するというものが最も簡単な解決策である.高度な方法もあるが,最善策は何を計算しているかによる.以下の例では,この問題とそれに対する2種類の解決策を示す.

これは より小さい正の整数の数(148個ある)を計算する.しかし,解x[t]が非常に単純なため,メソッドが大きい時間ステップを取り,それによりほとんどが見失われている:
すべてのイベントが見付かるよう,最大刻み幅を制限する:

この例から,終点が大きくなるにつれ,最大刻み幅を小さくしていく必要があるという可能性が高いということは明らかである.いたるところで非常に小さい刻み幅を取るのは,極めて非効率的である.イベント関数と同じ時間スケールで変化する変数を設定することにより,適応型の時間刻み制限を導入することが可能である.

以下で,積分すべきイベント関数である追加関数を導入する.以下のように修正して,メソッドで必要なだけのステップが取れるようにすると,までの正確な値を妥当な時間で見付けることができる:

オプションの要約

"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"線形補間を用いてイベント時間を特定する.この後,イベント時間での会を見付けるために,三次エルミート補間が使われる
"StepBegin"イベントはステップの最初の解で与えられる
"StepEnd"イベントはステップの最後の解で与えられる

"EventLocationMethod"オプションの設定

"Brent"オプション

オプション名
デフォルト値
"MaxIterations"100メソッドの1ステップ内でイベントを見付けるために使う最大反復回数
"AccuracyGoal"AutomaticFindRootに渡される目標確度設定.Automaticのとき,FindRootに渡される値は,NDSolveの局所誤差設定に基づく
"PrecisionGoal"AutomaticFindRootに渡される目標精度設定.Automaticのとき,FindRootに渡される値は,NDSolveの局所誤差設定に基づく
"SolutionApproximation"Automatic絞込みの過程でイベント関数を評価するための解を近似する方法.Automaticあるいは"CubicHermiteInterpolation"のいずれかを取る

イベント位置特定メソッド"Brent"法のオプション