NDSolveのメソッドプラグインフレームワーク

はじめに別のメソッドの開始
古典的ルンゲ・クッタ(Runge–Kutta)法アダムス法
クランク・ニコルソン（Crank–Nicolson）法

はじめに

NDSolveに設定されている制御メカニズムにより，ユーザが数値積分アルゴリズムを自分で定義し，それをNDSolveのMethodオプションの指定として使うことができる．

NDSolveはその数値アルゴリズム，およびその情報にオブジェクト指向型でアクセスする．数値積分の各ステップで，NDSolveは必要に応じてプライベートデータが保存できる形式でメソッドを保存する．

AlgorithmIdentifier [data]常微分方程式の特殊な数値積分アルゴリズムで使う必要のある任意のデータ data を含むアルゴリズムオブジェクト．そのデータは実質的に，アルゴリズムには見えない．AlgorithmIdentifier はWolframシステムのシンボルでなければならない．NDSolveはオプションMethod->AlgorithmIdentifier を使って，このアルゴリズムにアクセスする．

NDSolveで使われるメソッドデータの構造

NDSolveはアルゴリズムに関連するデータに直接アクセスしないので，必要な情報を使いやすく効率的な形式で保存することができる．このプライベートデータに保存されるアルゴリズムと情報は，アルゴリズムオブジェクトのメソッド関数によってのみアクセスされる．メソッドアルゴリズムオブジェクトとNDSolve間の主なインタラクションは，"Step"関数を介して行われる．

obj["Step"[inputs]]数値アルゴリズムを使って1つの時間ステップを取ることを試みる

NDSolveによって使われる，必要な"Step"関数

次は，メソッドの重要な特性と，定義しなければならないその"Step"関数のリストである．

obj["StepInputs"]ステップ関数に与えなければならない入力引数のリスト
obj["StepOutput"]ステップ関数が返す出力のリスト
obj["StepMode"]ステップモードに設定する
obj["DifferenceOrder"]アルゴリズムの現在の漸近的差分次数に設定する

NDSolveでメソッドアルゴリズムオブジェクト obj を使うために定義しなければならない特性

下の表はステップ入力および／または出力に使うことのできる指定である．

簡約名
完全名    
要素
"F""Function"常微分方程式系の右辺あるいは微分代数方程式系の残差に与えられる関数
"T""Time"現在の時間
"H""TimeStep"時間刻み
"X""DependentVariables"離散変数以外の従属変数の現行値
"XI""DependentVariablesIncrement" がステップ後に従属変数の値を与えるような．が$Failedとして与えられたら，ステップは棄却され再計算される（返される時間刻みは小さくなければならない）
"XP""TemporalDerivatives"従属変数の時間導関数
"D""DiscreteVariables"値の連続範囲を取る離散変数
"ID""IndexedDiscreteVariables"値の有限集合を取る離散変数
"P""Parameters"感度が計算されているパラメータ以外のパラメータ
"SP""SensitivityParameters"感度が計算されているパラメータ
"SD""SolutionData"解データのリスト
"Obj""MethodObject"次のステップで使われる obj の新しい設定
"State""StateData"NDSolve`StateDataオブジェクト

ステップ入力および出力の指定

解データは，さまざまな変数の型に関連する要素を持つ値のリストである．アクティブな方向の解データは，state@"SolutionData"["Active"]を使ってNDSolve`StateDataオブジェクトの state から得ることができる．上記の指定の多くは解データの要素に直接対応している．

入力の場合，関数の指定は，その関数に与えられる引数で与えられる．これらの引数は解データの要素でなければならない．AllはNDSolveによって使われる引数がすべて想定されていなければならないということを示すために使うことができる．

出力では，メソッドオブジェクトが各ステップを変更しない場合等．場合によってメソッドが値のいくつかしか返さないならば，すべての出力にNullを使うことができる．ただし，使ってはならないことを示す増分には使えない．

古典的ルンゲ・クッタ(Runge–Kutta)法

以下は，簡単な数値計算アルゴリズムを設定して，それにアクセスする方法の例である．

