ねこ騙し数学、ページビュー累計１００万突破！！ [ひとこと言わねば]
ねこ騙し数学、ページビュー累計１００万突破！！
境界値問題の前進積分による解法２ [数値解析]
境界値問題の前進積分による解法２
次の微分方程式があるとする。
この解は、
２階常微分方程式（１）の境界値問題を４次のルンゲ・クッタ法で解く方法を思いついたので、紹介する。
微分方程式
は、
とおくと、
となり、次の連立微分方程式に書き換えることが可能。
そして、（２）は、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ループは不要！！
