数値解析(概要)

一階の常微分方程式は, 一般に

\displaystyle \frac{d}{dx}u(x) = f(x, u(x))

の形で記述されます。 $t$は独立な変数, $f(x, u)$は既知の二変数関数であり, $u(x)$は求める未知関数です。

また, 複数の一階の常微分方程式の連立方程式を考えることもあります。

\displaystyle
  \begin{cases}
    \frac{d}{dx}u_i(x) = f_i(x, u_1(x), \ldots, u_m(x)),\,\,\,(i = 1, \ldots, m)
  \end{cases}

まずは, 1つの一階の常微分方程式の解軌跡を求める方法を紹介します。また, この解法は高階の微分方程式や高階の連立微分方程式の解放に拡張することができます[1][齊藤宣一, 数値解析入門 初版. 東京大学出版会, 東京, p209-p252(第8章), (2012)]。

***

本稿で紹介するのは, 離散変数法と呼ばれるもので,有限個の点に対して関数$u(x)$の値を予測するものです。まず, 関数$u$を予測する区間(定義域)$[a, b]$を$N$分割します。必ずしも等分割する必要はないですが, 今回は$[a, b]$を$N$等分割し, $h = \frac{b-a}{N},,x_i = a + ih,(i = 0, \ldots, N)$とします。

初期点$u(x_0=a) = u_a$を与えておき, この点の情報をもとに$u(x_1), u(x_2), ...$を逐次的に推測していく方法(い)と, 一度に$u(x_1)$から$u(x_N)$までの値を求める方法(ろ)があります。基本的に(ろ)は$N\times N$行列の逆行列を求める必要があるため, (い)と比べて計算機的に解くのが遅くなることがあります(その分, 精度の面では期待できます)。両端点の条件$u(a) = u_a, u(b) = u_b$が与えられている場合は(ろ)を用いて解くしかなさそうだと思われるかもしれませんが, (い)を工夫することで, 求める方法も提案されています(2)。

## 1階微分方程式これからいくつか手法を紹介しますが, 基本の考えはみな同じで

\displaystyle \frac{d}{dx}u(x) = f(x, u(x))

の両辺を極小区間$[x_i, x_{i+1}]$で積分することを考えます。左辺については, $\int_{x_i}^{x_{i+1}}\frac{d}{dx}u(x)dx = u(x_{i+1}) - u(x_i)$と計算できますが, 右辺 $\int_{x_i}^{x_{i+1}}f(x, u(x))dx$ をどう近似するか? というところで手法の違いが出てきます。

以下, 見やすくするために, $u_i = u(x_i)$と書くこともあります。

オイラー法(Eular)

\displaystyle \int_{x_i}^{x_{i+1}}f(x, u(x))dx \approx hf(x_i, u(x_i))

と近似します。

これより

u_{i+1} = u_i + hf(x_i, u_i)

という漸化式が得られます。これから初期点$u_0 = u_a$から逐次的に$u_N$まで求められることがわかりますね。またこれは前進オイラー法と呼ばれるもので, 端点$f(x, u_{i+1})$で面積を近似するのは後退オイラー法と呼ばれます。

u_{i+1} = u_i + hf(x, u_{i+1})

memo: 収束次数1

クランク-ニコルソン法(Crank-Nicolson)

