ODEの数値解法アプローチ

常微分方程式(Ordinary Differential Equation, ODE)の初期値問題(Initial Value Problem, IVP)に対する数値解法をアプローチごとに紹介する。網羅的ではないし、結果的に異なるアプローチから同じ数値解法が得られることもあるが、数値解法のアイデアを理解するために、アプローチごとに分類して説明する。

ODEの初期値問題の数値アプローチ

yを時間に対する関数とし、$y(t)\in\mathbb{R}^N$とする。微分方程式$\frac{dy}{dt} = f(t, y)$に対して初期値$\quad y(t_0) = y_0$が与えられたとき、$y(t)$を求める問題を初期値問題(IVP)というが、これの数値的な解法のアプローチを以下に挙げる。時間を離散化して, $t_n = t_0 + nh$とし$y_n(=y(t_n))$を使って$y_{n+1}$を求める。$h$は時間の刻み幅で、$10^{-3}$などの小さい値を取ることが多い。
以下、$y$は真の解を表し、$u$は数値解を表すものとする。

1. テイラー展開で近似

テイラー展開

$$y(t_{n+1}) = y_n + f(t_n, y_n)\underbrace{(t_{n+1} - t_n)}_{=h} + \frac{1}{2}f'(t_n, y_n)\underbrace{(t_{n+1} - t_n)^2}_{=h^2} + \cdots$$

を用いて、$y(t)$を近似する方法。高次の項を無視すると、微分方程式の右辺$f$の値やその微分を利用して、$y(t)$を近似的に求めることができる。

• 次数1で打ち切ると $u_{n+1} = u_n + h f(t_n, u_n).$ これはオイラー法と呼ばれる数値解法の一種である。単に$f(t_n, u_n)$を用いて、$u_{n+1}$を線形に近似しているので、実装やアルゴリズムは簡単だが、後に記述するように精度や安定性が低い。

• 次数2で近似すると $u_{n+1} = u_n + h f(t_n, u_n) + \frac{h^2}{2}f'(t_n, u_n)$ となり、

  $$f'(t_n, u_n) = \frac{\partial f}{\partial t}(t_n, u_n) + \frac{\partial f}{\partial u}(t_n, u_n)\underbrace{\frac{du}{dt}}_{=f(t_n, u_n)} \quad \text{(連鎖率)}$$

  を計算する必要がある。ここで、$Df(t_0, u_0)$は$f$の$u$に関するヤコビ行列である。まとめると

  $$u_{n+1} \approx u_n + h f(t_n, u_n) + \frac{h^2}{2}\left(\frac{\partial f}{\partial t}(t_n, u_n) + Df(t_n, u_n)f(t_n, u_n)\right)$$

  だが、ヤコビ行列の計算が計算コストが高いことが多いので、ポピュラーな方法ではない。

2. 積分方程式を求積法で近似

微分方程式を積分方程式に書き直すと

$$y_{n+1} = y_n + \int_{t_n}^{t_{n+1}} f(t, y(t)) dt$$

となる。この積分を求積法で近似する方法。

• Forward Euler法 $$\int^b_a g(t) dt \approx (b-a)g(a)$$ を使うと, $$u_{n+1} = u_n + h f(t_n, u_n)$$ で先程テイラー展開で1次近似したオイラー法と同じ式が得られる。
• Backward Euler法 $$\int^b_a g(t) dt \approx (b-a)g(b)$$ を使うと, $u_{n+1} = u_n + h f(t_{n+1}, u_{n+1})$ となる。これも一次近似であり、$u_{n+1}$が両辺に現れるため、非線形方程式を解く必要がある。
  このように、右辺に未知の値$u_{n+1}$が現れる方法を 陰解法(implicit method) と呼ぶ。Forward Eulerのような、右辺に未知の値が現れない方法を 陽解法(explicit method) と呼ぶ。陰解法は安定性が高いが、非線形方程式を解く必要があるため、計算コストが高いことが多い。 例えばNewton法を使って $$F(u_{n+1}) = u_{n+1} - u_n - h f(t_{n+1}, u_{n+1}) = 0$$ を解く。安定性が高いが、計算コストも高い。
• 台形則 台形則 $$\int^b_a g(t) dt \approx \frac{b-a}{2}(g(a) + g(b))$$ を使うと, $u_{n+1} = u_n + \frac{h}{2}(f(t_n, u_n) + f(t_{n+1}, u_{n+1}))$ となる。これは2次近似で、$u_{n+1}$が両辺に現れるため、Newton法などで非線形方程式を解く必要があるが、安定性が高い(ただしBackward Euler法よりは低い)。
• 中点法 $$\int^b_a g(t) dt \approx (b-a)g\left(\frac{a+b}{2}\right)$$ を使うと, $u_{n+1} = u_n + h f\left(t_n + \frac{h}{2}, \underbrace{u(t_n + \frac{h}{2})}_{?}\right)$ となる。$u(t_n + \frac{h}{2})$は未知なので、例えば$\frac{u_n + u_{n+1}}{2}$などで置き換える。したがって $$u_{n+1} = u_n + h f\left(t_n + \frac{h}{2}, \frac{u_n + u_{n+1}}{2}\right)$$ となり、これも非線形方程式を解く必要があるが、安定性が高い(ただしBackward Euler法よりは低い)。

