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は次の形式の指数1の微分代数方程式を解く．

微分代数方程式の「指数」とは，常微分方程式系を得るために微分代数方程式を微分するのに必要な回数のことである．微分代数方程式およびその起源に関する情報は，微分代数方程式の数値解法に記載されている．

IDAは，可変刻み形式で実装された一次から五次までの後退微分式(BDF)で系(2)を解く．時間のときの次のBDFは以下の公式で与えられる．

係数は次数と過去の刻み幅によって決まる．BDFをDAE (3)に適用すると，以下の非線形方程式系が与えられる．

この系の解はニュートン型の方法を使って求められる．一般的にヤコビアンの近似が使われる．

「IDAの特筆すべき機能は，各時間刻みに内在する非線形系の解法として，ニュートン／直接法あるいは非厳密ニュートン／クリロフ（反復）法のどちらかが選べることである．」[HT99]Wolfram言語では，メソッドのオプションを使ってこれらのソルバにアクセスすることもできれば，大規模問題では自動的に疎行列直接解法に切り換えるWolfram言語のデフォルトのLinearSolveを使うこともできる．

解法の各ステップで，IDAは局所打切り誤差の推定量を計算する．刻み幅と次数は，重み付きノルムNorm[E_n/w_n]が1未満となるよう選ばれる．の番目の要素は以下のように与えられる．

値 prec と acc はNDSolveのPrecisionGoal->prec とAccuracyGoal->acc の設定から取られる．

IDAは特に非線形方程式の解法において非常に柔軟性に富んでいるので，解法の制御を可能にするオプションが多数ある．NDSolveにMethod->{IDA,ida method options}というオプションを設定することで，IDAのメソッドオプションが使用ができる．

IDA法のオプションはNDSolve`コンテキストのシンボルIDAと関連付けられている：

Wolfram Language code: Options[NDSolve`IDA]

IDA法のオプション名	デフォルト値
"ImplicitSolver"	"Newton"	陰方程式の解法
"MaxDifferenceOrder"	5	BDFの最大次数

IDAのメソッドオプション

NDSolveから返されたInterpolatingFunctionオブジェクトで計算された中間値の厳密な確度が重要なときは，NDSolveメソッドオプション設定のInterpolationOrder->Allを使うとよい．この設定では，時間ステップ間の解を表すために，「密な出力」と呼ばれることもあるメソッドの次数に基づいた補間が使用される．デフォルトで，NDSolveは視覚的に満足できる解を表すために最小限のデータを保存する．データ量を最小限に保つことは，より複雑な解の場合に必要なメモリと時間を節約することになる．

最小限の出力とメソッドの完全な補間次数の差が顕著に現れる例として，厳密解がDSolveで計算できる簡単な線形方程式の解から導いた数量を追跡してみよう．

これは，解x[t]およびy[t]に基づいて，数量を時間の関数として与える関数fを定義する：

Wolfram Language code: f[t_] := x[t]^2 + y[t]^2;

線形方程式を初期条件とともに定義する：

Wolfram Language code:

eqns = {x'[t] == x[t] - 2 y[t], y'[t] == x[t] + y[t]};
ics = {x[0] == 1, y[0] == 1};

時間の関数としてのfの厳密値は，DSolveを使って記号的に計算することができる：

Wolfram Language code: Subscript[f, exact][t_] = First[f[t] /. DSolve[{eqns, ics}, {x, y}, t]]

厳密解は，密な出力で，あるいは密な出力なしで計算された解と比較される．

数量を追跡するには，微分方程式の数値解からその数量を導出する関数を生成するのが最も簡単である：

Wolfram Language code: f1[t_] = First[f[t] /. NDSolve[{eqns, ics}, {x, y}, {t, 0, 1}]]

これは密な出力で計算することもできる：

Wolfram Language code: Subscript[f1, dense][t_] = First[f[t] /. NDSolve[{eqns, ics}, {x, y}, {t, 0, 1}, InterpolationOrder -> All]]

以下は計算された2つの解の相対誤差を示すプロットである．時間ステップで計算された解は，黒い点で示してある．デフォルトの出力誤差は灰色で，密な出力誤差は黒で示してある：

Wolfram Language code: Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Wolfram Language code:

t1  = Cases[f1[t], (if_InterpolatingFunction)[t] -> 
	InterpolatingFunctionCoordinates[if], Infinity][[1, 1]];
pode  = Show[{ListPlot[Transpose[{t1, (Subscript[f, exact][t1] - f1[t1]) / Subscript[f, exact][t1]}], PlotStyle -> PointSize[.02]], Plot[(Subscript[f, exact][t] - f1[t]) / Subscript[f, exact][t], {t, 0, 1}, PlotStyle -> RGBColor[.8, .8, .8]], Plot[(Subscript[f, exact][t] - Subscript[f1, dense][t]) / Subscript[f, exact][t], {t, 0, 1}]}, PlotRange -> All]

