オイラー法

Euler's Method

このページの目標

前進オイラー法の漸化式と接線近似の幾何学的意味を理解し、誤差のオーダー安定性陰的方法の概念を学ぶ。

1. 問題設定

初期値問題(IVP: Initial Value Problem)を考える:

\begin{equation} \frac{dy}{dx} = f(x, y), \qquad y(x_0) = y_0 \label{eq:ivp} \end{equation}

両辺を $x_n$ から $x_{n+1}$ まで積分すると、

$$y(x_{n+1}) = y(x_n) + \int_{x_n}^{x_{n+1}} f\bigl(x,\, y(x)\bigr)\, dx$$

となる。右辺の積分を正確に求められれば厳密解が得られるが、一般には $f$ の中に未知関数 $y(x)$ が含まれるため解析的に計算できない。そこで、この積分を何らかの方法で近似することが数値解法の基本的な考え方である。

解析的に解を求められない場合、数値的に近似解を構成する必要がある。区間 $[x_0, x_0 + L]$ を等間隔に分割し

$$x_n = x_0 + nh, \quad n = 0, 1, 2, \dots, N, \qquad h = \frac{L}{N}$$

として、各格子点 $x_n$ における $y(x_n)$ の近似値 $y_n$ を逐次的に求める。

区間 [x0, x0+L] の等間隔分割と格子点の配置
図1: 区間 $[x_0, x_0 + L]$ を刻み幅 $h$ で $N$ 等分し、格子点 $x_0, x_1, \dots, x_N$ を配置する

2. 前進オイラー法(陽的オイラー法)

定義:前進オイラー法

初期値 $y_0$ から出発し、以下の漸化式で $y_1, y_2, \dots$ を順次計算する:

$$y_{n+1} = y_n + h \cdot f(x_n, y_n)$$

2.1 導出

テイラー展開から導出できる。厳密解 $y(x)$ に対して

$$y(x_{n+1}) = y(x_n) + h \, y'(x_n) + \frac{h^2}{2} y''(x_n) + O(h^3) = y(x_n) + h \, f(x_n, y(x_n)) + O(h^2)$$

$O(h^2)$ の項を無視すれば前進オイラー法の漸化式が得られる。

積分の観点からも同じ結果が得られる。第1節の積分形において、被積分関数 $f(x, y(x))$ を区間 $[x_n, x_{n+1}]$ で左端の値 $f(x_n, y_n)$ で一定と近似する(左端矩形公式)と、

$$\int_{x_n}^{x_{n+1}} f(x,\, y(x))\, dx \approx h \cdot f(x_n, y_n)$$

となり、前進オイラー法の漸化式が得られる。つまり前進オイラー法は、$f$ を各小区間で階段状に近似して積分していることに相当する。

左端矩形公式による階段近似の図。厳密解 y(x) の曲線と、各区間で左端の値を一定とする階段状の近似を重ねて表示
図2: 左端矩形公式による階段近似 — 各区間 $[x_n, x_{n+1}]$ で $f$ を左端の値 $f(x_n, y_n)$ で一定とみなし、オレンジの矩形の面積で積分を近似する

2.2 幾何学的解釈

上の積分の観点を幾何学的に言い換えると、次のようになる。式 $\eqref{eq:ivp}$ より $f(x, y) = dy/dx$ であったから、$f(x_n, y_n)$ は点 $(x_n, y_n)$ における解曲線の傾き、すなわち接線の傾きそのものである。したがって「$f$ を左端の値で一定とみなす」とは、解曲線を接線で置き換えて直進することにほかならない。

つまり前進オイラー法は、点 $(x_n, y_n)$ における接線に沿って刻み幅 $h$ だけ進み、到達した点を次の近似値 $y_{n+1}$ とする方法である。$x$ 方向に $h$ 進むと $y$ 方向に $h \cdot f(x_n, y_n)$ だけ変化する。