これは，四次の古典的ルンゲ・クッタを使って，常微分方程式の積分の1ステップを取るメソッド関数を定義する．このメソッドは簡単なので，プライベートデータを保存する必要はない：
以下はNDSolveのステップ関数の入力と出力を指定する．___テンプレートは，CRK4メソッドオブジェクトのすべてのインスタンスで利用できるようにするために使われている：
同じ指定を完全名を使って書き表してみる：
次はNDSolveがメソッドに使用するのに適切な差分次数を得ることができるようメソッド関数を定義する：
NDSolveに時間刻みの制御方法が分かるように，ステップモードに対するメソッド関数を定義する．このアルゴリズムメソッドには，刻み幅制御がないので，ステップモードをFixedと定義する：
固定刻み幅の簡単な調和振動方程式を積分する：

一般に固定刻み幅を使うことは，積分の局所的な問題に応じて刻み幅を変化させることよりも効率が悪い．現在の陽的ルンゲ・クッタ法（NDSolveでMethod->"ExplicitRungeKutta"としてアクセスできる）には，いわゆる組込みの誤差推定法が含まれているため，適切な刻み幅を非常に効率的に決定することができる．別の方法として，外挿を使う組込みの刻み幅制御法を使うこともできる．"DoubleStep"メソッドは，刻み幅 の1ステップと刻み幅 の2ステップによる時間刻みの積分に基づく外挿を使う．メソッド"Extrapolation"はより高度な外挿を行い，積分が実行されるに従って自動的に外挿の程度を変更する．しかし，このメソッドは，通常，差分次数が1および2の基本メソッドで使われる．

"DoubleStep"メソッドで制御されているステップで四次の古典的ルンゲ・クッタを使って，簡単な調和振動を積分する：
計算された解の誤差をステップの最後で比較するプロットを作成する．"DoubleStep"メソッドの誤差は，青で示される：

全体的に見て，固定刻み幅は可変刻み幅より誤差が小さくなっている．これは，固定刻み幅のステップの方が格段に小さく，3倍以上のステップ数を必要とするためである．局所解の構造が顕著に変化するような問題の場合は，その差はさらに大きいものになり得る．

硬さ検出の機能については「NDSolveのDoubleStep法」に記載されている．

より高度なメソッドでは，メソッドで使用するデータを設定することが必要であったり，その方が効率的であったりすることがある．NDSolveがある数値計算アルゴリズムを初めて使うとき，初期化関数を呼び出す．自分のメソッドに適切なデータを設定する初期化規則を定義することができる．

NDSolve`InitializeMethod[Algorithm Identifier,stepmode,sd,f,state,Algorithm Options]
NDSolveが特定の積分アルゴリズムを初めて使うときに，初期化のために評価する式．stepmode は，メソッドが刻み幅制御器あるいは別のメソッドのフレームワーク内から呼び出されるかどうかによって，AutomaticかFixedのいずれかになる．sd は現行の解データで，state はNDSolveが使うNDSolveStateオブジェクトである．Algorithm Options はあるアルゴリズムを使う指定で特に与えられるすべてのオプションを含むリストである．例えばMethod->{Algorithm Identifier,opts}の{opts}がこれに当たる
Algorithm Identifier/:NDSolve`InitializeMethod[Algorithm Identifier,stepmode_,sd_,f_NumericalFunction,state_NDSolveState,{opts___?OptionQ}]:=initialization
規則がアルゴリズムと関連付けられるような初期化の定義．initialization はアルゴリズムオブジェクトを Algorithm Identifier[data]という形式で返す

NDSolveのメソッドの初期化

