常微分方程式の数値解法

常微分方程式の数値解法 (じょうびぶんほうていしきのすうちかいほう、英: Numerical methods for ODEs) は、数値解析において常微分方程式を数値的に解く技術の総称である^[1][2]。

数値解法の必要性

これまで様々な自然現象 (物理現象など) を記述するために多くの常微分方程式が作られ、多くの数学者たちがその解法を探求してきたが、フックス型微分方程式^[3][4]などを除いて、式変形による計算だけで厳密に解ける常微分方程式は多くない。そのため多くの研究者たちが常微分方程式を数値的に解く技術について研究をしてきた^[1][2]。最も標準的な手法はルンゲ・クッタ法であり^{[1][2][5][A 1]}、MATLABにはode45として搭載されている。しかしこれは万能なソルバーとは言えない。例えばパンルヴェ方程式^[6][7][8]やリッカチ方程式^[9]などは非線形性によって精度の良い計算ができず、数値実験結果だけを見ていると間違った結論 (幻影解) にたどり着く危険がある。そのため^[要出典]、

• 線型多段法^[1][2]
• リープ・フロッグ法
• オイラー法の進化系
  • 後退オイラー法 (en)
  • en:semi-implicit Euler method
  • en:Euler–Maruyama method (確率微分方程式に特化した解法)
• en:exponential integrator (行列指数関数を使う解法)^{[A 2][A 3]}
• 狙い撃ち法
• en:Bulirsch–Stoer algorithm^{[A 4][A 5]}
• シンプレクティック数値積分法^{[10][11][A 6][A 7][A 8]}
• テイラー級数を用いる方法^{[A 9][A 10][A 11][A 12][A 13]}

などの新しい解法に関する研究が進められている。

初期値問題

${\displaystyle N}$ 個の変数 ${\displaystyle y=(y_{1},\cdots ,y_{N})\in \mathbb {R} ^{N}}$ に関する1階常微分方程式、すなわち ${\displaystyle [t_{0},b]\times \mathbb {R} ^{N}}$ （またはその開集合）で定義されたベクトル値連続関数 ${\displaystyle f(t,y)=(f_{1}(t,y),\cdots ,f_{N}(t,y))}$ により定まる次の方程式

${\displaystyle {\frac {dy}{dt}}=f(t,y)}$

について考える。その初期値問題 (initial-value problem) とは、初期条件 ${\displaystyle y(t_{0})=y_{0}}$ を満たす関数 ${\displaystyle y(t)}$ を求めることである^[12]。関数 ${\displaystyle f}$ が第2引数についてリプシッツ連続であるとき、すなわち定数 ${\displaystyle L}$ が存在し

${\displaystyle \left\|f(t,y_{1})-f(t,y_{2})\right\|\leq L\|y_{1}-y_{2}\|}$

（${\displaystyle \|\cdot \|}$ は ${\displaystyle \mathbb {R} ^{N}}$ のノルム）を満たすとき、その初期値問題には解が一意に存在する^[13]（ピカール・リンデレフの定理）。本節では常微分方程式の初期値問題を数値的に解くことについて議論する。

初期値問題の例

例として、電気回路の研究から導かれたファン・デル・ポール振動子について考える。その運動方程式は

${\displaystyle {\frac {d^{2}x}{dt^{2}}}-\mu (1-x^{2}){\frac {dx}{dt}}+x=0}$

であり、これは2階の常微分方程式であるものの、${\displaystyle v=dx/dt}$ とおくと

${\displaystyle {\frac {dx}{dt}}=v,\ \ {\frac {dv}{dt}}=-x+\mu (1-x^{2})v}$

という上述の形に帰着できる^[14]。

一段法

区間 ${\displaystyle [t_{0},t]}$ の厳密解を ${\displaystyle y(t)}$ とする。一段法^[15]として知られるクラスの数値解法では、離散化した時刻 ${\displaystyle t_{i}=t_{0}+ih}$ (${\displaystyle h=(t-t_{0})/n}$) での厳密解 ${\displaystyle y(t_{i})}$ の近似値 ${\displaystyle \eta _{i}}$ を

${\displaystyle \eta _{i+1}=\eta _{i}+h\Phi (t_{i},\eta _{i};h)}$

という漸化式によって定める^[16]。関数 ${\displaystyle \Phi }$ の選択が数値積分スキームを選択することに対応する。離散化した時刻の差分 ${\displaystyle h=t_{i+1}-t_{i}}$ を刻み幅あるいはステップサイズと呼ぶ^[17]。なお、ここでは時間ステップ ${\displaystyle h=t_{i+1}-t_{i}}$ は一定としたが、これを動的に決定する適応刻み（英語版）という手法もある^[18][19]。

厳密解 ${\displaystyle y(t)}$ から差分商

${\displaystyle \Delta (t,\eta ;h)={\begin{cases}{\frac {y(t+h)-\eta }{h}}&h\neq 0\\f(t,\eta )&h=0\end{cases}}}$

