境界値問題の前進積分による解法２

境界値問題の前進積分による解法２

 

次の微分方程式があるとする。

　　

この解は、

　　

 

２階常微分方程式（１）の境界値問題を４次のルンゲ・クッタ法で解く方法を思いついたので、紹介する。

微分方程式

　　

は、

　　

とおくと、

　　

となり、次の連立微分方程式に書き換えることが可能。

　　

そして、（２）は、x=x₀におけるy、pの初期値y₀、p₀が与えられれば、（２）は（４次の）ルンゲ・クッタ法を用いて前進的に解くことができる。

しかし、（１）の境界値問題では、pの初期値p₀が与えられていない。

そこで、y₀=0、p₀=1と推測し、n=10、h=0.1として、ルンゲ・クッタ法で解くと、次のような計算結果が得られる。

 

 

x=1におけるyの値はy(0)=1だから、として計算したyの誤差は

　　

次に、とし、ルンゲ・クッタ法を用いて解くと、次のようになる。

 

 

このとき、x=1におけるyの計算値は1.425591だから、誤差は

　　

である。

この計算結果を用いてx=0におけるpを次のように推測する。

　　

になる。

ところで、この微分方程式の厳密解は

　　

であり、これを用いてy'(0)を求めると、

　　

だから、何と、小数点５位まで正確な値を推測できる。

先に求めたを小数点７位で四捨五入した

　　

として計算した結果は次の通り。

 

 

x=1でのyの計算値は1になっており、この境界値問題を（４次の）ルンゲ・クッタ法を用いて見事解くことができた！！

 

３回も計算するのは面倒なので、推測機能を備えたプログラムを新たにC言語で作った。

　――最近、Fortranでずっとプログラムを書いていたから、Cでプログラムを書くことに、結構、手こずった(^^ゞ――

それを用いて解いた結果は次の通り。

 

 

計算に使用したプログラムは次の通り。

 

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

double f(double , double, double);
double g(double , double, double);
double exact(double);

void Runge_Kutta(double *x, double *y, double *z, double h, int n) {
    double dk[2][4];
    int i;
   
    for (i=0; i < n; i++) {
        dk[0][0] = h*f(x[i],y[i],z[i]);
        dk[1][0] = h*g(x[i],y[i],z[i]);
        dk[0][1] = h*f(x[i]+h/2, y[i]+dk[0][0]/2, z[i]+dk[1][0]/2.);
        dk[1][1] = h*g(x[i]+h/2, y[i]+dk[0][0]/2, z[i]+dk[1][0]/2.);
        dk[0][2] = h*f(x[i]+h/2, y[i]+dk[0][1]/2, z[i]+dk[1][1]/2.);
        dk[1][2] = h*g(x[i]+h/2, y[i]+dk[0][1]/2, z[i]+dk[1][1]/2.);
        dk[0][3] = h*f(x[i]+h, y[i]+dk[0][2], z[i]+dk[1][2]);
        dk[1][3] = h*g(x[i]+h, y[i]+dk[0][2], z[i]+dk[1][2]);
       
        x[i+1]=x[i]+h;
        y[i+1]=y[i]+(dk[0][0]+2*dk[0][1]+2.*dk[0][2]+dk[0][3])/6.;
        z[i+1]=z[i]+(dk[1][0]+2*dk[1][1]+2.*dk[1][2]+dk[1][3])/6.;
    }
}

main() {
    double x[N+1],y[N+1],p[N+1];
    double h=0.1;
    double p0=1., p1=2., w;
    double eps0, eps1, eps=0.000001;
    double yn;
    int i, cnt, n = N;


    x[0]=0.; y[0]=0.; p[0]=p0; yn=1.0; // p0はx=0におけるp=dy/dxの推測値、yn=y(1)はx=1におけるyの値
    Runge_Kutta(x, y, p, h, n);  // y0=1、p0=1としてルンゲ・クッタ法で計算
    eps0=y[n]-yn; // x=1におけるyの誤差を計算
    p[0]=p1; //p₀の推測値を更新
    Runge_Kutta(x, y, p, h, n);  // y0=1、p0=2としてルンゲ・クッタ法で計算
    eps1=y[n]-yn; // x=1におけるyの誤差を計算

    w=(eps1*p0-eps0*p1)/(eps1-eps0); // 補間法を用いてp₀を推測
    p0=p1; p1=w;
   
    cnt=1;
    while(cnt <= 20) {
        eps0=eps1;
        p[0]=p1;
        Runge_Kutta(x, y, p, h, n); // 推測値を用いて再計算
        eps1=y[n]-yn;  // 誤差を計算
        if (fabs(eps1)<eps) break; // 収束条件を満たしていたら計算終了
        w=(eps1*p0-eps0*p1)/(eps1-eps0); // 補間を用いてp₀を推測
        p0=p1; p1=w;
        cnt++;
    }

    printf("    x       y      厳密解\n");
    for(i=0;i<=n;i++) {
        printf("%f %f %f\n",x[i],y[i],exact(x[i]));
    }
}

double f(double x,double y, double z) {
    return z;
}

double g(double x, double y, double z) {
    return -2.*x*z - x*x*y;
}

double exact(double x) {
    double e;
    e=exp(1);
    return pow(e,3./2.)/(e*e-1)*(exp(x)-exp(-x))*exp(-x*x/2);
}

 

そして、「どうだ」と、胸を張るネムネコであった。

 

 

画像元：YouTubeの上の動画

 

さらに、自画自賛ソングを！！

 

 

なお、ループ（while）が存在するのは、非同次線形微分方程式だけではなく、非線形の微分方程式にも対応できるようにしたため。

線形微分方程式ならば、whileループは不要！！