前進オイラー法の接線近似の図解
図3: 前進オイラー法の接線近似($y' = y$, $y(0)=1$, $h=0.5$)— 各点で接線に沿って進むが、ステップごとに誤差が蓄積する

図3で注意すべき点がある。青い各線分の傾きは、厳密解(緑の曲線)の傾きではなく、近似値 $y_n$ から計算した傾き $f(x_n, y_n)$ である。最初の区間 $[x_0, x_1]$ では $y_0 = y(x_0)$ なので両者は一致するが、2区間目以降は近似値 $y_n$ が厳密解 $y(x_n)$ からずれているため、青い線分の傾きも厳密解の傾きとは異なる。この「ずれた傾きでさらに進む」ことが誤差の蓄積を引き起こす原因である。

2.3 計算例

$y' = y$, $y(0) = 1$(厳密解 $y = e^x$)を $h = 0.5$ で $x = 2$ まで計算する。

$n$ $x_n$ $y_n$(オイラー) $y(x_n)$(厳密) 誤差
0 0.0 1.0000 1.0000 0.0000
1 0.5 1.5000 1.6487 0.1487
2 1.0 2.2500 2.7183 0.4683
3 1.5 3.3750 4.4817 1.1067
4 2.0 5.0625 7.3891 2.3266

$h = 0.5$ では誤差が大きいが、$h$ を小さくすれば改善する。ステップ幅を変えたときの $x = 2$ での誤差を比較すると:

ステップ幅 $h$ と $x=2$ での誤差の関係
$h$ ステップ数 $y_N$(オイラー) 誤差
0.545.06252.3266
0.1206.7275$6.6 \times 10^{-1}$
0.012007.3160$7.3 \times 10^{-2}$
0.0012,0007.3817$7.4 \times 10^{-3}$
0.000120,0007.3883$7.4 \times 10^{-4}$

$h$ を10分の1にすると誤差もほぼ10分の1になる。これはオイラー法が1次精度($O(h)$)であることを反映している。

$h$ は小さいほどよいか? 打ち切り誤差の観点では $h$ を小さくするほど精度が上がるが、実際には限界がある。

  • 丸め誤差の蓄積:$h$ が極端に小さいとステップ数が膨大になり、浮動小数点の丸め誤差が蓄積する。ある点で打ち切り誤差の減少と丸め誤差の増大が釣り合い、それ以上 $h$ を小さくしても精度は改善しない。
  • 計算時間:$h$ を半分にすると計算量は2倍になる。精度と計算コストのバランスが重要である。

高い精度が必要な場合は、$h$ をむやみに小さくするのではなく、ルンゲ=クッタ法のような高次精度の方法を使う方が効率的である。

3. 誤差解析

3.1 局所打ち切り誤差

局所打ち切り誤差(local truncation error)とは、1ステップで導入される誤差であり、厳密解からの出発を仮定して定義する:

$$\tau_n = y(x_{n+1}) - \bigl[y(x_n) + h f(x_n, y(x_n))\bigr] = \frac{h^2}{2} y''(\xi_n) = O(h^2)$$

テイラー展開の $O(h^2)$ の項が打ち切られるため、局所打ち切り誤差は $O(h^2)$ である。

3.2 大域誤差

大域誤差(global error)は $e_n = y(x_n) - y_n$ で定義される。$N = L/h$ ステップの蓄積により

$$|e_N| \le C \cdot h$$

すなわち大域誤差は $O(h)$ である。これをオイラー法は1次精度の方法であるという。直観的には、$O(h^2)$ の局所誤差が $N = O(1/h)$ ステップ蓄積して $O(h^2) \times O(1/h) = O(h)$ となる。

精度を $p$ 次にするには、1ステップの局所打ち切り誤差が $O(h^{p+1})$ である必要がある。

4. 後退オイラー法(陰的オイラー法)

定義:後退オイラー法

$y_{n+1} = y_n + h \cdot f(x_{n+1}, y_{n+1})$

積分の観点では、前進オイラー法が左端矩形公式(区間の左端の値で $f$ を一定と近似)であったのに対し、後退オイラー法は右端矩形公式(区間の右端の値 $f(x_{n+1}, y_{n+1})$ で一定と近似)に相当する:

$$\int_{x_n}^{x_{n+1}} f(x,\, y(x))\, dx \approx h \cdot f(x_{n+1}, y_{n+1})$$
右端矩形公式による階段近似の図。厳密解 y(x) の曲線と、各区間で右端の値を一定とする階段状の近似を重ねて表示
図4: 右端矩形公式による階段近似 — 各区間 $[x_n, x_{n+1}]$ で $f$ を右端の値 $f(x_{n+1}, y_{n+1})$ で一定とみなし、紫の矩形の面積で積分を近似する

右辺に未知の $y_{n+1}$ が現れるため、各ステップで非線形方程式を解く必要がある(通常はニュートン法を用いる)。計算コストは前進オイラー法より高いが、安定性に優れる。

4.1 なぜ安定性に優れるのか?

直観的な理由を、減衰する解 $y' = \lambda y$($\lambda < 0$)で見てみよう。厳密解は $y = e^{\lambda x} \to 0$ と減衰する。

前進オイラー法は「古い情報」$y_n$ に基づいて次の値を決める:

$$y_{n+1} = (1 + h\lambda)\, y_n$$

$\lambda < 0$ のとき、$h$ が大きすぎると $|1 + h\lambda| > 1$ となり、減衰するはずの解が振動・発散してしまう。例えば $\lambda = -10$ で $h = 0.3$ とすると $1 + h\lambda = 1 - 3 = -2$ であり、$|y_n|$ は毎ステップ2倍に増大する。

後退オイラー法は「行き先の情報」$y_{n+1}$ を使って自分自身を決める:

$$y_{n+1} = y_n + h\lambda\, y_{n+1} \implies y_{n+1} = \frac{y_n}{1 - h\lambda}$$

$\lambda < 0$, $h > 0$ のとき、$1 - h\lambda = 1 + h|\lambda| > 1$ であるから、$|y_{n+1}| < |y_n|$ が常に成り立つ。刻み幅 $h$ をどんなに大きくしても解は減衰し、発散しない。これが後退オイラー法の安定性の本質である。

言い換えれば、前進オイラー法は急激に変化する場面で「古い傾き」のまま大きく飛び出してしまうが、後退オイラー法は「到着点での傾き」で自動的にブレーキがかかる。この性質が特に重要になる硬い(stiff)方程式については第5.3節で述べる。

5. 安定性解析

安定性を調べるためにテスト方程式 $y' = \lambda y$($\mathrm{Re}(\lambda) < 0$)を用いる。厳密解は $y = e^{\lambda x} \to 0$($x \to \infty$)であり、数値解もこの減衰を再現することが望ましい。

5.1 前進オイラー法

$y_{n+1} = (1 + h\lambda) y_n$ であるから、増幅因子は $R(z) = 1 + z$($z = h\lambda$)。$|y_n| \to 0$ となるためには

$$|1 + z| < 1$$

が必要である。これは複素平面で中心 $(-1, 0)$、半径 $1$ の円の内部であり、安定領域と呼ばれる。$\lambda$ が実の大きな負の値の場合、$h < 2/|\lambda|$ という刻み幅の制限を受ける。

5.2 後退オイラー法

$y_{n+1} = \frac{1}{1 - z} y_n$ であるから、安定条件は $|1/(1-z)| < 1$ すなわち $|1 - z| > 1$ である。$\mathrm{Re}(\lambda) < 0$ ならばこの条件は $h > 0$ のすべてで満たされ、後退オイラー法はA安定(任意の刻み幅で安定)である。

5.3 硬い(stiff)方程式と後退オイラー法の利点

時定数の大きく異なる成分が混在する方程式を硬い(stiff)方程式と呼ぶ。例えば $y' = -1000y$ では $\lambda = -1000$ であり、前進オイラー法では $h < 2/|\lambda| = 0.002$ でないと発散してしまう。一方、後退オイラー法はどんな大きな $h$ でも安定に計算できる。計算コストは1ステップあたり高くても、刻み幅を大きく取れるため総ステップ数を大幅に削減でき、結果的に効率が良くなる。これが硬い方程式に対する後退オイラー法(および陰的方法一般)の利点である。

6. 改良オイラー法(ホイン法)

定義:改良オイラー法(ホイン法)

以下の2段階で $y_{n+1}$ を計算する:

$$\tilde{y}_{n+1} = y_n + h \, f(x_n, y_n) \qquad \text{(予測子:前進オイラー)}$$ $$y_{n+1} = y_n + \frac{h}{2} \bigl[f(x_n, y_n) + f(x_{n+1}, \tilde{y}_{n+1})\bigr] \qquad \text{(修正子:台形則)}$$

積分の観点では、$f$ を左端と右端の平均で近似する台形公式に相当する:

$$\int_{x_n}^{x_{n+1}} f(x,\, y(x))\, dx \approx \frac{h}{2}\bigl[f(x_n, y_n) + f(x_{n+1}, y_{n+1})\bigr]$$
台形公式による近似の図。厳密解 y(x) の曲線と、各区間で左端と右端を直線で結ぶ台形近似を重ねて表示
図5: 台形公式による近似 — 各区間 $[x_n, x_{n+1}]$ で $f$ を左端と右端の値の平均で近似し、緑の台形の面積で積分を近似する

矩形近似より曲線によく沿うため、精度が向上する:

  • 局所打ち切り誤差:$O(h^3)$
  • 大域誤差:$O(h^2)$(2次精度

改良オイラー法は2段2次のルンゲ=クッタ法(RK2)の一種であり、陽的方法でありながら前進オイラー法の2倍の精度次数を持つ。計算コストは1ステップあたり $f$ の評価が2回(前進オイラー法の2倍)であるが、同じ精度を達成するために必要な刻み幅を大幅に粗くできるため、総合的な効率は改良オイラー法が優れることが多い。

7. 各方法の比較

ここに述べた3つの陽的方法(前進オイラー法、改良オイラー法=ホイン法)と、別ページで解説する ルンゲ=クッタ法 の特性を比較すると次表のようになる。

方法 局所誤差 大域誤差 $f$ 評価/step A安定
前進オイラー $O(h^2)$ $O(h)$ 1 No
後退オイラー $O(h^2)$ $O(h)$ 1 + 求解 Yes
改良オイラー $O(h^3)$ $O(h^2)$ 2 No
ルンゲ=クッタ法 $O(h^5)$ $O(h^4)$ 4 No

上の表にある4つの方法について、刻み幅 $h$ と大域誤差の関係を対数スケールで比較すると、各手法の精度次数の違いが明確にわかる。

刻み幅hと大域誤差の関係(対数スケール)— 4手法の比較
図6: 刻み幅 $h$ と大域誤差の関係($y' = y$, $y(0)=1$, $x=2$ での誤差を実測)— 前進・後退オイラー法はともに1次精度 $O(h)$ で傾きが近いが、改良オイラー法は2次 $O(h^2)$、ルンゲ=クッタ法は4次 $O(h^4)$ と精度が大幅に向上する

この図では $y' = y$($\lambda = 1$、解が増大する方程式)を用いているため、前進・後退オイラー法の精度はほぼ同等に見える。後退オイラー法の真の利点は精度ではなく安定性にあり、硬い方程式で威力を発揮する(第5.3節参照)。

参考文献

  • J. C. Butcher, Numerical Methods for Ordinary Differential Equations, 3rd ed., Wiley, 2016.
  • E. Hairer, S. P. Nørsett, G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer, 1993.
  • U. M. Ascher, L. R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, SIAM, 1998.
  • Euler method — Wikipedia(英語)
  • オイラー法 — Wikipedia(日本語)

まとめ

この記事のポイント

  • 前進オイラー法:$y_{n+1} = y_n + hf(x_n, y_n)$(接線近似)
  • 局所打ち切り誤差 $O(h^2)$、大域誤差 $O(h)$(1次精度
  • 後退オイラー法(陰的):各ステップで方程式を解く必要があるがA安定
  • 安定性はテスト方程式 $y' = \lambda y$ で解析し、安定領域は $|R(z)| < 1$
  • 改良オイラー法(ホイン法):2次精度で前進オイラー法より高精度