を導入するとき、数値積分スキーム ${\displaystyle \Phi }$ の局所離散化誤差^[20] (local discretization error) は

${\displaystyle \Delta (t,\eta ;h)-\Phi (t,\eta ;h)}$

により定義される^[21][22]。整合性のために ${\displaystyle h\to 0}$ の極限で局所離散化誤差（厳密にはその上限）は 0 に収束することが要求される^[23][24]。さらに、この極限で局所離散化誤差が ${\displaystyle \Delta (t,\eta ;h)-\Phi (t,\eta ;h)=O(h^{p})}$ を満足するとき、この積分スキーム ${\displaystyle \Phi }$ は ${\displaystyle p}$ 次精度であるという^[25][23]。

一方、${\displaystyle t}$ を固定して ${\displaystyle n\to \infty }$ とするとき、大域離散化誤差^[26] (global discretization error)

${\displaystyle e(t;h_{n})=\eta (t;h_{n})-y(t),\ \ h_{n}={\frac {t-t_{0}}{n}}}$

が 0 に収束するならば、その積分スキームは収束するという^[27][24]。${\displaystyle p}$ 次精度 (${\displaystyle p>0}$) の1段法を十分に滑らかな関数 ${\displaystyle f}$ に適用するとき、そのスキームは収束し大域離散化誤差は

${\displaystyle e(t;h_{n})=O(h_{n}^{p})}$

のように振る舞うことが保証されている、すなわち大域離散化誤差のオーダーは局所離散化誤差のオーダーに等しい^[28][29]。この結果はすべての整合的な一段法が ${\displaystyle h\to 0}$ で漸近的に安定であることを意味するものの、ただし現実的に可能な ${\displaystyle h}$ で1段法が安定であることは必ずしも保証されない（#安定性節を参照）^[30]。

オイラー法、ホイン法、古典的ルンゲ＝クッタ法 (RK4) の相対誤差の比較。初期値 ${\displaystyle y(0)=0}$ のもとでの常微分方程式 ${\displaystyle {\frac {dy}{dt}}=\cos(y)}$ の数値解の ${\displaystyle t=1}$ での値をステップサイズの関数として対数プロットした。各手法がそれぞれ1次精度、2次精度、4次精度であることに対応して、傾き 1, 2, 4 で誤差が減少している^[31]。

積分スキーム ${\displaystyle \Phi }$ としては以下のものが知られている。

• オイラー法^[32]: ${\displaystyle \Phi (t,\eta ;h)=f(t,\eta )}$ とするもの。これは1次精度のスキームである。
• ホイン法（英語版）^[33]: ${\displaystyle \Phi (t,\eta ;h)=[f(t,\eta )+f(t+h,\eta +hf(t,\eta ))]/2}$ とするもの。これは2次精度のスキームであり、関数 ${\displaystyle f}$ の評価を2回必要とする。
• 古典的ルンゲ＝クッタ法: これは4次精度のスキームであり、関数 ${\displaystyle f}$ の評価を4回必要とする^[34]。

このうち古典的ルンゲ＝クッタ法は適用可能範囲の広さ^[18]やプログラミングの容易さのために広く用いられている^[35]（ただし Press らは著書において計算速度の観点から古典的ルンゲ＝クッタ法に否定的に言及している^[36]）。

→「ルンゲ＝クッタ法」も参照

多段法

一段法は ${\displaystyle \eta _{i+1}}$ の値を ${\displaystyle \eta _{i}}$ だけから定めるものであったが、より多くのステップでの値 ${\displaystyle \eta _{i-1}}$, ..., ${\displaystyle \eta _{i-r+1}}$ を使う積分スキームは多段法^[15] (multistep method) と呼ばれる^[37]。この場合、最初の ${\displaystyle \eta _{0}}$, ..., ${\displaystyle \eta _{r-1}}$ はこのスキームでは定めることができず、1段法などの他の方法を用いる必要がある^[37]。多段法としてはアダムス・バッシュフォース法や、それを応用する予測子修正子法などがある^[38]（これはどちらも ${\displaystyle \eta _{i+1}}$ が ${\displaystyle \eta _{i-r+1}}$, ..., ${\displaystyle \eta _{i}}$ の線型関数として定まるため、線型多段法と呼ばれる^[39]）。

→詳細は「線型多段法」を参照

安定性

数値積分スキームの安定性はしばしばそれを初期値問題

${\displaystyle {\frac {dy}{dt}}=\lambda y,\ \ y(0)=1}$

（${\displaystyle \lambda }$ は複素数の定数）に適用することで定量化される^[40]。問題のスキームを一定の刻み幅 ${\displaystyle h}$ で適用することで得られる数列 ${\displaystyle \{y_{n}\}_{n=0,1,\cdots }}$ について、それが極限 ${\displaystyle n\to \infty }$ で 0 に収束するような ${\displaystyle h\lambda }$ がなす複素平面上の領域を絶対安定領域^[41]と呼ぶ^[40]。もとの初期値問題の厳密解は ${\displaystyle y(t)=e^{\lambda t}}$ であり、これは ${\displaystyle \mathrm {Re} \{\lambda \}<0}$ のとき ${\displaystyle t\to \infty }$ で 0 に収束する^[42]。そこである積分スキームの絶対安定領域が左半平面を含むとき、そのスキームはA-安定 (A-stable) であるという^[43][44]。

