1. はじめに

N次方程式とは、未知数の最高次数が N である方程式である。

例として、2次方程式を考える。

$$
x^2 - 2x - 63 = 0
$$

これは次のように因数分解できる。

\begin{align}
x^2 - 2x - 63 &= (x-9)(x+7)
\end{align}

したがって解は

$$
x = -7, 9
$$

となる。

このように通常の方程式では、解は数値として得られる。

一方で、微分方程式では解は数値ではなく、関数となる。

例として次の微分方程式を考える。

$$
\frac{dy}{dx} = y
$$

この方程式の解は、

$$
y = Ce^x
$$

のような関数となる。

今回は微分方程式のうち、1つの変数だけで微分を行う常微分方程式を取り上げる。また、常微分方程式の解法として知られる次の2種類の解法を取り上げ、比較を Python で行う。

Euler 法（オイラー法）
Runge-Kutta 法（ルンゲクッタ法）

本記事で使用したコードは、次の GitHub リポジトリに掲載している。

2. 常微分方程式の解き方

ここでは、常微分方程式の基本的な解き方について説明する。常微分方程式には様々な解法が存在するが、ここでは最も基本的な例を扱う。

例として、はじめにで取り上げた次の問題を扱う。

$$
\frac{dy}{dx} = y
$$

変数分離を行うと、

$$
\frac{1}{y}dy = dx
$$

となる。

両辺を積分すると

$$
\int \frac{1}{y}dy = \int dx
$$

となる。

これを計算すると、

\begin{align}
\int \frac{1}{y}dy &= \int dx \\
\ln |y| &= x + C \quad \text{(C は積分定数)} \\
y &= e^{x + C} \\
  &= e^C e^x
\end{align}

ここで、$e^C$ は $x$ に依存しない定数である。そこで、$e^C$ を新しい定数 $C$ を使って書き直すと、

$$
y = Ce^x
$$

となる。

このように微分方程式の中には解析的に解くことができるものも存在する。

しかし実際には、解析的に解くことが難しい微分方程式も多い。
そのような場合に用いられるのが数値解法¹ である。

3. Euler 法

3.1. 概要

ここでは、常微分方程式の数値解法の一つである Euler 法（オイラー法）について説明する。

常微分方程式は一般に、

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

と表される。

ここで、$f(x, y)$ は「その点 $(x, y)$ における傾き」を表している。

例えば、現在の位置が $(x_n, y_n)$ であるとする。
このとき、微分方程式からその点での傾き

$$
f(x_n, y_n)
$$

が分かる。

Euler 法では、この傾きを使って現在の点 $(x_{n}, y_{n})$ から次の点 $(x_{n+1}, y_{n+1})$ を近似する。

この「少し先」に進む幅を $h$ とおく。

つまり、

$$
x_{n+1} = x_n + h
$$

である。

そして、傾きが $f(x_n, y_n)$ なので、

$$
\text{変化量} = h f(x_n, y_n)
$$

となる。

したがって、次の値 $y_{n+1}$ は

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

で求められる。

図3-1: Euler 法による1ステップの近似

これを初期値 $(x_0, y_0)$ から順に繰り返していく方法が Euler 法である。

3.2. Euler 法の実装

ここでは Euler 法を Python で実装していく。

例として次の問題を考える。

$$
\dfrac{dy}{dx} = y, \ y(0) = 1
$$

ここでは、

ステップ幅 $h$： 0.1
Euler 法の繰り返し回数 $N$： 50
$(x_0, y_0)$： $(0, 1)$

として、実装を行う。この実装で $(x_0,y_0)$ から $(x_N,y_N)$ までの点列を求める。

3.2.1. Python による実装

使用した Python コードは GitHub に掲載している。

3.2.2. Euler 法実行結果

Python コードの実行結果を次に示す。

ここで赤い曲線は微分方程式を解析的に解いて得られる真の解であり、

$$
y = e^x
$$

で表される。

途中までは常微分方程式の解き方で説明したように、解を求めていく。

ここで得た解