NDSolve`InitializeMethodはプロテクトされているため，これに規則を加えるためには，まずプロテクトを解除する必要がある．規則はメソッドに関連付けられたままにしておいた方がよい．これを整然と行うためには，上記のようにTagSetを使って初期化定義を作る．

例として，上に示したルンゲ・クッタ法を再定義して，厳密な係数2，1/2，1/6を使う代りに適切な精度の数値を使い，計算を少し速くしたいとする．

必要な係数に保存されている数値を使い，四次の古典的ルンゲ・クッタを使って，常微分方程式の積分の1ステップを取るメソッド関数を定義する：
これは後で使うデータを含むアルゴリズムオブジェクトを初期化する規則を定義する：

数の数値を保存すると，"DoubleStep"を使った長い積分で5％から10％のスピードアップが可能である．

クランク・ニコルソン（Crank–Nicolson）法

これは陰解法を設定する方法の例である．常微分方程式系 を解くためには，時間刻み幅 のスキームは である．ここで であり である．一般に非線形の のとき，方程式はニュートンの反復法で解く必要がある．

各ステップでヤコビアンが更新される本物のニュートン反復法を使う代りに，ヤコビアンが1度だけ計算され，そのLU分解を繰返し使用して線形方程式を解く．時間ステップが大きすぎない限り，この方法は通常完全なニュートン法と同様によく収束する．また，2回目以降の各反復は，実質的に安価になる．

ニュートン反復法が収束しないという問題を避けるために，このメソッドは最大反復回数と収束誤差を制御するメソッドオプションを取るよう設定されている．

クランク・ニコルソン法のオプションとメソッドの初期化を定義する：
ステップ関数を定義する：
NDSolveのクランク・ニコルソン法の入力，出力，差分次数，ステップモードを指定する：
オプション設定でメソッドを呼び出すことができる方法である：

線形の問題では，解は最初の反復で与えられる．しかし，非線形問題では，収束にもっと時間がかかる．

非線形問題は，このような大きい時間ステップで1回の反復で収束することは想定できないため，以下は失敗する：
反復回数を増やすと，解くのに成功する：

クランク・ニコルソン法は通常偏微分方程式を解く場合に考慮されるので，以下は波動方程式を解く方法の例である．これは線形方程式なので，1回の反復で収束し，メソッドは非常に速い．

クランク・ニコルソン法を使って，周期境界条件を持つ波動方程式を解く：
解の密度プロットである：

別のメソッドの開始

NDSolveフレームワークの強みの一つは，いろいろな方法でメソッドをネストすることができるという点にある．

NDSolve`InitializeSubmethod[caller,submspec,stepmode,sd,f,state]メソッド caller の初期化において，Method->submspec によって指定されたサブメソッドを初期化する．stepmode，sd，f，state はNDSolve`InitializeMethodに含まれるものと同じである．
NDSolve`InvokeMethod[submobj,f,h,sd,state]関数 f，時間ステップ h，解データ sd，NDSolve`StateDataオブジェクト state を持つメソッドオブジェクト submobj でサブメソッドを開始して，リスト{hnew,sdnew,submobjnew}を返す．ここで hnew は新しい時間ステップ，sdnew は変化しない要素をNullと設定したステップの終点における解データ，submobjnew は新しいサブメソッドオブジェクトである．

サブメソッドを使う

NDSolve`InitializeSubmethodは実質上，与えられたすべてのサブオプションからメソッド名を分離した後にNDSolve`InitializeMethodを呼ぶ．成功したら，後でNDSolve`InvokeMethodで使用できるメソッドオブジェクトを返す．

NDSolve`InvokeMethodは実質的に，サブメソッドに対する"Step"関数を構築する．関数の引数はサブメソッドの"StepInput"指定により指定されるので，NDSolve`InvokeMethodに与えられる関数 f はNDSolveで使われる引数すべてを持ち，"Function"[All]による呼出し"StepInput"指定で指定することができる．組込みメソッドでは，明示的な"Step"関数はない場合があり，NDSolve`InvokeMethodが有効な内部コードを使ってそれを操作する．NDSolve`InvokeMethodによって返される結果は，{"TimeStep","SolutionData","MethodObject"}の"StepOutput"に適合するよう，常にサブメソッドの実際の出力から構築される．

次は，ステップと差分次数の簡単な監視を提供するサブメソッドを呼び出す例である．

メソッドのメソッドオプションと初期化を定義する：
メソッドのステップ関数は単純にサブメソッドを起動し，結果と差分次数に監視関数を適用する：
入力と出力を設定する：
サブメソッドから他のすべての特性を引き出す：
メソッドを使う：

より複雑な積分の場合，このメソッドはデータを収集するためにReapと一緒に使うことができる．

ファン・デル・ポール振動子を積分するための時間の関数として，刻み幅と差分次数の対数を収集し，解をプロットする：
急激な変化があった辺りのステップの刻み幅と差分の次数の対数を示す：

アダムス法

NDSolveフレームワークでは，自動的に刻み幅を制御するアルゴリズムを書くことはたいして難しいことではない．しかし，数学的，プログラミング的見地からすると，誤差推定および適切な刻み幅の決定に必要とされる事項があるため，格段に難しいこととなっている．以下の例は，ShampineとWattsのFortran DEABMコードを部分的に適用して，NDSolveフレームワークにフィットさせたものである．アルゴリズムは，[SG75]に記述されている基準に基づいて，刻み幅と次数の両方を適宜選択する．

第1段階は，係数を定義することである．この積分法は，可変刻み幅係数を用いる．刻み幅のシーケンス （ここで は現在のステップ）があるとすると，次数 のAdams–Bashforth予測子と次数 のAdams–Moulton修正子を持つメソッドの係数 は，以下が成り立つ．