3. 予測子修正子法

先の方法は, 右辺にも$u_{n+1}$が現れるため、非線形方程式を解く必要があったが、予測子修正子法は右辺の$u_{n+1}$を$u_n$などの既知の値を使ったEuler法等で予測し(予測子)、それを陰的な更新式(修正子)に代入することで、非線形方程式を解く必要をなくしたものである。例えば以下のようなものがある。

Euler予測子を用いた中点法

Euler法で中点の点を

$$u_{n+\frac{1}{2}} = u_n + \frac{h}{2}f(t_n, u_n)$$

で近似することで、更新式が

$$u_{n+1} = u_n + h f\left(t_n + \frac{h}{2}, u_{n+\frac{1}{2}}\right)$$

とexplicitになる。

Euler予測子を用いた台形則

Euler法で右辺の$u_{n+1}$を

$$u^*_{n+1} = u_n + h f(t_n, u_n)$$

と近似することで、更新式が

$$u_{n+1} = u_n + \frac{h}{2}(f(t_n, u_n) + f(t_{n+1}, u^*_{n+1}))$$

とexplicitになる。

4段4次Runge-Kutta法

求積法のSimpson則

$$\int^b_a g(t) dt \approx \frac{b-a}{6}(g(a) + 4g\left(\frac{a+b}{2}\right) + g(b))$$

を用いて、ステップ$n+1/2$と$n+1$の両方を予測子で近似してexplicitにする。更新式は

$$\begin{aligned} k_1 &= f(t_n, u_n) \\ k_2 &= f\left(t_n + \frac{h}{2}, u_n + \frac{h}{2}k_1\right) \quad \left(\approx f\left(u\left(t_n + \frac{h}{2}\right)\right)\right)\\ k_3 &= f\left(t_n + \frac{h}{2}, u_n + \frac{h}{2}k_2\right) \quad \left(\approx f\left(u\left(t_n + \frac{h}{2}\right)\right)\right)\\ k_4 &= f(t_n + h, u_n + hk_3) \quad \left(\approx f\left(u\left(t_n + h\right)\right)\right)\\ u_{n+1} &= u_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}$$

となる。これは4次近似で、安定性も高いが、計算コストも今までの低次の方法より高い。

4. 格子点の多項式補間で近似

上までは常に$u_n$を使って$u_{n+1}$を求める方法であった(中間の点を補間することで次数を上げてた)。しかし、$u_{n-1}, u_{n-2},\cdots$などのより過去の点を使った多段的な解法が存在し、 ここでは多項式補間を用いて導出される代表的な線形多段法として、Adams系とBDFを紹介する。

Adams系

Adams系は$t_n, t_{n-1}, \cdots$における微分方程式の右辺$f$の値$f_n, f_{n-1}, \cdots$を使って$f$を補間する多項式$p(t)$を作り、

$$y_{n+1} = y_n + \int_{t_n}^{t_{n+1}} p(t) dt$$

と積分する。陽解法がAdams-Bashforth法、$f_{n+1}$を使う陰解法がAdams-Moulton法である。

• Adams-Bashforth法は、Newton補間多項式を使う。まずは差分商を定義すると、 $$\begin{aligned} f[t_n, t_{n-1}] &= \frac{f_n - f_{n-1}}{t_n - t_{n-1}} \\ f[t_n, t_{n-1}, t_{n-2}] &= \frac{f[t_n, t_{n-1}] - f[t_{n-1}, t_{n-2}]}{t_n - t_{n-2}}. \end{aligned}$$ つまり差分商$f[t_n, t_{n-1}]$は ($f$の変化量) / ($t$の変化量)、$f[t_n, t_{n-1}, t_{n-2}]$は$f[t_n, t_{n-1}]$と$f[t_{n-1}, t_{n-2}]$の差分商である。(それぞれ、幾何的には傾きと曲率に対応する。) これを使って、$f$のNewton補間多項式を作ると¹ $$\begin{aligned} p_k(t) = f_n + f[t_n, t_{n-1}](t - t_n) + f[t_n, t_{n-1}, t_{n-2}](t - t_n)(t - t_{n-1}) + \cdots\\ + f[t_n, t_{n-1}, \cdots, t_{n-k+1}](t - t_n)(t - t_{n-1})\cdots(t - t_{n-k+1}). \end{aligned}$$ これを、step幅が定数$h$であることを使って積分すると、$u_{n+1}$を求める式が得られる。 例えば、
  • 2段 ($k=2$) の場合 $$u_{n+1} = u_n + \frac{h}{2}(3f_n - f_{n-1}).$$
  • 3段 ($k=3$) の場合 $$u_{n+1} = u_n + \frac{h}{12}(23f_n - 16f_{n-1} + 5f_{n-2}).$$