$$
y = C e^{x}
$$

に対して、$y(0) = 1$ より、

$$
y(0) = C e^{0} = 1
$$

これを解いて、

$$
C = 1
$$

よって、求める $y$ の値は、

$$
y = e^{x}
$$

である。

これと Euler 法によって得られた数値解を比較すると、$x$ が大きくなるにつれて両者の差が徐々に広がっていることが分かる。

しかし、全体としては解曲線の形をおおよそ近似できていることが確認できる。

Euler法は実装が簡単であるが，精度が低いという欠点がある。

この問題を改善するために考えられた方法が Runge-Kutta 法である。

4. Runge-Kutta 法

4.1. 概要

ここでは常微分方程式の数値解法の一つである Runge-Kutta 法（ルンゲクッタ法）について説明する。

Runge-Kutta 法は Euler 法を拡張した方法であり、複数の傾きを組み合わせることでより精度の高い近似を行うことができる。

Euler 法では現在の点 $(x_n, y_n)$ における傾き $f(x_n,y_n)$ のみを用いて次の点を計算した。

一方 Runge-Kutta 法では、区間内の複数の点で傾きを計算し、それらを組み合わせて次の値を求める。

この考え方により、Euler 法よりも高い精度で解を近似することが可能になる。

最も基本的な Runge-Kutta 法として、2次の Runge-Kutta 法（RK2）がある。

4.2. 2次 Runge-Kutta 法 (RK2)

RK2では次のように2つの傾きを計算する。

\begin{align}
k_1 &= f(x_n, y_n) \\
k_2 &= f(x_n + hk, \ y_n + hk k_1) 
\end{align}

これらを組み合わせて次の値を計算する。

$$
y_{n+1}
= y_n + h \left( (1-\frac{1}{2k})k_1 + \frac{1}{2k}k_2 \right)
$$

ここまでがRK2の一般形である。

4.3. RK2の係数による違い

この式にはパラメータ $k$ が含まれている。$k$ の値によって、いくつかの代表的な方法が得られる。

$k = \dfrac{1}{2}$ のとき

$$
1 - \dfrac{1}{2k} = 0
$$

よって、

$$
y_{n+1} = y_n + h k_2
$$

ここで、$k = \frac{1}{2}$ と $k_1 = f(x_n, y_n)$ を $k_2$ に代入して、

$$
k_2 = f(x_n + \dfrac{h}{2}, \ y_n + \dfrac{h}{2} f(x_n, y_n))
$$

以上より、

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

これは 修正オイラー法 (Modified Euler method) や 中点法 (Midpoint method) と呼ばれる。

$k=1$ のとき

$$
1 - \dfrac{1}{2} = \dfrac{1}{2}
$$

なので、$k=\frac{1}{2}$ の時と同様に代入して計算すると、

$$
y_{n+1} = y_n + \dfrac{h}{2}(k_1 + k_2)
$$

となる。これは ホイン法 (Heun's method) または、 改良オイラー法 (Improved Euler method) と呼ばれる。

4.4. 4次Runge-Kutta法 (RK4)

ここまでは2次の Runge-Kutta 法について解説をした。これをさらに改良したものとして、4次の Runge-Kutta 法がある。

一般的に Runge-Kutta 法というと、この4次の Runge-Kutta 法がある。

4次の Runge-Kutta 法は次の4つの傾きを計算する。

\begin{align}
k_1 &= f(x_n, y_n) \\
k_2 &= f(x_n + \frac{h}{2},\ y_n + \frac{h}{2}k_1) \\
k_3 &= f(x_n + \frac{h}{2},\ y_n + \frac{h}{2}k_2) \\
k_4 &= f(x_n + h,\ y_n + h k_3)
\end{align}

これらを組み合わせて、次の値を計算する。

$$
y_{n+1} = y_n + \dfrac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)
$$

4.5. Runge-Kutta 法の実装

ここでは Runge-Kutta 法を Python で実装していく。

例として Euler 法と同じ問題を考える。

$$
\dfrac{dy}{dx} = y, \ y(0) = 1
$$