ここでは差分商である．

誤差推定で使われる の他，係数 と を計算する関数を定義する．この公式は[HNW93]からのもので，原則的に同じ表記法を使っている：

hlist は過去のステップからの刻み幅のリスト である．定数係数のアダムス係数は1度，さらに簡単に計算することができる．固定刻み幅のAdams–Moulton係数は，メソッドの次数の変更に対する誤差予測で使われるので，値を保存する規則で1度定義するということは理にかなっている．

固定刻み幅のAdams–Moulton係数の値を計算し，保存する関数を定義する：

次にステップ間で必要な情報を維持するデータ構造を設定して，そのデータをどのように初期化するかを定義する．保存しなければならない主要な情報は，過去の刻み幅のリスト hlist および差分商である．このメソッドは誤差推定を行うので，NDSolve`StateDataオブジェクトから正しいノルムを得る必要がある．精度等，この他のデータは，初期化と便宜のために保存される．

NDSolveからAdamsBMメソッドを初期化する規則を設定する：

初期化時にはステップが1つもないので，hlist は{}に初期化される．第0差分商は関数値なので，は初期条件の解のベクトルの導関数に初期化される．は行列であることに注意する．第3要素は，次のステップの最適次数を決定するために使われるもので，零ベクトルに初期化される．これは実質的には付加的な差分商である．データの他の部分の用途は，ステップ関数の定義で明確にされる．

初期化は，使用するメソッドの最大次数を制限するために使えるオプション値を得るようにも設定できる．コードDEABMでは，次数は12に制限される．これが機械精度計算の実際の限界であるためである．しかし，Wolframシステムでは高精度で計算され，高次のメソッドには非常に有利である．そのためこの種の絶対的な限界は必要ない．従ってオプションのデフォルトはと設定する．

AdamsBMメソッドのMaxDifferenceOrderオプションのデフォルトを設定する：

より高度な"Step"メソッド関数に進む前に，簡単なメソッド関数"StepMode"と"DifferenceOrder"を定義する．

AdamsBMメソッドのステップモードを常にAutomaticを定義する．つまり，このメソッドは，固定刻み幅に可変刻み幅を要求する可能性のある制御器からは呼び出せないということである：
AdamsBMメソッドの差分次数を定義する．これは保存されている過去の値の数により異なる：

最後に"Step"関数を定義する．実際のステップの過程は，ほんの数行である．この他のコードは，ShampineとWattsのDEABMコードに従って，自動刻み幅および次数選択を扱う．

前述のテンプレートに従ったステップデータを返す，AdamsBMの"Step"メソッド関数を定義する：

ここのコードにはDEABMから逸脱したものが少し含まれている．最も顕著なのは，DEABMがアップデートが必要な係数だけしか計算しないのに対して，このコードでは各ステップで係数が再計算されていることである．これはコードの簡潔さを維持するための変更であるが，特に小型から中型程度のシステムでは，明らかにパフォーマンスが劣る．2つ目の変更点は，拒絶するステップを制限するコードのほとんどはNDSolveに任せられているため，このコードには刻み幅が小さすぎるか，あるいは許容誤差が大きすぎるかを調べる術がないという点である．硬さ検出のヒューリスティックも欠けている．次数と刻み幅の決定コードはもモジュール方式で別の関数になっている．

，， で誤差推定 を構築し，低次数化されるべきか否かを決定する関数を定義する
ステップが成功した後で，最適の次数を決定する関数を定義する：
刻み幅 のステップが成功した後で，最適の刻み幅を決定する関数を定義する：

これらの定義を入力したら，Method->AdamsBMを使うだけでNDSolveのメソッドにアクセスできる．

既に定義したアダムス法で調和振動方程式を解く：
計算された解の誤差を表示する．誤差は妥当な範囲内にあることが分かる．最初の数点の後，刻み幅が増大していることに注目のこと：

厳密な許容誤差で高精度の計算を行うと，このメソッドは組込みメソッドのいくつかを超えるパフォーマンスを見せる可能性がある．これは，組込みのメソッドが次数12の制限付きのコードから適用されているためである．

係数の計算には長時間かかる．

これは，予想されたほど高い次数を使っていない．

いかなる場合でも，かかる時間の半分は係数の生成に費やされるので，係数のアップデートを把握する必要がある．