修正オイラー法

修正オイラー法

１階の常微分方程式の初期値問題の一般形は

　　

だが、より簡単な次のものを考える。

　　

（２）の解は

　　

g(x)>0の場合、(3)の右辺

　　

はg(x)とx軸と２直線で囲まれた部分の面積Sに等しい。

オイラー法の場合、(4)の積分は

　　

で近似される。

右図で正しい面積はピンク色の部分、(5)式で与えられる面積は斜線部で示されている。

そして、このとき、（３）は

　　

一方、式（４）のこの定積分を台形公式で近似すると

　　

（右図参照）

この図から明らかなように、（５）より（7）の方がより面積を精度よく近似しており、（３）を

　　

と近似したものの方が精度よく計算できるものと考えられる。

（１）と（２）はまったく同じ問題ではないけれど、yはxの関数であり、したがって、f(x,y)はxの関数であり、f(x,y)=g(x)と考えれば、（８）の結果を使って、

　　

が得られる

この近似計算法が修正オイラー法と呼ばれるものである。

⑨の右辺にが出てくるが、これは左辺のと、本来、同一のもので未知数。したがって、⑨を使ってを求めるためには、何らかの方法で右辺のを推測する必要がある。
そこで、

　　

と近似することによって、⑨の右辺を計算することができる。

修正オイラー法をもちいて微分方程式の初期値問題を近似的に計算するには次のようにすればよい。

　　

微分方程式

　　

を、分割幅h=0.1として、x=0からx=3まで、オイラー法と修正オイラー法を用いて数値的に解いたものが下図である。

赤の＋はオイラー法、青の＊が修正オイラー法を使って求めた数値解、そして、実線がこの微分方程式の解、すなわち、

　　

である。

図を見れば分かるように、修正オイラー法の方が厳密解とよく一致している。

計算する区間[a,b]の分割数nと計算精度の関係を右の図に示す。

精度の比較のために使用した微分方程式は

　　

a=0、b=2とし、計算する領域を[0,2]としている。

このグラフは、縦軸に厳密解との誤差、横軸に分割数をとったもので、縦軸、横軸ともに対数目盛りを使っている。

このグラフを見ると、オイラー法では分割数が１０倍になると誤差がおよそ1/１０に減少し、修正オイラー法では分割数が１０倍になると誤差が1/10²になっていることがわかると思う。
つまり、オイラー法の計算精度が分割数nの逆数（分割幅h）に比例しているのに対し、修正オイラー法は分割数の逆数（分割幅h）の２乗に逆比例するのであった。

参考までに、この計算に使用したプログラムを以下に示す。

#include <stdio.h>
#include <math.h>
#define N 100000

double f(double x, double y) { // dy/dx=f(x,y)=x-y
return x-y;
}

double exact(double x) { // dy/dx=x-yの厳密解
return 2.*exp(-x)+x-1;
}


// オイラー法
void Euler(double x, double *y, double h, int n) {
int i;

for (i=0; i <n; i++) {
y[i+1] = y[i]+f(x,y[i])*h;
x=x+h;
}
}

// 修正オイラー法
void IEuler(double x, double *y, double h, int n) {
double dy1, dy2;
int i;

for (i=0; i <n; i++) {
dy1 = h*f(x,y[i]);
dy2 = h*f(x+h,y[i]+dy1);
y[i+1] = y[i]+0.5*(dy1+dy2);
x=x+h;}
}


main() {
double y[N+1];
int n[] = {10, 50, 100, 500, 1000, 5000,10000};
int i;
double x0, y0,h;
double a,b;
double x;
double err1, err2;

a=0; b = 2.;
x0=a; y0=1.;
for (i=0; i < 7; i++) {
h=(b-a)/n[i];
x = x0; y[0] =y0;
Euler(x0,y,h,n[i]);
err1=fabs(y[n[i]]-exact(b));
IEuler(x0,y,h,n[i]);
err2=fabs(y[n[i]]-exact(b));
printf("%d %f %12.8f %12.8f1 \n",n[i],h,err1, err2);}
}

タグ：微分積分 数値解析