このプロットから，時間ステップが大きくなるにつれ，デフォルトの解の出力の誤差が時間ステップ間で格段に大きくなることが明確である．蜜な出力の解は，時間ステップの間でもソルバによって計算された解を表す．しかし，密な出力の解が取るスペースはずっと大きくなることを覚えておかねばならない．

解から導きたい数量が複雑な場合，その数量を代数的量として与えることで許容誤差内に局所的に収め，ソルバに強制的に誤差制御させるようにすることができる．

従属変数を興味のある数量に等しく設定する代数方程式で，従属変数を加える：

Wolfram Language code: f2[t_] = First[g[t] /. NDSolve[{eqns, ics, g[t] == f[t]}, {x, y, g}, {t, 0, 1}]]

微分代数方程式の解の相対誤差を比較するプロットを作成する．IDAの時間刻みは青い点で示されている：

Wolfram Language code:

t2  = InterpolatingFunctionCoordinates[Head[f2[t]]][[1]];
Show[{ListPlot[Transpose[{t2, (Subscript[f, exact][t2] - f2[t2]) / Subscript[f, exact][t2]}], PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.02]}], Plot[(Subscript[f, exact][t] - f2[t]) / Subscript[f, exact][t], {t, 0, 1}, PlotStyle -> RGBColor[.7, .7, 1]]}, PlotRange -> All]

微分代数方程式の解は，関心数量の誤差を制御するためにずっと多くのステップを必要とするが，密な出力を使うよりもメモリ使用量はずっと少なくてすむ．

次はデフォルトのサイズ，密な出力のサイズ，微分代数方程式の解のサイズを比較する：

Wolfram Language code: ByteCount /@ {f1[t], Subscript[f1, dense][t], f2[t]}

これから先は，"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を参照のこと

"GMRES"メソッドオプション

例題として，二次元のバーガーズ(Burgers)方程式を挙げる．

これは通常常微分方程式ソルバで解かれるが，この例ではDAEソルバを使って2つのことが行われる．まず境界条件が代数条件として実行される．次にNDSolveが代数項を使って保存差分法を使うよう強制的される．比較のため，初期条件，境界条件に既知の厳密解を使う．次の計算は，この2つのアプローチの違いを示す．

以下でバーガーズ方程式を満足する関数を定義する：

Wolfram Language code: Bsol[t_, x_, y_] = 1 / (1 + Exp[(x + y - t) / (2 ν)]);

厳密解に矛盾しない単位正方形の初期条件と境界条件を定義する：

Wolfram Language code:

ic = u[0, x, y] == Bsol[0, x, y];
bc = {
	u[t, 0, y] == Bsol[t, 0, y], u[t, 1, y] == Bsol[t, 1, y], 
	u[t, x, 0] == Bsol[t, x, 0], u[t, x, 1] == Bsol[t, x, 1]};

微分方程式を定義する：

Wolfram Language code:

de = D[u[t, x, y], t] == ν( D[u[t, x, y], x, x]  + D[u[t, x, y], y, y])  - u[t, x, y](D[u[t, x, y], x] + D[u[t, x, y], y]);

妥当な時間で解が見付けられるような値に拡散定数

を割り当てて，

における解のプロットを表示する：

Wolfram Language code:

ν = 0.025;
Plot3D[Bsol[1, x, y], {x, 0, 1}, {y, 0, 1}]

プロットから，のとき前面が非常に急な傾斜が見られる．これは一定の速度で移動する．

NDSolveのデフォルト設定とIDA法を使って問題を解く．ただし例外として"MethodOfLines"には"DifferentiateBoundaryConditions"オプションを設定し，NDSolveが境界条件を代数的に扱うようにする：

Wolfram Language code:

Timing[sol = NDSolve[{de, ic, bc}, u, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}, Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> False}]]

比較できる厳密解があるので，全体的な解の確度も比較できる．

すべての時間刻みにおけるグリッド点で，厳密解と計算された解との最大偏差を求める関数を定義する：

Wolfram Language code:

errfun[sol_] := Module[{ifun = First[u /. sol], grid, exvals, gvals}, 
	grid = InterpolatingFunctionGrid[ifun];
	gvals = InterpolatingFunctionValuesOnGrid[ifun];
	exvals = Apply[Bsol, Transpose[grid, RotateLeft[Range[ArrayDepth[grid]], 1]]];
	Max[Abs[exvals - gvals]]]

計算された解の最大誤差を計算する：

Wolfram Language code: errfun[sol]

次では，IDA法の異なるオプション設定による差を比較する．オプション設定を強調するために，解の計算時間を測定して最大誤差を与える関数を定義する．