Euler 法の時と同様に、

ステップ幅 $h$： 0.1
繰り返し回数 $N$： 50
$(x_0, y_0)$： $(0, 1)$

である。この実装で、Euler 法の時と同様に $(x_0, y_0)$ から $(x_n, y_n)$ までの点列を求める。

4.5.1. Python による実装

使用した Python コードは GitHub に掲載している。

4.5.2. Runge-Kutta 法実行結果

Python コードの実行結果を次に示す。

先ほどの Euler 法の実行結果と比較すると、RK4で実装した今回はより精度高く近似できていることがわかる。

5. Euler 法と Runge-Kutta 法の比較

ここまで2つの数値解析の手法である

Euler 法
Runge-Kutta 法

について説明をしてきた。

ここではそれぞれの手法を条件を変えながらシミュレーションを行うことで、それぞれの手法の特徴を調べる。

5.1. 実験概要

今回の実験では、空気抵抗を考慮した重力運動についてシミュレーションを行う。

自由落下運動では重力のみが働くと考えることもできるが、実際の物体の運動では空気抵抗の影響を受ける。そのため速度が大きくなるにつれて、抵抗力によって加速が抑えられる。

ここでは空気抵抗が速度に比例するモデルを考え、Euler 法と Runge-Kutta 法を用いて数値的に速度の変化を求める。

5.2. 運動方程式

質量 $m$ の物体が鉛直方向に落下するとき、物体には次の2つの力が働く。

重力
空気抵抗

重力は

$$
F_g = mg
$$

で表される。

一方、空気抵抗は速度に比例すると仮定すると

$$
F_d = -kv
$$

と書くことができる。

ここで

$k$：空気抵抗係数
$v$：速度

である。

ニュートンの運動方程式

$$
ma = F
$$

を用いると、物体に働く力のつり合いから

$$
ma = mg - kv
$$

となる。

図5-2: 自由落下する物体に働く力²

この式は、物体に働く重力と空気抵抗の関係を表している。

ここで加速度 $a$ は

$$
a = \frac{dv}{dt}
$$

であるため、

$$
m\frac{dv}{dt} = mg - kv
$$

となる。

両辺を $m$ で割ると、

$$
\frac{dv}{dt} = g - \frac{k}{m}v
$$

となる。

これは速度 $v$ に関する常微分方程式であり、時間とともに速度がどのように変化するかを表している。

5.3. 数値計算の設定

今回微分方程式

$$
\frac{dv}{dt} = g - \frac{k}{m}v
$$

を数値的に解くことで、時間とともに変化する速度を求める。

本実験では、次の条件でシミュレーションを行う。

変数	値	意味
`initial_y`	0.0	初期速度。今回は静止状態から落下を開始する
`g`	9.8	重力加速度
`k`	0.5	空気抵抗係数。値が大きいほど空気抵抗が強い
`m`	1.0	物体の質量
`dt`	各値	時間の刻み幅。数式では $h$ に対応する
`dt_list`	`[0.05, 0.5, 1.5]`	比較する刻み幅の一覧
`t_end`	10	シミュレーション終了時刻
`N`	`int(t_end / dt)`	`t_end` まで進めるための計算回数
`K`	0.5	RK2の係数（中点法に対応）

ここで dt は、1回の計算で時間をどれだけ進めるかを表す。数式で書いた Euler 法や Runge-Kutta 法のステップ幅 $h$ に対応する。

一方、t_end はシミュレーションをどこまでの時間で終了するかを表す。今回は t_end = 10 として、10秒後までの速度変化を計算する。

また N は計算回数であり、刻み幅 dt に応じて

$$
N = \frac{t_{\text{end}}}{dt}
$$

によって決定される。

これらの条件のもとで、Euler 法および Runge-Kutta 法による数値計算を行う。

5.3.1. Python による実装

使用した Python コードは GitHub に掲載している。

5.3.2 実験結果

Python コードの実行結果を次に示す。

図から、時間の刻み幅 dt が小さい場合には、いずれの手法でも精度高く、理論値に近似できていることがわかる。

