コンポーネントとデータ構造
はじめに
NDSolveはいくつかの基本的なステップに分けることができる.高度な使い方では,これらのステップを別々に実行するためのコンポーネントにアクセスした方が便利なこともある.
NDSolveはこれらのステップをすべて内部的に実行するため,通常ユーザはその詳細を見ることはない.
NDSolve`ProcessEquationsは,微分方程式系を初期値問題,境界値問題,微分代数問題,偏微分問題等に分ける.また適切なデフォルトの積分メソッドを選択し,メインのNDSolve`StateDataデータ構造を構築する.
NDSolve`Iterateは数値解法を進める.最初の起動(複数の場合もある)で数値積分法が初期化される.
NDSolve`ProcessSolutionsは数値データをInterpolatingFunctionに変換し,それぞれの解を表すようにする.
NDSolve`ProcessEquationsは微分方程式系を解くのに,全体の時間の大部分を使用することがある.そのような場合は,このステップを1度だけ実行し,NDSolve`Reinitializeを使って異なるオプションや初期条件で繰返しとくようにするとよい.
例題
ndssdata =
First[NDSolve`ProcessEquations[{y''[t] + y[t] == 0, y[0] == 1, y'[0] == 0}, {y, y'}, t, Method -> "ExplicitRungeKutta"]]NDSolve`Iterate[ndssdata, 10]ndsol = NDSolve`ProcessSolutions[ndssdata]ParametricPlot[{y[t], y'[t]} /. ndsol, {t, 0, 10}]ProcessEquations
NDSolveを使うあらゆる解法の第一段階は,指定された方程式を,実際の積分アルゴリズムが効率よくアクセスできる形式に処理することである.この段階には最低限でも,各変数の微分階数の決定,一階の系を得るために必要な代替の生成,関数についての関数の時間微分の解決,その結果の"NumericalFunction"オブジェクトへの変換が含まれる.同じ方程式の集合にこの処理を繰り返す時間を節約したい場合,または,数値積分処理をもっと自分でコントロールしたい場合は,NDSolve`ProcessEquationsで処理段階を別々に実行することができる.
| NDSolve`ProcessEquations[{eqn1,eqn2,…},{u1,u2,…},t] | |
| 関数{u1,u2,…}の微分方程式{eqn1,eqn2,…}を標準形に処理する.関数についての関数の時間微分の解と,各解に関連するデータを含むNDSolve`StateDataオブジェクトのリストを返す.t は,NDSolveのように値の範囲のリストで指定することができる | |
| NDSolve`ProcessEquations[{eqn1,eqn2,…},{u1,u2,…},{x1,x1min,x1max},{x2,x2min,x2max},…] | |
| 関数{u1,u2,…}の偏微分方程式{eqn1,eqn2,…}を標準形に処理する.関数について,その関数の時間微分の解と,各解に関連するデータを含むNDSolve`StateDataオブジェクトのリストを返す.xj が時間変数ならば,これは境界 xj min,xj max で指定する必要はない | |
NDSolve方程式の処理
NDSolve`ProcessEquations[{y'[x] ^ 2 == y[x] + x, y[0] == 1}, y, x]Reinitialize
より高度な問題の解を求める際,同じ微分方程式を異なる初期条件で繰返し解くことは珍しくない.方程式を処理するのに,微分方程式を数値的に積分するのと同じくらいの時間がかかることもある.このような場合,簡単に新しい初期値を与えることができると非常に便利である.
| NDSolve`Reinitialize[state,conditions] | 方程式と変数を,NDSolve`StateDataオブジェクト state を作成するのに使われたものと同じと仮定して,方程式 conditions の関数の初期値に対する可能な各解について1つの新しいNDSolve`StateDataオブジェクトを作成し,そのリストを形成する |
state = First[NDSolve`ProcessEquations[{x''[t] + x[t] == 0, x[0] == 0, x'[0] == 1}, x, t]]newstate = NDSolve`Reinitialize[state, {x[1] ^ 3 == 1, x'[1] == 0}]NDSolve`Reinitializeを使うと,同じ微分方程式を多くの異なる初期条件で解く必要があるときに,境界値問題の狙い撃ち法のように,計算時間を短縮することができる.
NDSolveオプションの部分集合は,NDSolve`Reinitializeのオプションとして指定することができる.
newstate = NDSolve`Reinitialize[state, {x[0] == 0, x'[0] == 1}, StartingStepSize -> 1 / 10]解の反復
NDSolve`StateDataオブジェクトの使い方で大切なもののひとつに,積分をよりコントロールできるというものがある.問題によっては,解をチェックしてやり直すことや,状況に応じてパラメータを変更することが適切なものもある.
state = First[NDSolve`ProcessEquations[{x''[t] + (1 + 4 UnitStep[Sin[t]])x[t] == 0, x[0] == 1, x'[0] == 0}, x, t, Method -> "ExplicitRungeKutta"]]NDSolve`ProcessEquationsを使うときは,変数
の範囲を明示的に与える必要はない.その情報は,方程式を解くための形式に変更するのに必要ではないからである.しかし,PDEでは空間変数の範囲は適切な打切りを決定するのに必須であるため,その情報を与えなければならない.
NDSolve`Iterate[state, 1]NDSolve`Iterateは,変数状態に割り当てられたNDSolve`StateDataオブジェクトを変更するので,値を返さない.従って,このコマンドは「指標によるリストの操作」で説明してあるように,リストの部分を設定するのと類似の方法で変数の値に影響を及ぼす.state の値は積分の結果現在時間を表示するようになったので,その値が変化したということが分かる.
stateさらに積分したい場合は,時間の値を大きくして再度NDSolve`Iterateを呼び出してもよい.
NDSolve`Iterate[state, 3]最初の現在時間よりも早い時間を指定することができる.この場合,積分は時間について逆に進行する.
NDSolve`Iterate[state, -Pi / 2]NDSolve`Iterateを使うと,中止する中間時間を指定することができる.これは,不連続を避けたい場合等に便利である.通常この手法は,この例で使用されている陽的ルンゲ・クッタ(Runge-Kutta)法のような,いわゆる1ステップ法で使った方が効果的である.しかし,これは通常デフォルトメソッドのNDSolveでも動作する.
NDSolve`Iterate[state, π Range[10]]解の関数の取得
ある時間まで系を積分したら,現在の解の値が調べられ,それまでに計算された解を表わす近似関数が生成できることが通常求められる.NDSolve`ProcessSolutionsコマンドを使うと,両方行うことができる.
| NDSolve`ProcessSolutions[state] | state で計算された解をInterpolatingFunctionオブジェクトの規則のリストとして得る |
InterpolatingFunctionオブジェクトとして解を得る
sol = NDSolve`ProcessSolutions[state]Plot[Evaluate[x[t] /. sol], {t, 0, 10 Pi}]NDSolveを直接使った場合と同様に,NDSolve`ProcessEquationsの第2引数で指定した各関数に規則がある.解の指定された要素だけが,InterpolatingFunctionオブジェクトが生成できる様式で保存される.
sol = NDSolve`ProcessSolutions[state, "Forward"]方向 dir に指定できるのは,"Forward"か"Backward"である.これらはそれぞれ初期条件から正方向の積分と逆方向の積分を意味する.
| "Forward" | 一時変数の値を増加させる方向での積分 |
| "Backward" | 一時変数の値を減少させる方向での積分 |
| "Active" | 現在積分が行われている方向での積分.通常この値はアクティブな積分中に使われるメソッドの初期化からのみ呼び出される |
NDSolve`ProcessSolutionの出力は,常に独立変数の指定値における従属変数,あるいは,保存された値全体に補間された従属変数で与えられる.つまり,偏微分方程式が積分されているとき,空間変数上の従属変数を表す結果が得られるということである.
state = First[NDSolve`ProcessEquations[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == Cos[π / 2 x], u[t, 0] == 1, u[t, 1] == 0}, u, t, {x, 0, 1}]];
NDSolve`Iterate[state, {-1 / 8, 2}]NDSolve`ProcessSolutions[state, "Forward"]解は空間変数
を補間するInterpolatingFunctionオブジェクトとして与えられる.
NDSolve`ProcessSolutions[state, "Backward"]
偏微分方程式の現在の解を処理するとき,空間誤差推定がチェックされる(これは非常に時間がかかることがあるので,通常,解の生成時以外には行われない).空間誤差推定が過剰であるため,NDSolve::eerrというメッセージが表示される.熱方程式と「backward」という言葉との関連付けで不安定を意味するものにより,この例題で何が誤っているかのヒントが得られる.
Plot[Evaluate[u[-0.125, x] /. %], {x, 0, 1}]解のプロットにより,実際に不安定性が問題であることが分かる.
熱方程式の例題は簡単なので,逆方向時間の解に問題があるということが分かるが,NDSolve`IterateとNDSolve`ProcessSolutionsを使って偏微分方程式の解を監視すると,結果的に希望通りの正確さにはならない解を計算する手間を省くことができる.別の簡単な監視形式は以下の通りである.
L = -10;
state = First[NDSolve`ProcessEquations[{D[u[t, x, y], t, t] == D[u[t, x, y], x, x] + D[u[t, x, y], y, y] - Sin[u[t, x, y]],
u[0, x, y] == Exp[-(x^2 + y^2)], Derivative[1, 0, 0][u][0, x, y] == 0,
u[t, -L, y] == u[t, L, y], u[t, x, -L] == u[t, x, L]}, u, t, {x, -L, L}, {y, -L, L}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"}}]];
GraphicsGrid[Partition[Table[
NDSolve`Iterate[state, τ];
Plot3D[Evaluate[u[τ, x, y] /. NDSolve`ProcessSolutions[state, "Forward"]], {x, -L, L}, {y, -L, L}, PlotRange -> {-1 / 4, 1 / 4}],
{τ, 0., 20., 5.}], 2]]このように解を監視すると,見付かった解で十分であると思ったときに計算を中断することが通常可能である.それでもまだNDSolve`StateDataオブジェクトを使うと,計算された解を得ることができる.
NDSolve`StateDataのプロパティ
NDSolve`StateDataオブジェクトには多くの情報が含まれているが,どこに情報が保存されているかを分かりやすくするような配列ではなく,解を反復するのが簡単になるように配列されている.しかし,状態データオブジェクトから情報が得たい場合もあるだろう.この目的のため,情報へのアクセスが簡単にできるようなプロパティがいくつか定義されている.
NDSolve`StateDataオブジェクトからの情報の最も重要な用途の一つに,積分メソッドの設定というものがある.この例は「NDSolveのメソッドプラグインフレームワーク」に示してある.
| state@"VariableData" | 解データに従って構造化されたすべての変数のリストのリストを与える |
| state@"VariableDimensions" | 変数それぞれの次元のリストを与える |
| state@"VariablePositions" | 解データに含まれる,各変数の位置のリストを与える |
| state@"NumericalFunction" | 一時変数 t について,解ベクトルの導関数を評価するのに使われるNumericalFunctionオブジェクトを与える |
| state@"ProcessExpression"[args,expr,dims] | |
| 数値的に expr を評価するためのNumericalFunctionオブジェクトを与える state を生成するのにNDSolve`ProcessEquationsにより使われた変数変換と同じものを使って,式 expr を処理する.args は数値関数の引数であり,All,あるいは系の従属変数である引数のリストのどちらかである.dims はAutomatic,あるいは数値的関数の結果に想定される次元を与える明示的なリストである | |
| state@"SystemSize" | 解かれようとしている一階常微分方程式の有効な数を与える |
| state@"MaxSteps" | 微分方程式を反復するのに許される最大回数 |
| state@"WorkingPrecision" | 方程式を解くために使われる作業精度を与える |
| state@"Norm" | 誤差を測定するために使われる,スケールされたノルム |
NDSolve`StateDataオブジェクト state のプロパティ
state = First[NDSolve`ProcessEquations[{x''[t] + Sin[x[t]] == 0, x[0] == 3, x'[0] == 0}, x, t]]nf = state@"SystemSize"2つのスカラー常微分方程式があるので,系の大きさは2になる.
利用できる情報の多くは,現在の解の値に依存する. NDSolveは前進と後退(初期条件時から)の積分方向で異なるデータを保存する.このデータには次の方向性によってアクセスすることができる.
| state@"SolutionData"[dir] | 解データの現在の値を積分方向 dir で与える |
| state@"TimeStep"[dir] | 次のステップの時間刻みを積分方向 dir で与える |
| state@"TimeStepsUsed"[dir] | 現在の時間を得るために使う時間刻みの数を積分方向 dir で与える |
| state@"MethodData"[dir] | 使用されるメソッドデータオブジェクトを積分方向 dir で与える |
NDSolve`StateDataオブジェクト state の方向性
方向引数が省略されると,両方向のデータを含むリスト(初期条件において1つの要素を持つリスト)が返される.省略しない場合は,方向は前のセクションで指定されているように"Forward","Backward","Active"のいずれかである.
state@"TimeStep"NDSolve`Iterate[state, 10];
state@"TimeStep"["Forward"]state@"TimeStepsUsed"state@"SolutionData"SolutionData
解データは,それぞれの型の解の成分の値のリストのリストで構成されている.解データの成分は,時間,従属変数,離散変数等の異なる型の変数に対応する.特定の問題に使われていない成分はすべて空のリスト({})かNoneとして与えられる.
| NDSolve`SolutionDataComponent[sd,c] | 解データ sd から解の成分 c を得る |
| NDSolve`SetSolutionDataComponent[sd,c,v] | シンボル sd に割り当てられた解データにおいて解の成分 c を v に設定する |
| nf@"SolutionDataComponents" | NDSolveで使われるNumericalFunction nf の引数として使われる解データの成分を与える |
簡略名
| 完全名 |
成分
|
| "T" | "Time" | 現在の時間 |
| "S" | "Space" | 空間離散化 |
| "X" | "DependentVariables" | 離散変数以外の従属変数の現行値 |
| "X'" | "TemporalDerivatives" | 従属変数の時間導関数 |
| "D" | "DiscreteVariables" | 値の連続範囲を取る離散変数 |
| "ID" | "IndexedDiscreteVariables" | 値の有限集合を取る離散変数 |
| "P" | "Parameters" | 感度が計算されるパラメータ以外のパラメータ |
| "SP" | "SensitivityParameters" | 感度が計算されるパラメータ |
state = First[NDSolve`ProcessEquations[{x''[t] + a[t] Sin[x[t]] == 0, x[0] == 3, x'[0] == 0, a[0] == 1, WhenEvent[x[t] == 0, a[t] -> a[t] + 1]}, {a, x}, t, DiscreteVariables -> a]];
NDSolve`Iterate[state, 10];
statevars = state@"Variables"vars = state@"VariableData"変数データと解データの両方からのコンポーネントを抽出するためには,NDSolve`SolutionDataComponentを使うことができる.
NDSolve`SolutionDataComponent[vars, "DependentVariables"]{sdb, sdf} = state@"SolutionData"NDSolve`SolutionDataComponent[sdf, "DependentVariables"]NDSolve`SolutionDataComponent[sdb, "DiscreteVariables"]NDSolve`SolutionDataComponent[sdb, "IndexedDiscreteVariables"]解データは,実質的にNDSolveが解を計算するために使う生データである.上の例題では,二階方程式は自動的にで一階の系に簡約されるので,従属変数の成分は長さ2のリストである.現行の解の部分をはっきり知りたいときは,NDSolve`ProcessSolutionsを使うことができる.
NDSolve`ProcessSolutions[state, "Forward"]NumericalFunction
時間積分を行うときにNDSolveが使用するメソッドは,
となるような常微分方程式の右辺の関数,あるいは
となるような微分代数方程式の残差関数で使うことができる.ここで
は実数であり,
,
,その他使用されている解データの要素はベクトルである.ベクトル
および
によって表される実際の変数は,もっと複雑な構造をしている可能性があるので,NDSolveは,ベクトル引数を取り,根底の関数に基づいてベクトルを返す中間段階のNumericalFunctionオブジェクトを使う.
state = First[NDSolve`ProcessEquations[{y'[t] == {{.1, .2}, {.3, .4}}.y[t] - 10 Sign[z[t]], z'[t] == y[t].y[t], y[0] == {2, 1}, z[0] == -1}, {y, z}, t]]nf = state@"NumericalFunction"評価のためにNDSolveが構築するNumericalFunctionは,問題が実際に使われる解データの成分を引数としてのみ使用する."SolutionDataComponents"プロパティを使うと,実際にどの成分が使われるかが分かる.
comp = nf@"SolutionDataComponents"sd = state@"SolutionData"[0]args = NDSolve`SolutionDataComponent[sd, comp]指標付きの離散変数は,Sign関数の非連続性自動処理によって,"DiscontinuitySignature"変数として導入された(以下に示す解を参照のこと).
Apply[nf, args]% == NDSolve`SolutionDataComponent[sd, "X'"]与えられた問題の解データで評価することが簡単になるように,解データで評価する2つの便利な関数が含まれた.
| NDSolve`EvaluateWithSolutionData[nf,sd] | 解データのリスト sd の引数でNumericalFunction nf を評価する |
| NDSolve`EvaluateJacobianWithSolutionData[nf,sd,c] | 解データ成分 c のすべての変数について,ヤコビアン導関数を評価する |
NDSolve`EvaluateWithSolutionData[nf, sd]NumericalFunctionが実行することの一つに,導関数を自動的に処理するというものがある.デフォルトでは,記号的導関数が見付かったらそれが使われる.見付からなかったら,有限差分が使われる.
NDSolve`EvaluateJacobianWithSolutionData[nf, sd, "X"]次の例では,ベクトルの引数とDot積が使われているため有限差分が使われる.
NDSolve`Iterate[state, 1];
Plot[Evaluate[{z[t], z'[t]} /. NDSolve`ProcessSolutions[state]], {t, 0, .2}, PlotRange -> All]NDSolveの変数のある関数を評価するために,NumericalFunctionを設定しておくと便利なことがある.これは"ProcessExpression"を使って設定できる.
組込みメソッドの多くは,この過程をメソッドの初期化中に行うので,関数は何回でも効率よく評価することができる.例えば,"projection"メソッドは,その初期値から不変数の偏差を評価する関数を設定し,各ステップの後で投射を実装するために導関数を使用する.
state = First[NDSolve`ProcessEquations[{x''[t] + a[t] ^ 2 x[t] == 0, WhenEvent[x[t] > 0, a[t] -> a[t] + 1], x[0] == 1, x'[0] == 0, a[0] == 1}, x, t, DiscreteVariables -> a, MaxSteps -> Infinity]];efun = state@"ProcessExpression"[Automatic, (1/2)(x'[t]^2 + (a[t] x[t])^2), {}]einit = NDSolve`EvaluateWithSolutionData[efun, state@"SolutionData"["Forward"]]ListPlot[edata = Table[NDSolve`Iterate[state, τ];{τ, NDSolve`EvaluateWithSolutionData[efun, state@"SolutionData"["Forward"]] - einit}, {τ, 0, 50}]]Plot[Evaluate[x[t] /. NDSolve`ProcessSolutions[state]], {t, 49.99, 50}]