数値解析におけるEuler法とRunge-Kutta法(RK4)の精度について実験してみる（なにこのタイトルこわい）

今年になってEuler法やRunge-Kutta法という言葉を授業で頻繁に耳にする。
これらは微分方程式を数値計算で解く方法の一種である。

それぞれについて簡単に説明すると、

Euler法

Euler法は、微分方程式 $\frac{dy(t)}{dt}=f(t,y)$ と初期値 $y(t_0)=y_0$ が与えられたとき、
$y(t_0)=y_0$
$y(t+\Delta~t)=y(t)+f(t,y)\Delta~t$
を用いて漸化的に $y(t)$ を求める方法である。

Runge-Kutta法

Runge-Kutta法(RK4)は、Euler法に似た方法である。
同様に、微分方程式 $\frac{dy(t)}{dt}=f(t,y)$ と初期値 $y(t_0)=y_0$ が与えられたとき、
$y(t_0)=y_0$
$y(t+\Delta~t)=y(t)+\frac{\Delta~t}{6}(k_1+2k_2+2k_3+k_4)$
ただし、
$k_1=f(t,y)$
$k_2=f(t+\frac{\Delta~t}{2},y+\frac{\Delta~t}{2}k_1)$
$k_3=f(t+\frac{\Delta~t}{2},y+\frac{\Delta~t}{2}k_2)$
$k_4=f(t+\Delta~t,y+\Delta~tk_3)$
を用いて漸化的に $y(t)$ を求める方法である。

Euler法とRunge-Kutta法ではRunge-Kutta法の方が解の精度が高い。（見た目の複雑さからして一目瞭然である）
また、当然のことながら精度は $\Delta~t$ の値にも影響される。

今回は、 $\Delta~t$ の値を変えながらEuler法とRunge-Kutta法を用い、精度がどの程度変わるのかを実験してみようと思う。

実験方法

RC回路に一定の電圧 $v_i$ を入力した場合の出力電圧 $v_o(t)$ をシミュレートする。
$v_o(t)$ は以下の微分方程式で表せる。

$\frac{dv_o(t)}{dt}=\frac{1}{RC}\{v_i-v_o(t)\}$

簡単の為に $R=1,C=1,v_i=1$ とすると、

$\frac{dv_o(t)}{dt}=1-v_o(t)$ 　… (1)

この微分方程式の解は

$v_o(t)=1-e^t$ 　… (2)

である。

(1)式にEuler法やRunge-Kutta法を適用し、v_o(t)のグラフを描く。
それを(2)式のグラフと比較することで精度を視覚的に確かめることができる。

実験結果

まずはEuler法から。

赤線が(2)式から直接求めたグラフで、緑線がEuler法により求めたグラフである。

Euler法 dt=2.0（何の近似にもなっていない）

Euler法 dt=1.0

Euler法 dt=0.5

Euler法 dt=0.2

Euler法 dt=0.1

Euler法 dt=0.05

Euler法 dt=0.02

Euler法 dt=0.01

Euler法ではdt=0.01,0.02辺りでグラフがほぼ一致することが確認できる。

続いてRunge-Kutta法。

赤線が(2)式から直接求めたグラフで、緑線がRunge-Kutta法により求めたグラフである。

Runge-Kutta法 dt=2.0

Runge-Kutta法 dt=1.0

Runge-Kutta法 dt=0.5

Runge-Kutta法 dt=0.2

Runge-Kutta法 dt=0.1

Euler法とは比べ物にならない驚きの精度である。
dt=1.0の状態ですら既にグラフがほぼ一致している。
（間隔が広く、線がガクガクなので見にくいが、頂点はほぼ(2)式のグラフ上に位置している）

まとめ

Runge-Kutta法の効率の良さが良く分かる結果であった。
最後に本実験に使用したC++のプログラムを書いて本記事を終わる。

#include<iostream>
#include<fstream>
#include<cmath>
#include<iomanip>

using namespace std;

const double R(1.0);
const double C(1.0);
const double vi(1.0);
const double dt(1.0);

double f(double t,double vo){
	return (vi-vo)/R/C;
}

double dvo_euler(double t, double vo){
	return f(t,vo)*dt;
}

double dvo_rk4(double t, double vo){
	double k1,k2,k3,k4;
	k1=f(t,vo);
	k2=f(t+dt/2.0,vo+dt/2.0*k1);
	k3=f(t+dt/2.0,vo+dt/2.0*k2);
	k4=f(t+dt,vo+dt*k3);
	return dt/6.0*(k1+k2*2.0+k3*2.0+k4);
}

int main(){
double vo=0.0,t=0.0;
	ofstream ofs("sim.dat");
	ofs << setprecision(7) << t << " " << setprecision(16) << vo << endl;
	while(1){
		vo+=dvo_euler(t,vo);
		//vo+=dvo_rk4(t,vo);
		t+=dt;
		if(t>20.0)break;
		ofs << setprecision(7) << t << " " << setprecision(16) << vo << endl;
	}
	return 0;
}