境界値問題の数値解法

狙い撃ち法

狙い撃ち("Shooting")法は,境界条件をある点における多変量関数として考え,根を与える初期条件を見付けということに境界値問題を還元することにより作用する.狙い撃ち法の利点は,初期値問題のためのメソッドのスピードと適応性を利用するという点である.しかし,有限差分法や選点法ほどロバストなメソッドではないという欠点もある.成長モードの初期値問題の中には,たとえ境界値問題自体が非常に適切で安定していても,本質的に不安定なものもある.

次の境界値問題

を考える.狙い撃ち法では となるような初期条件 が探索される.初期条件を変化させているので, は初期条件の関数として考えるのが妥当であり,狙い撃ちは下が成り立つような を見付けるものと考えることができる.

の関数を設定した後で,問題は実際にFindRootに渡され,根を与える初期値 が探索される.デフォルトメソッドではヤコビアンの計算を含むニュートン法が使われる.ヤコビアンは有限差分を使って計算することができるが,初期値問題の解はその初期条件に敏感で適当に正確な導関数値を得るには影響が大きすぎる可能性があるので,ヤコビアンは常微分方程式の解として計算した方がよい.

線形化とニュートン法

線形問題は以下の式

で記述することができる.ここで は行列,はベクトルでどちらも に依存する可能性がある.は定数ベクトル, は定数行列である.

とする.ここで, について初期値条件および境界条件両方を微分すると下が得られる.

の関数と考えると線形なので,が得られる.従って となるような の値は,どのような特定の初期条件 の場合についても

を満足する.

非線形問題の場合,を非線形常微分方程式系のヤコビアンとし, 番目の境界条件のヤコビアンとする.線形化された系に対する の計算により,特定の初期条件に対する非線形系のヤコビアンが得られ,ニュートン反復へとつながる.

"StartingInitialConditions"

境界値問題では,初期値問題の場合のように一意である保証はない.は1つの解だけを見付ける.FindRootが非線形代数方程式系に対して見付ける特定の解に,初期値を変更することで影響が与えられるように,で見付かった解を反復を開始する異なる初期条件を与えることで変更することができる.

メソッドのオプションであり,狙い撃ちプロセスを開始する初期条件の値と位置を指定することができる.

狙い撃ち法はデフォルトでは初期条件なしで開始されるので,解がない場合はそれが返される.

次では のときの境界値問題 の非常に単純な解を計算する.
In[1]:=
Click for copyable input
Out[2]=

デフォルトで,は区間の左側から始まり,時間を追うごとに前進する.後退した方がよい場合,あるいは区間の中間付近の点から行った方がよい場合もある.

線形境界値問題を考える.

これは以下の解を持つ.

が適切な値の場合,から始まる初期値問題は および 項が成長するために不安定になる.同様に から始めると,前進方向の項ほどは大きくないが, 項から不安定性が生じる. のある値を超えると,どちらかの方向の感度が大きくなりすぎるために,狙い撃ち法ではよい解が得られなくなる.しかし,その点までなら,両方向での成長がバランスの取れている区間内の点を選ぶことにより,最適な解が得られる.

次は方程式,境界条件,厳密解をWolfram言語入力として返す.
In[3]:=
Click for copyable input
次はデフォルトの から狙い撃ちを開始して の系を解く.
In[6]:=
Click for copyable input
Out[6]=

からを行うと,不正な行列について,さらに境界条件が想定されるように満足されないとの警告が出力される.これは における小さい誤差が により増幅されるためである.これの逆数は,局所的切り捨て誤差の大きさの次数と同じなので,このプロットに見られるような可視の誤差は驚くほどのものではない.逆方向では,大きさはずっと小さい であるため,解はずっとよいものである.

からの狙い撃ち法を使って解を計算する.
In[7]:=
Click for copyable input
Out[7]=