IDAの異なるオプション設定を比較する関数を定義する：

Wolfram Language code:

TimeSolution[idaopts___] := Module[{time, err, steps}, 
	time = First[Timing[sol = NDSolve[{de, ic, bc}, u, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}, Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> False, Method -> {IDA, idaopts}}]]];
	err = errfun[sol];
	steps = Length[First[InterpolatingFunctionCoordinates[First[u /. sol]]]] "Steps";
	{"Time" -> time, "Error"  -> err, steps}]

前のデフォルトを使うオプションはひとつもない：

Wolfram Language code: TimeSolution[]

"Band"法を使う：

Wolfram Language code: TimeSolution["ImplicitSolver" -> {"Newton", "LinearSolveMethod" -> "Band"}]

ここでは，ヤコビアンのバンド幅が比較的大きい四階導関数が使われている，境界付近で使われている片側のステンシルにより上と下に幅が加わっている等の理由から，"Band"法はあまり有効ではない．バンド幅は明示的に指定することができる．

幅に1方向のみの差分のステンシルを含むよう設定して"Band"法を使う：

Wolfram Language code:

TimeSolution["ImplicitSolver" -> {"Newton", 
	"LinearSolveMethod" -> {"Band", "BandWidth" -> 3}}]

解法の計算時間は短いが，誤差はわずかに大きく，時間刻みの合計数はかなり多くなっている．問題がさらに硬ければ，もう一方からの情報が欠けているため反復は収束しなかったであろう．理想的には，バンド幅は情報を次元のどこからも削除すべきではない．

と

方向で使われたグリッドを計算し，使われた点の数を表示する：

Wolfram Language code:

{X, Y} = InterpolatingFunctionCoordinates[First[u /. sol]][[{2, 3}]];
{nx, ny} = {Length[X], Length[Y]}

幅を両方向のステンシルの少なくとも一部を含むよう設定して"Band"法を使う：

Wolfram Language code:

TimeSolution["ImplicitSolver" -> {"Newton", 
	 "LinearSolveMethod" -> {"Band", "BandWidth" -> 51}}]

バンド幅の設定をより適切なものにしても，解の計算にかかる時間はデフォルトのバンド幅の場合よりもわずかに速いが，デフォルトの線形ソルバを使うよりも遅い．"Band"法は二次元問題で有効な場合もあるが，通常一次元問題で最も有効である．

前処理を使わずに"GMRES"の陰的ソルバを使って解を計算する：

Wolfram Language code: TimeSolution["ImplicitSolver" -> "GMRES"]

前処理をせずに"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に含まれるような種類の因数分解を使うのはよいことである．

対角行列の場合は，単に相互行列を使うだけで反転が可能である．対角前処理を設定する上で最も難しい部分は，境界の値がそれによって影響を受けてはならないという点である．

前処理を計算するための微分行列の対角要素を見付ける：

Wolfram Language code:

DM = NDSolve`FiniteDifferenceDerivative[{2, 0}, {X, Y}]@"DifferentiationMatrix" + NDSolve`FiniteDifferenceDerivative[{0, 2}, {X, Y}]@"DifferentiationMatrix";
Short[diag = Tr[DM, List]]

これにより，境界で簡単な代数状態を満足する要素が，平坦化された解のべクトルとなっている位置を得る：

Wolfram Language code:

bound = SparseArray[{{i_, 1} -> 1., {i_, ny} -> 1., {1, i_} -> 1., {nx, i_} -> 1.}, {nx, ny}, 0.];
Short[pos = Drop[ArrayRules[Flatten[bound]], -1][[All, 1, 1]]]

前処理を有効に反転するために呼び出される関数を設定する関数を定義する．対角行列の場合は，相互行列を取るだけで逆になる：

Wolfram Language code:

pfun[t_, x_, xp_, c_] := Module[{d, dd}, 
	d = 1. / (c - ν diag);
	d[[pos]] = 1.;
	Function[# dd] /. dd -> d]

前処理された"GMRES"法を使って解を求める：

Wolfram Language code: TimeSolution["ImplicitSolver" -> {"GMRES", "Preconditioner" -> pfun}]

未完成の前処理を使っても，"GMRES"法は疎行列直接解法を使うよりも速く解を計算する．

高階の一次微分を伴う偏微分方程式の離散化や偏微分方程式系では，対応するNDSolve`StateDataオブジェクトを見て，前処理の構造形式が正確に得られるように，変数の順序を決めなければならない．

Top

その他のラーニングリソース

テクニカルサポート

Wolframソリューション

教育のためのWolframソリューション

使い始める

Grow Your Skills

Wolframと繋がる

大人用の教育プログラム

若者のための教育プログラム

読む

NDSolveのIDAメソッド