NDSolveのIDAメソッド
IDAパッケージは米国ローレンスリバモア(Lawrence Livermore)国立研究所のCASC(Center for Applied Scientific Computing)で開発されたSUNDIALS (SUite of Nonlinear and DIfferential/ALgebraic equation Solvers)の一部である.IDAのユーザガイド[HT99]には,「IDAは微分代数方程式(DAE)系の初期値問題のための汎用ソルバである.IDAとう名称はImplicit Differential-Algebraicソルバの頭文字を取ったものである.IDAはDASPKに基づいており...」と記述されている.DASPK [BHP94],[1]は大規模の微分代数系を解くためのFortranコードである.
Wolfram言語ではIDAパッケージにインターフェースが提供されているので,残差を評価したりプログラムをコンパイルしたりするのにC言語で関数を書くのではなく,NDSolveに入力した方程式からWolfram言語が自動的に関数を生成する.
微分代数方程式の「指数」とは,常微分方程式系を得るために微分代数方程式を微分するのに必要な回数のことである.微分代数方程式およびその起源に関する情報は,微分代数方程式の数値解法に記載されている.
IDAは,可変刻み形式で実装された一次から五次までの後退微分式(BDF)で系(2)を解く.時間 のときの 次のBDFは以下の公式で与えられる.
係数 は次数 と過去の刻み幅によって決まる.BDFをDAE (3)に適用すると,以下の非線形方程式系が与えられる.
この系の解はニュートン型の方法を使って求められる.一般的にヤコビアンの近似が使われる.
「IDAの特筆すべき機能は,各時間刻みに内在する非線形系の解法として,ニュートン/直接法あるいは非厳密ニュートン/クリロフ(反復)法のどちらかが選べることである.」[HT99]Wolfram言語では,メソッドのオプションを使ってこれらのソルバにアクセスすることもできれば,大規模問題では自動的に疎行列直接解法に切り換えるWolfram言語のデフォルトのLinearSolveを使うこともできる.
解法の各ステップで,IDAは局所打切り誤差の推定量 を計算する.刻み幅と次数は,重み付きノルムNorm[En/wn]が1未満となるよう選ばれる. の 番目の要素 は以下のように与えられる.
値 prec と acc はNDSolveのPrecisionGoal->prec とAccuracyGoal->acc の設定から取られる.
IDAは特に非線形方程式の解法において非常に柔軟性に富んでいるので,解法の制御を可能にするオプションが多数ある.NDSolveにMethod->{IDA,ida method options}というオプションを設定することで,IDAのメソッドオプションが使用ができる.
NDSolveから返されたInterpolatingFunctionオブジェクトで計算された中間値の厳密な確度が重要なときは,NDSolveメソッドオプション設定のInterpolationOrder->Allを使うとよい.この設定では,時間ステップ間の解を表すために,「密な出力」と呼ばれることもあるメソッドの次数に基づいた補間が使用される.デフォルトで,NDSolveは視覚的に満足できる解を表すために最小限のデータを保存する.データ量を最小限に保つことは,より複雑な解の場合に必要なメモリと時間を節約することになる.
最小限の出力とメソッドの完全な補間次数の差が顕著に現れる例として,厳密解がDSolveで計算できる簡単な線形方程式の解から導いた数量 を追跡してみよう.
厳密解は,密な出力で,あるいは密な出力なしで計算された解と比較される.
このプロットから,時間ステップが大きくなるにつれ,デフォルトの解の出力の誤差が時間ステップ間で格段に大きくなることが明確である.蜜な出力の解は,時間ステップの間でもソルバによって計算された解を表す.しかし,密な出力の解が取るスペースはずっと大きくなることを覚えておかねばならない.
解から導きたい数量が複雑な場合,その数量を代数的量として与えることで許容誤差内に局所的に収め,ソルバに強制的に誤差制御させるようにすることができる.
微分代数方程式の解は,関心数量の誤差を制御するためにずっと多くのステップを必要とするが,密な出力を使うよりもメモリ使用量はずっと少なくてすむ.
これから先は,"ImplicitSolver"オプションの2つの設定"Newton"と"GMRES"に焦点を当てて論を進める.設定を"Newton"にすると,ヤコビアンあるいはその近似が計算され,ニュートンステップを見付けるよう線形系が解かれる.一方,"GMRES"は反復非線形ソルバ法で,ヤコビアン全体を計算するのではなく各反復ステップの方向微分係数を計算するものである.
"Newton"法には"LinearSolveMethod"というオプションが1つあり,ニュートンステップを見付けるのに必要な線形方程式の解き方をWolfram言語に指示するために使うことができる.このオプションは以下の値を取ることができる.
Automatic | これがデフォルト.ヤコビアンのバンド幅に応じて,Wolfram言語のメソッドLinearSolveとBandとを自動的に切り替える.バンド幅の大きい系では,疎なヤコビアンを持つ大規模システムのための疎行列直接解法に自動的に切り替えられる |
"Band" | IDAバンド法を使う.詳細はIDAユーザマニュアルを参照のこと |
"Dense" | IDAデンス法を使う.詳細はIDAユーザマニュアルを参照のこと |
"LinearSolveMethod"オプションで可能な設定
"GMRES"法の方が実質的に高速であるが,通常使用する際に面倒なことがある.GMRESが実際に有効となるためには,オプションで提供される前処理が必要な場合が多いためである.クリロフ(Krylov)の部分空間処理を制御する他のメソッドオプションもある.これらの使い方については,IDAのユーザマニュアル[HT99]を参照するとよい.
GMRESメソッドオプション名 |
デフォルト値
| |
"Preconditioner" | Automatic | 前処理を実行する別の関数を返すWolfram言語関数 |
"OrthogonalizationType" | "ModifiedGramSchmidt" | "ClassicalGramSchmidt"とすることもできる.IDAユーザガイドで変数gstypeを参照のこと |
"MaxKrylovSubspaceDimension" | Automatic | 部分空間の最大次元.IDAユーザガイドで変数maxlを参照のこと |
"MaxKrylovRestarts" | Automatic | 再実行の最大回数.IDAユーザガイドで変数maxrsを参照のこと |
例題として,二次元のバーガーズ(Burgers)方程式を挙げる.
これは通常常微分方程式ソルバで解かれるが,この例ではDAEソルバを使って2つのことが行われる.まず境界条件が代数条件として実行される.次にNDSolveが代数項を使って保存差分法を使うよう強制的される.比較のため,初期条件,境界条件に既知の厳密解を使う.次の計算は,この2つのアプローチの違いを示す.
プロットから,のとき前面が非常に急な傾斜が見られる.これは一定の速度で移動する.
次では,IDA法の異なるオプション設定による差を比較する.オプション設定を強調するために,解の計算時間を測定して最大誤差を与える関数を定義する.
ここでは,ヤコビアンのバンド幅が比較的大きい四階導関数が使われている,境界付近で使われている片側のステンシルにより上と下に幅が加わっている等の理由から,"Band"法はあまり有効ではない.バンド幅は明示的に指定することができる.
解法の計算時間は短いが,誤差はわずかに大きく,時間刻みの合計数はかなり多くなっている.問題がさらに硬ければ,もう一方からの情報が欠けているため反復は収束しなかったであろう.理想的には,バンド幅は情報を次元のどこからも削除すべきではない.
バンド幅の設定をより適切なものにしても,解の計算にかかる時間はデフォルトのバンド幅の場合よりもわずかに速いが,デフォルトの線形ソルバを使うよりも遅い."Band"法は二次元問題で有効な場合もあるが,通常一次元問題で最も有効である.
前処理をせずに"GMRES"法を使うと,デフォルトメソッドと比較して,計算時間が長くなったりステップ数が多くなったりする.このような理由で,この方法は前処理なしではお勧めしない.しかし,よい前処理を見付けることは通常自明なことではない.この例では対角前処理を使う.
"Preconditioner"オプションの設定は,NDSolveによって与えられる4つの引数を取る関数 f でなければならない.その際,f[t,x,x′,c]が残差ベクトルに前処理を適用する別の関数を返す(前処理の使い方についての詳細はIDAユーザガイド[HT99]を参照のこと).引数 t,x,x′,c はそれぞれ現時間,解のべクトル,解の導関数べクトル,上記の公式(4)の定数 c である.例えば,これらの引数の関数として,適切な前処理行列 を生成する手続きを決めることができれば
"Preconditioner"->Function[{t,x,xp,c},LinearSolve[P[t,x,xp,c]]]
を使って,前処理行列 を実質的に反転させるLinearSolveFunctionオブジェクトが生成できる.通常前処理関数は設定される度に残差ベクトルに数回適用されるので,LinearSolveFunctionに含まれるような種類の因数分解を使うのはよいことである.
対角行列の場合は,単に相互行列を使うだけで反転が可能である.対角前処理を設定する上で最も難しい部分は,境界の値がそれによって影響を受けてはならないという点である.
未完成の前処理を使っても,"GMRES"法は疎行列直接解法を使うよりも速く解を計算する.
高階の一次微分を伴う偏微分方程式の離散化や偏微分方程式系では,対応するNDSolve`StateDataオブジェクトを見て,前処理の構造形式が正確に得られるように,変数の順序を決めなければならない.