• Adams-Moulton法は、$f_{n+1}$を使う陰解法である。多項式は $$\begin{aligned} p_k(t) = f_{n+1} + f[t_{n+1}, t_n](t - t_{n+1}) + f[t_{n+1}, t_n, t_{n-1}](t - t_{n+1})(t - t_n) + \cdots\\ + f[t_{n+1}, t_n, \cdots, t_{n-k+1}](t - t_{n+1})(t - t_n)\cdots(t - t_{n-k+2}). \end{aligned}$$ これを刻み幅$h$として積分すると、$u_{n+1}$を求める式が得られる。 例えば、
  • 2段 ($k=2$) の場合 $$u_{n+1} = u_n + \frac{h}{12}(5f_{n+1} + 8f_n - f_{n-1}).$$
  • 3段 ($k=3$) の場合 $$u_{n+1} = u_n + \frac{h}{24}(9f_{n+1} + 19f_n - 5f_{n-1} + f_{n-2}).$$

Backward Differentiation Formula, BDF法 (後退差分法)

これも多項式補間だが、微分方程式の右辺の$f$を補間するのではなく、数値解$u(t)$を補間する($u_{n+1}$を節点に含める)。 補間した多項式$p(t)$を微分して、$t_{n+1}$で評価する、つまり

$$p'(t_{n+1}) = f(t_{n+1}, u_{n+1})$$

という方程式を作り、これを解くことで$u_{n+1}$を求める。 例えば2段($k=2$)のBDF2の場合、$u_{n+1}, u_n, u_{n-1}$を通る多項式$p(t)$を作り、これを微分して

$$p'(t_{n+1}) = \frac{3u_{n+1} - 4u_n + u_{n-1}}{2h}$$

を得る。したがって

$$\frac{3u_{n+1} - 4u_n + u_{n-1}}{2h} = f(t_{n+1}, u_{n+1})$$

を解くことで$u_{n+1}$を求める。BDF法は安定性が高いが、非線形方程式を解く必要がある。

5. 未定係数法

Adams系やBDFを含む線形多段法 (Linear Multistep Method, LMM) は、一般に

$$u_{n+1} + \alpha_0 u_n + \cdots + \alpha_{k-1} u_{n-k+1} = h (\beta_{-1} f_{n+1} + \beta_0 f_n + \cdots + \beta_{k-1} f_{n-k+1})$$

の形で表される。ここで$\alpha_i, \beta_i$は各メソッド固有の係数である。 Approach 4では、補間多項式を積分 or 微分して係数を求めたが、未定係数法では、まずこの形を仮定し、$u(t)$をテイラー展開して、$u_{n+1}$を$t_n$の周りで展開することで、係数$\alpha_i, \beta_i$を決定する。

手法の分類

以上でアプローチを挙げたが、複数の観点からそれらをグループ分けできる。

1. 陽解法と陰解法:
   • 陽解法: $u_{n+1}= \ldots$ の右辺に$u_{n+1}$が現れない陽解法は計算コストが低いが安定性が低い。
   • 陰解法: 方程式をNewton法などで毎回解く必要があるので計算コストが高いが安定性が高い。
2. 単段法と多段法:
   • 単段法: $u_{n+1}$を求めるのに、$u_n$のみを使う。安定性が高い。
   • 多段法: $u_{n+1}$を求めるのに、さらに$u_{n-1}, u_{n-2}, \cdots$などの過去の点を使う。効率が良いが安定性に欠ける。余談だが、多段法は最初の数ステップは単段法で計算する必要がある。
3. 線形と非線形:
   • 線形: $u_{n+1}$が過去の$u$の線形結合で表される。Adams系やBDF法といった線形多段法など。
   • 非線形: 予測子を使うRunge-Kutta法などは$f$の値を$f$の引数に入れるため、一般に非線形である。

References

• Maria Cameron, Numerical Methods for Solving Ordinary Differential Equations, AMSC 661 Lecture Notes, University of Maryland.

Footnotes

1. Lagrange補間多項式も同値だが、点数$k$が増えても、Lagrangeの場合は係数を全て更新する必要があるが、Newtonの場合は既存の補間多項式に新しい高次の項を追加するだけで済む。 ↩