選ぶのに適した点は,各方向での感度がバランスの取れている点であり,ここではおおよそ である.この点では,の誤差はまだ妥当に制御下におかれる.

からの狙い撃ち法を使って の解を計算する.
In[8]:=
Click for copyable input
Out[8]=

オプションの要約

オプション名
デフォルト値
"StartingInitialConditions"Automatic狙い撃ち法を開始する初期座標
"ImplicitSolver"Automatic境界条件により定義された陰的方程式を解くために使われるメソッド.これはFindRootMethodオプションで取ることのできる値でなければならない
"MaxIterations"Automatic陰的ソルバメソッドで使用する反復の回数
"Method"Automatic常微分方程式系の積分に使用するメソッド

法のオプション

追跡法

追跡("Chasing")法という方法はゲルファント(Gelfand)とロクチェフスキー(Lokutsiyevskii)の原稿に記述されていたもので,初めて英語で[BZ65]に発表され,[Na79]でさらに詳しく説明された.追跡法では,まず1組の補助問題を立て,それを解いて境界のひとつにおける初期条件を求める.一旦初期条件が決まったら,初期値問題を解く通常の方法が適用できる.追跡法は,実質的に,問題の線形性を十分に活用する狙い撃ち法なのである.

以下の線形常微分方程式について考えてみよう.

この方程式で,であり,は係数行列,は以下のような 個の線形境界条件を持つ非同次係数ベクトルである.

ここで は係数ベクトルである.これから,拡大された同一系を構築する.

ここで以下が成り立つ.

追跡法とは,が成り立つようなベクトル関数を求めることである.一旦関数が分かり,十分な境界条件があれば,以下の式

を解くことで,通常の初期値問題のソルバに用いることができる初期条件を定めることができる.の最初の要素が常に1なので,系(3)の解は零ではない点に注意のこと.

このように,境界条件問題を解くことは,についての補助問題を解くことに還元される.についての式を微分すると,

が得られる.これに微分方程式を代入し

転置する.

これはすべての解 に当てはまるはずなので, の初期値問題が得られる.

すべての境界値問題の解が求めたいところ を とすると,Wolfram言語はNDSolveを使うだけで の補助問題を まで積分することによって解く.結果は(3)の行列に組み込まれる.これを について解くと,初期値問題が得られる.この初期値問題をNDSolveで積分すると,返された解が与えられる.

この方法のこの変形については MathSource パッケージ[R98]でさらに説明が加えられ,使われている.このパッケージを使うと線形固有値問題を解くこともできる.

補助問題を設定するには,これ以外にもGel'fandとLokutsiyevskiiが最初に提唱した方法により近い非線形の方法がある.境界条件が線形に独立であると仮定する(そうでなければ,問題が十分に特定されない).すると,各 に少なくとも1つの非零の構成要素があることになる.一般性を失うことなく,と仮定する.次に, の他の構成要素 の観点から について解く.ここでであり である.(5) を使い, を置き換え,の他の構成要素の観点から を考えると,以下の非線形方程式が得られる.

上の式で から 番目の列を除いたものであり, 番目の列である.非線形の方法は線形の方法よりも数値的に安定していることがあるが,実数直線上の積分が特異値に至る可能性があるという不都合もある.この問題は複素平面の等高線上で積分することで除外することができるが,複素平面上の積分は実数直線上の単純な積分よりも数値誤差が大きくなることが多いので,実際には非線形の方法の方が線形の方法よりも常によい結果を与えるとはいえない.この理由と,一般的な速さという理由で,Wolfram言語ではデフォルトで線形の方法を使っている.

これで二階常微分方程式の二点境界値問題を解く.
In[9]:=
Click for copyable input
Out[9]=
解のプロットである.
In[10]:=
Click for copyable input
Out[10]=
このソルバは線形方程式系の複数点の境界条件問題を解くことができる.各境界方程式は,一意的な値 になければならない点に注意のこと.
In[11]:=
Click for copyable input
Out[12]=

