常微分方程式の数値解法1

今回の講義は常微分方程式（の初期値問題） \begin{equation} \frac{\mathrm d}{\mathrm dt} y(t) = f(y(t)), \quad y(0)=y_0 \in \mathbb R^d \end{equation} の数値計算についてです．（常微分方程式の境界値問題というものもありますが，これはむしろ偏微分方程式の特殊ケースとみなす方が見通しがよく，本講義では扱いません）． ここで，\(f:\mathbb R^d\to \mathbb R^d\) は与えられた（滑らかな）関数，\(y_0\) は与えられたベクトル（これを初期値と言います）です． 常微分方程式の厳密解は形式的には \begin{equation} \tag{1} \label{eq:ode1} y(t) = y_0 + \int_0^t f(y(s))\,\mathrm d s \end{equation} と書けますが，非常に特殊な場合を除いて右辺の積分を解析的に求めることはできません． 解を具体的に書くことができなくても，数学的には，解の存在・定性的性質など考えるべきことはたくさんあります． 一方で，物理など現代科学の多くの分野では，微分方程式の解を具体的に数値として構成する必要があります．そこで，常微分方程式の近似（数値）解を構成する手法が必要になります． 特に，離散化した時刻\( 0=t_0 < t_1< t_2 <\cdots \)における近似解を組織的に生成する方法を「離散変数法」といいます．

一言で離散変数法（離散変分法と書いていたのを修正しました；「離散変分法」は偏微分方程式の数値解法の一つで，阪大の降籏大介先生が創始したものです）といっても実に多様な方法がありますが，一つのアイデアは，\eqref{eq:ode1}の右辺の積分を近似するというものです．時間刻み幅を\(h\)と表すことにして，\(y(h),y(2h),y(3h),\dots\)の近似解を構成することを考えましょう．以後， \(y(nh)\) の近似解を \(y_n\) と表します．

例えば，数値積分を扱った以前の講義のAlgorithm1を用いると \[ y(h)=y_0 + \int_0^h f(y(s)) \,\mathrm ds \approx y_0 + h f(y_0) \] ですから，\( y(h)\) の近似 \( y_1\) を \[ y_1 = y_0 + h f(y_0)\] によって構成することが考えられます．この方法は \[\frac{y_1-y_0}{h} = f(y_0) \] とみれば，もとの常微分方程式そのものの近似と解釈することもできます． さらに，\( y(nh)\) の近似 \( y_n \) は \( y_0\) と \(y_1\) をそれぞれ \(y_n\) と \(y_{n+1}\) に置き換えて \[ y_{n+1} = y_n + h f(y_n), \quad n = 0,1,2,\dots\] によって構成することにしましょう．この方法を陽的Euler (オイラー) 法といいます．

積分の近似を変えると，例えば \[ y_1 = y_0 + h f(y_1)\] といった解法が考えられます．これを陰的Euler法といいます．\(f\) の中が \(y_1\) になっていますから，\(y_0\) から \(y_1\) を計算するときに， 陽的Euler法とは異なり，方程式を解かなければなりません． このように，先の時間の近似解を構成するために方程式を解かなければならない解法を 「陰的」な解法といい，その必要がない解法を「陽的」な解法といいます． \(f\) が線形でかつ \(d=1\) の場合は，陰的であってもただの（一変数の）１次方程式ですから，それを解くのは極めて簡単ですが, \(d>2\) の場合は，連立１次方程式になり，連立１次方程式の数値解法（アルゴリズム）は自明ではありません（クラメールの公式は \(d\) に対して計算量が指数的に増えるので実用上は使えません．また，「逆行列を使えばよいのでは？」と思うかもしれませんが，一回生のときの線形代数を思い出してもらうと分かるように，逆行列は複数の連立一次方程式を解いてはじめて求められるので，通常，連立一次方程式を解くために逆行列を用いることはありません．）．\(f\) が非線形の場合はさらに大変で，陽的な解法と比べて非常に大きな計算コストがかかることがあります．

天下り的ですが， \begin{align} Y_1 &= y_0 \cr Y_2 &= y_0 + \frac{h}{2} f(Y_1) \cr Y_3 &= y_0 + \frac{h}{2} f(Y_2) \cr Y_4 &= y_0 + h f(Y_3) \cr y_1 &= y_0 + \frac{h}{6} \big( f(Y_1) + 2f(Y_2) + 2f(Y_3) + f(Y_4) \big) \end{align} という解法も有名で，これを（4段4次の）陽的Runge–Kutta (ルンゲ・クッタ；以下RK) 法といいます． 4つの中間変数 \(Y_1,Y_2,Y_3,Y_4\) 用いることから4段の解法といいます（陽解法の場合，全体の計算コストはおおよそ関数 \(f\) の評価回数に比例します． その意味でRK法は陽的Euler法のおおよそ4倍のコストがかかります）．「4次」の意味については後述します．

解法を暗記する必要は全くありませんが，陽的Euler法，陰的Euler法，RK法は（前回のFourier級数と同様に）現代科学のリテラシーとして（少なくともその存在は）知っておくべき解法です．

（どのような解法の場合も）一般に \( y_1 \neq y(h)\) ですが，\(y_1\) がより \(y(h)\) に近い解法ほど精度のよい解法と言ってよいでしょう． ここで \(y(h)\) を形式的にTaylor展開してみます： \begin{equation} y(h) = y(0) + h \dot{y}(0) + \frac{h^2}{2} \ddot{y}(0) + \cdots = y_0 + h f(y_0) + \frac{h^2}{2} f^\prime (y_0) f(y_0) \cdots. \end{equation} これと陽的Euler法を比べてみると，\(h\) に関して1次の項まで一致していることが分かります． 一般に，厳密解 \(y(h)\) のTaylor展開と近似解 \(y_1\) のTaylor展開（\(y_1\) を \(h\) の関数とみなして \(h\) に関してTaylor展開）が \(h\) の \(p\) 次の項まで一致しているとき，その解法を「 \(p\) 次精度の解法」といいます． \(p\) が大きいほど， \(h\) を小さくしたときにより早く厳密解に収束しますし，ある固定の \(h\) で異なる解法を比較したとき， \(p\) が大きい方がより厳密解に近い数値解になることが多いです． しかし，陽的Euler法とRK法の比較から分かるように，通常， \(p\) を大きくするにはより大きな計算コストが必要になり，その意味で，精度とコストの間にはトレードオフがあります． 実際の応用で適切な数値解法を選択するためには，解きたい微分方程式の数学的性質と数値解法の様々な性質をよく理解しておくことが重要です． 陽的Euler法と陰的Euler法は1次精度解法で，RK法は4次精度解法です．

詳細は述べませんが，計算コストと精度の他にも，例えば「安定性」といった概念があり，一般に陽解法に比べて陰解法の方が安定に計算できることが多いです．ただし，この講義では，プログラムを書く際には，（簡単に解ける場合を除いて）陰解法は扱いません（連立一次方程式などの解き方も理解する必要があるため．．．）．

例題として，\(d=1\) の場合の次の微分方程式を考えましょう（\(d>1\) の場合は来週扱います）： \begin{equation} \frac{\mathrm d}{\mathrm dt}y = -y, \quad y(0) = 1 \end{equation} （この方程式の厳密解は簡単に書けて， \( y(t) = \mathrm{e}^{-t}\) です）．

では，プログラムを見て課題を行っていきましょう．

サンプルファイルはこちら