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

はじめに

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

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

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

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

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

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

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

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

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

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

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

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

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

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

In[1]:=

In[2]:=

In[4]:=

In[6]:=

In[7]:=

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

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

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

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

NDSolveのメソッドの初期化

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

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

In[12]:=

In[13]:=

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

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

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

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

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

In[14]:=

In[15]:=

In[16]:=

In[17]:=

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

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

別のメソッドの開始

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

サブメソッドを使う

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

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

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

In[1]:=

In[2]:=

In[3]:=

In[4]:=

In[6]:=

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

アダムス法

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

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

ここでは差分商である．

In[29]:=

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

In[30]:=

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

In[32]:=

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

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

In[33]:=

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

In[34]:=

In[35]:=

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

In[36]:=

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

In[37]:=

In[38]:=

In[40]:=

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

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

In[43]:=

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

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

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