通常,境界値方程式がEqualの公差限度を満足することは期待できない.

以下で境界条件が「満足されている」かどうか検証する.
In[13]:=
Click for copyable input
Out[13]=

通常境界条件は,NDSolveのオプションであるAccuracyGoalPrecisionGoalの最大許容度のときしか満足されない.一般に,確度と精度は大域的ではなく局所的な誤差制御に使われるので,実際の確度と精度はこれよりも劣る.

以下は,境界条件のそれぞれにおいて残差誤差を検証する.
In[14]:=
Click for copyable input
Out[14]=

NDSolveに解のない問題を与えると,数値誤差によってそれが解き得る問題に見えることがある.通常NDSolveは警告メッセージを発する.

これは解をもたない境界値問題である.
In[15]:=
Click for copyable input
Out[15]=

この場合,解がないために全区間に渡る積分ができない.

方程式が悪条件となる別の例として,境界条件が一意解を与えない場合がある.

次は一意解を持たない境界値問題である.この一般解がDSolveで記号的に計算されて表示される.
In[16]:=
Click for copyable input
Out[16]=
初期条件について解くべき行列が特異であるにも関わらず解を持つため,NDSolveは警告メッセージを発する.
In[17]:=
Click for copyable input
Out[17]=
補間点にフィットすると,どの解が見付かったのかが分かる.これにより,誤差のプロットが最もフィットする実際の解と相対的になる.
In[18]:=
Click for copyable input
Out[22]=

通常,Wolfram言語が使うデフォルト値でうまくいくが,NDSolveにオプションMethod->{"Chasing",chasing options}を与えて追跡法をコントロールすることもできる.次は が取り得る値の表である.

オプション名
デフォルト値
MethodAutomatic追跡法が生成した初期値問題の計算に使う数値メソッド
"ExtraPrecision"0補助初期値問題を解くために使う追加精度の桁数
"ChasingType""LinearChasing"使用する追跡のタイプ.または

NDSolveメソッドのオプション

メソッドには2つのオプションがある.

オプション名
デフォルト値
"ContourType""Ellipse"複素平面における積分が必要なときに使用する等高線の形.あるいは
"ContourRatio"1/10等高線の長さに対する幅の割合.一般にこの数が小さいほど結果が正確になるが,そうすると方程式を解くのが数値的により困難になる.

メソッドオプションのオプション

これらのオプション,特には,計算された初期条件に対する感度が強い場合に極めて有効である.

以下はDSolveを使って記号的に計算した単純な解を持つ境界値問題である.
In[23]:=
Click for copyable input
Out[24]=
これはNDSolveで追跡法を使って計算した解誤差を示している.
In[25]:=
Click for copyable input
Out[26]=
初期条件を解くために精度を上げると誤差が著しく小さくなる.
In[27]:=
Click for copyable input
Out[28]=

一旦初期条件が見付かると,誤差の大部分は解の計算によるものなので,これ以上精度を上げてもあまり役に立たない.誤差を小さくするためには,NDSolveAccuracyGoalオプションとPrecisionGoalオプションの値をより厳しくする必要がある.

以下では初期条件の計算の精度を上げ,AccuracyGoalオプションとPrecisionGoalオプションの設定もより厳しくしてある.
In[29]:=
Click for copyable input
Out[30]=

パラメータを持つ境界値問題

境界値問題が生じるような分野の多くでは,問題自体に,所望の解の一部となり得る固有値等の未定のパラメータがある場合がある.このようなパラメータを従属変数として導入することにより,問題が標準形式の境界値問題として書ける場合がよくある.

例えば,チャネルの流れは次のようにモデル化することができる.

ここで (レイノルズ数)は与えられるが は未定である.

の値を見付けるためには,方程式 を加えるだけである.

以下で について を持つ流れ問題を解き,解 をプロットし, の値を返す.
In[31]:=
Click for copyable input
Out[31]=