So-net無料ブログ作成
検索選択

同僚を鉄棒から宙づり、暴言も 消防士13人処分 福岡 朝日 [ネムネコの呟き]





水素水「痩せる」根拠なし 消費者庁が再発防止命令 産経 [ネムネコの呟き]





修正オイラー法 [ネコ騙し数学]

修正オイラー法


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

  

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

  


s-graph-01.png(2)の解は

  

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

  

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


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

  

で近似される。

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

そして、このとき、(3)は

  


s-graph-02.png一方、式(4)のこの定積分を台形公式で近似すると

  

(右図参照)

この図から明らかなように、(5)より(7)の方がより面積を精度よく近似しており、(3)を

  

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

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

  

が得られる

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


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

  

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


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

  


微分方程式

  

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

s-graph-03.png

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

  

である。

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


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

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

  

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

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

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

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

#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);}
}


スウェーデン、徴兵制復活へ ロシアに対抗、女性も対象 朝日 [ネムネコの呟き]





アサド政権軍、ISからパルミラ再奪還 シリア 朝日 [ネムネコの呟き]





保有の検討進むも「配備まで5年」 情報収集や米国の理解にも課題 産経 [ネムネコの呟き]






【三井美奈の国際情報ファイル】トランプ政権の中東和平「1国家」案に困ったのはイスラエル? 産経 [ネムネコの呟き]