しかし dt を大きくすると数値誤差が目立つようになり、手法ごとの差が現れる。

例えば dt = 0.5 の場合、Euler 法では理論解から大きくずれた結果となっている。一方、Runge-Kutta 法では比較的理論解に近い結果が得られている。

さらに dt = 1.5 の場合には、この差がより顕著となる。RK4 は依然として理論解に近い結果を示しているのに対し、Euler 法や RK2 では誤差が大きくなっている。

この結果から、Euler 法は刻み幅に強く依存しやすいのに対し、Runge-Kutta 法、特に RK4 は比較的安定して高い精度を保つことが確認できた。

ここまではそれぞれの手法を視覚化して、理論値と比較をした。

今度は数値解の精度を定量的に比較するため、最大誤差と二乗平均誤差（RMSE）を用いる。

最大誤差は、数値解と理論解の差の絶対値のうち、最も大きいものを表す。

E_{\max}
=
\max_{1 \le i \le N}
\left|
v_i - v_i^{exact}
\right|

ここで

$v_i$：数値解
$v_i^{exact}$：理論解

である。

最大誤差は、計算区間の中で最も誤差が大きくなった点を示す。
したがって、数値解がどの程度大きく理論解から外れる可能性があるかを評価することができる。

RMSE が全体的な平均的誤差を表すのに対し、最大誤差は最悪の場合の誤差を表す指標である。
この2つを併用することで、数値解の精度をより多面的に評価することができる。

$$
RMSE =
\sqrt{
\frac{1}{N}
\sum_{i=1}^{N}
(v_i - v_i^{exact})^2
}
$$

ここで $v_i$ は数値解、$v_i^{exact}$ は理論解である。

2つの誤差について、計算した結果を次に示す。

RMSE

dt	Euler	RK2	RK4
0.05	0.054980	0.000463	0.000000
0.5	0.579072	0.054159	0.000171
1.5	2.172982	0.752476	0.020999

最大誤差

dt	Euler	RK2	RK4
0.05	0.091081	0.000765	0.000000
0.5	1.008875	0.091132	0.000289
1.5	4.358384	1.158289	0.034389

表から、RMSEと最大誤差の両方で、刻み幅 $dt$ が小さいほど誤差が小さくなることが確認できる。

これは数値解法において、刻み幅を細かくするほど近似の精度が向上するためである。

また、同じ刻み幅で比較すると

$$
\text{Euler 法} < \text{RK2} < \text{RK4}
$$

の順に誤差が小さくなっていることが分かる。

特に $dt = 1.5$ のように刻み幅を大きくした場合、Euler 法では大きな誤差が生じているのに対し、RK4 法では比較的小さい誤差に抑えられている。

このことから、Runge-Kutta 法は Euler 法に比べて高い精度を持つ数値解法であることが確認できる。

この結果は、Euler 法が一次精度、RK2 法が二次精度、RK4 法が四次精度を持つという理論的な性質とも一致している。

6. まとめ

本記事では常微分方程式の数値解法として

Euler 法
Runge-Kutta 法

について説明し、Python を用いて数値シミュレーションを行った。

実験の結果、刻み幅を小さくするほど誤差が小さくなること、また同じ刻み幅で比較した場合には

Euler 法 < RK2 < RK4

の順に精度が高くなることが確認できた。

特に RK4 法は刻み幅が比較的大きい場合でも高い精度を保つことが分かった。

参考資料

[R1]. うさぎでもわかる微分方程式　Part00　微分方程式ってなに？

[R2]. 二乗平均平方根誤差【RMSE】Root Mean Squared Error

[R3]. 回帰分析の評価指標

常微分方程式の数値解法
https://ja.wikipedia.org/wiki/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%E5%BC%8F%E3%81%AE%E6%95%B0%E5%80%A4%E8%A7%A3%E6%B3%95 ↩
【物理】空気抵抗がある自由落下をラプラス変換を用いて解説
https://study-satellite.com/expertise/physics/free-fall/ ↩