\displaystyle \int_{x_i}^{x_{i+1}}f(x, u(x))dx \approx \frac{h}{2}\left(f(x_i, u(x_i)) + f(x_{i+1}, u(x_{i+1})\right)

と近似します。

いわゆる, 台形近似ですね。これより

u_{i+1} = u_i + \frac{h}{2}\left(f(x_i, u_i) + f(x_{i+1}, u_{i+1})\right)

という漸化式が得られます。右辺に$u_{i+1}$が出てきているので, 等式を満たす$u_{i+1}$を求める必要があります(後退オイラー法もそうです)。$g(t) = t - u_i - \frac{h}{2}\left(f(x_i, u_i) + f(x_{i+1}, t)\right)$とおいて, $g(t^*) = 0$を満たす, $t^*$ を見つけてそれを$u_{i+1}$とすることになります。もし, $f$の導関数が計算できるのであればニュートン-ラフソン法(Newton-Raphson)3を, そうでなければその他の球根アルゴリズム(二分法, 割線法, ブレント法,...)4を使用するのが良いでしょう。

オイラー法のように$u_{i+1}$を直接計算できるものは陽的な方法, クランク-ニコルソン法のように直接計算できないものは陰的な方法と呼ばれます。

ホイン法(Heun)

この手法では, クランク-ニコルソン法による後進式の右辺の$u_{i+1}$を$u_i+hf(x_i, u_i)$に置き換えることで陽的な方法にします。

u_{i+1} = u_i + \frac{h}{2}\left(f(x_i, u_i) + f(x_{i+1}, u_i+hf(x_i, u_i))\right)

これは

u'|_{x=x_i} \approx \frac{u_{i+1}-u_{i}}{h} = f(x_i, u_i)

を解くことで得られます。

memo: 収束次数2

### アダムス-バッシュフォース法(Adams-Bashforth) この手法では, クランク-ニコルソン法による後進式の右辺の$f(x_{i+1}, u_{i+1})$を$f(x_{i-1}, u_{i-1})$と$f(x_i, u_i)$から直線近似させて陽的な方法にします。

u_{i+1} = u_i + \frac{h}{2}\left(3f(x_i, u_i) - f(x_{i-1}, u_{i-1})\right)

$f(x_{i+1}, u_{i+1}) + f(x_{i-1}, u_{i-1}) = 2f(x_i, u_i)$から上の更新式が導かれています。更新式に$u_n$と$u_{n-1}$が含まれているため, 二段法の一つになります。

ルンゲ-クッタ法(Runge-Kutta)

まず積分値を

\displaystyle \int_{x_i}^{x_{i+1}}f(x, u(x))dx \approx \frac{h}{6}\left(f(x_i, u(x_i)) + 4f(x_{i+1/2}, u_{i+1/2}) + f(x_{i+1}, u(x_{i+1})\right)

と近似します。これはシンプソン則という方法です(d)。端点と中点の3点を通る二次関数を用いて面積を近似します。

これをホイン法のように陽的な方法に変換します。

ただ, $u_{i+1/2}$の近似には

\begin{align}
u'|_{x=x_{i+1/2}} &\approx \frac{u_{i+1/2}-u_{i}}{h} = f(x_i, u_i)\cdots(い)\\
u'|_{x=x_{i+1/2}} &\approx \frac{u_{i+1/2}-u_{i}}{h} = f(x_{i+1/2}, u_{i+1/2}) \cdots(ろ)
\end{align}

の2通りを用いることにします。更新式は次の通りです。

\begin{align}
k_1 &= f(x_i, u(x_i)\\
k_2 &= f(x_{i+1/2}, u_{i}+\frac{1}{2}hk_1)\cdots(い)\\
k_3 &= f(x_{i+1/2}, u_{i}+\frac{1}{2}hk_2)\cdots(ろ)\\
k_4 &= f(x_{i+1}, u_{i} + \frac{1}{2}hk_3)\\
u_{i+1} &= \frac{h}{6}(k_1+2k_2+2k_3+k_4)
\end{align}

この形のルンゲ-クッタ法は特に, 4段階4次精度の古典的ルンゲ-クッタ法と呼ばれます。

memo: 収束次数4

連立微分方程式の場合

\displaystyle
  \begin{cases}
    \frac{d}{dx}u_i(x) = f_i(x, u_1(x), \ldots, u_m(x)),(i = 1, \ldots, m)
  \end{cases}

連立方程式への拡張は自然な方法で行います。

u(x) = \begin{bmatrix} u_1(x)\\ u_2(x)\\ \vdots\\ u_m(x) \end{bmatrix}\,\,\,\,
f(x, u) = \begin{bmatrix} f_1(x, u_1(x))\\ f_2(x, u_2(x))\\ \vdots\\ f_m(x, u_m(x))\end{bmatrix}

とおいて, 初期点$u(x_0) = a$を用いて

$u(x_{n+1}) = u(x_n) + (方法ごとに定まる関数)(x, u(x_n), u(x_{n-1}))$

を適用すれば, 解軌跡を求められます。

高階微分方程式の場合

例えば, 2階微分方程式の場合

u''(x) = f(x, u(x), u'(x))

1. 1階の連立微分方程式に変換する

u_1 \leftarrow u, \,\,\,\,\, u_2 \leftarrow u'

として

  \begin{cases}
    u_1' = u_2\\
    u_2' = f(x, u_1, u_2)
  \end{cases}

という連立微分方程式に落とし込めます。当然, これは2階以上の高階微分方程式にも適用できます。

さらに関数$f$が$x$のみで表せられるときは次の方法も使えます。

2. $u''$や$u'$を$u$を使って表す
Taylor展開から分かるように
$u''(x)$は$\cfrac{u(x-h)-2u(x)+u(x+h)}{h^2}$と2次近似できるため

\displaystyle u''(x) = f(x) \iff u(x-h)-2u(x)+u(x+h) = h^2f(x)

とすることができます。これより

u_{n-1}-2u_n+u_{n+1} = h^2f(x_n)\,\,\,(n = 1, \ldots, N-1)

という1次線形連立方程式が得られます。
あとはこれに端点条件 $u_0 = u_a$ と $u_N = u_b$ を用いれば行列の逆行列を解くことで求めることができます。

\begin{bmatrix}
1  & 0 & 0 & \cdots & 0 & 0\\
-1 & 2 & 1 & \cdots & 0 & 0\\
0  & -1& 2 & \cdots & 0 & 0\\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 2 & -1 \\
0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix}
\begin{bmatrix} u_0 \\ u_1 \\ u_2 \\ \vdots\\ u_{N-1} \\ u_N \end{bmatrix}
= \begin{bmatrix} u_a \\ h^2f(x_1) \\ h^2f(x_2) \\ \vdots\\ h^2f(x_{N-1}) \\u_b \end{bmatrix}

左の行列は対角優位行列なので逆行列が必ず存在し, またガウスの消去法や反復法によって逆行列が求めやすいような性質の良い行列になっています6。

また終端条件$u_N = u_b$が与えられていない場合でも, オイラー法 $u_N = u_{N-1} + hf(x_{N-1})$などを用いれば計算することができます。

\begin{bmatrix}
1  & 0 & 0 & \cdots & 0 & 0\\
-1 & 2 & 1 & \cdots & 0 & 0\\
0  & -1& 2 & \cdots & 0 & 0\\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 2 & -1 \\
0 & 0 & 0 & \cdots & -1 & 1 \end{bmatrix}
\begin{bmatrix} u_0 \\ u_1 \\ u_2 \\ \vdots\\ u_{N-1} \\ u_N \end{bmatrix}
= \begin{bmatrix} u_a \\ h^2f(x_1) \\ h^2f(x_2) \\ \vdots\\ h^2f(x_{N-1}) \\hf(x_N) \end{bmatrix}

終端条件

$u'' = f(x, u, u'')$の形の微分方程式の場合, 上記のような特殊な場合を除いて端点条件を初期点と終点の両方を入れて数値解析で解を求めることは一般にはできません。
そこで初期点の情報として$u'_0$に制約がない場合, 終端条件を満たすような$u'_0$を見つけるという方法があります2。

この場合, $g(u'_0) = $(終点条件におけるx=bでの値) - (数値解析解におけるx=bでの値)とすると, $g$の根を見つける問題になるため, 何かしらの求根アルゴリズムにより根を求めることになります。

大抵の求根アルゴリズムではアルゴリズムの最初に$g(x)g(y) < 0$となるような$x, y$を求めておく必要がありますが, 割線法ではそれを求めておく必要がありません。この論文でも割線法を用いていますが, ブレント法などの方が収束速度が良いため, まずは割線法を用いて探索を進め, $g(x)g(y) < 0$となるような$x, y$が求められたら他の手法に切り替えるという方法も効果があると思われます。

参考

[1] 齊藤宣一, 数値解析入門 初版. 東京大学出版会, 東京, p209-p252(第8章), (2012)
[2] 釜国男, 変分問題の数値解法 (計算経済学の研究その9). https://soka.repo.nii.ac.jp/?action=repository_uri&item_id=35712&file_id=15&file_no=1
[3] ニュートン法, https://ja.wikipedia.org/wiki/ニュートン法
[4] 根の探索アルゴリズム, https://qiita.com/nariaki3551/items/8694eee8793f88e89f41
[5] シンプソンの公式, https://ja.wikipedia.org/wiki/シンプソンの公式
[6] 渡部善隆, 連立1次方程式の基礎知識, http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/MANUSCRIPT/KOHO/GEPP/GEPP.pdf