関数 ${\displaystyle f(t,y)}$ が大きな値を取り積分スキームによっては極端に時間刻み幅 ${\displaystyle h}$ を小さくする必要がある場合、その常微分方程式は硬い方程式と呼ばれる^[45]。この場合には十分に安定なスキームを用いる必要が生じる^[44]。これはしばしば ${\displaystyle y_{i+1}}$ を陰関数的に ${\displaystyle y_{i}}$ 等から定める陰解法によって実現される^[46]。

→詳細は「硬い方程式」を参照

幾何学的解法

解析対象となる微分方程式が何らかの特別な性質を持つとき、汎用の数値積分法ではなく、その性質を尊重するように構成された数値解法を用いることがあり、そのような手法は構造保存型解法または幾何学的数値解法と呼ばれる^[47][48]。例えば古典力学（ハミルトン力学）において時間発展はシンプレクティック写像であり、エネルギー（ハミルトニアン）等の保存量が存在する^[48]。この場合にはシンプレクティック数値積分法というハミルトン系にのみ適用可能な数値解法が存在し、良好な性質を持つことが知られている^[49]。

→詳細は「シンプレクティック数値積分法」を参照

境界値問題

常微分方程式

${\displaystyle {\frac {dy}{dx}}=f(x,y)}$

の解 ${\displaystyle y(x)}$ で, ${\displaystyle a\neq b}$ に対する境界条件

${\displaystyle r(y(a),y(b))=c}$

を満足するものを求める問題を境界値問題 (boundary-value problem) と呼ぶ^[50]。初期値問題と異なり、境界値問題では複数の解が存在すること、あるいは解が存在しないことがあり得る。また、スツルム＝リウヴィル型微分方程式のように常微分方程式にパラメータ ${\displaystyle \lambda }$ が含まれ、境界値問題に解が存在するようにパラメータ ${\displaystyle \lambda }$ を同時に定める問題もある^[51]。

これらの問題を数値的に解く最も単純な方法が狙い撃ち法 (shooting method) である^[52]。一方の端点（例えば ${\displaystyle x=a}$）において初期条件 ${\displaystyle y(a)}$ を適当に定めて微分方程式をもう一方の端点 ${\displaystyle x=b}$ まで解き、境界条件を満たすように未定の初期条件（および固有値）を適切に選ぶ^[53][54]。これはニュートン法などの関数の根を求めるアルゴリズム（これはしばしば逐次反復を伴う）を常微分方程式の数値解法と組み合わせることを意味する^[55][56]。ただし初期条件によっては区間 ${\displaystyle [a,b]}$ 全体で定義された解が存在しないことがあり、そのために改良された手法がある^[57]。あるいは有限要素法 (finite element method) などの手法も用いられる^[58]。

解の存在検証

→「精度保証付き数値計算」および「計算機援用証明」も参照

高精度に解く技術が追求されている一方で、「計算機で解の存在を検証する」という研究もおこなわれている^[59]。このような研究が必要となるのは、近似解が求まったとしてもそれが幻影解である危険性があるからである。偏微分方程式ではすでに幻影解が報告されているので^{[59][B 1][B 2]}、常微分方程式でも警戒が必要である。偏微分方程式の時と同様に関数解析学的な手法^{[59][60][61][62]}も考えられるが、関数解析学に頼らない手法 (例えば狙い撃ち法、スペクトル法やアフィン演算など) に基づく研究が主流であり^[59]、欧米などの海外^{[B 3][B 4][B 5][B 6][B 7][B 8][B 9][B 10][B 11]}のみならず日本国内でも研究されている^{[B 12][B 13][B 14][B 15]}。また、爆発解 (英: Blow-up solution) に特化した精度保証付き解法も探求されている^{[B 16][B 17]}。

解の存在検証・計算機援用証明が行われた方程式

• ローレンツ方程式^{[B 18][B 19][B 20][B 21][B 22][B 23]}
• レスラー方程式^{[B 24]}
• ダフィング方程式^{[63][B 25][B 26]}

関連ソフトウェア・ライブラリ

• COSY INFINITY^{[B 27][B 28]}
• INTLAB: MATLAB/GNU Octaveを使用する区間演算ライブラリであり、version 11 からODEソルバー (AWA toolbox と Taylor model toolbox の二つ) が搭載された^{[B 3][B 29]}。
• kv C++による区間演算ライブラリであり、ODEソルバーを搭載している。^{[B 3][B 30]}
  • kv - GitHub
• CAPD C++による区間演算ライブラリであり、力学系への応用のために開発された
• en:Chebfun