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

今日のお休みソング、「偽物語」から『Marshmallow Justice』 [今日のアニソン]

今日のお休みソングは、「偽物語」から『Marshmallow Justice』です。



さらに、「偽物語」からこの曲を。




nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:音楽

修正オイラー法 [数値解析]

修正オイラー法


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

今日のアニソン、「みつどもえ」から『あんたなんか!』 [今日のアニソン]

今日のアニソンは、「みつどもえ」から『あんたなんか!』です。



「天上天下、唯我独尊」を地で行くネムネコは、
他人(ひと)に対して対抗意識を持ったり、ライバル意識をもったりしないから、こういう感情が正直わからないにゃ。

ネムネコは、基本的に他者に興味・関心がないんだケロ(笑)。
他者に興味・関心がないから、他人(ひと)から褒められても、感謝の言葉をもらっても、まったく、嬉しいと思わないにゃ。



これ↑は、まだまだ、青いにゃ。
杉崎さんは、親から与えられもので自尊感情をもとうとしているにゃ。
この境地↓に達してナンボだにゃ。



ヒトを「ダンゴムシ扱い」でき、さらに「自虐」すら自然と「自尊」にるようになって、はじめて、ネムネコと同じ境地に立てると思うにゃ(^^)


nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:音楽