SSブログ

4次精度のルンゲ=クッタ法で問題を解いてみた [ネコ騙し数学]

4次精度のルンゲ=クッタ法で問題を解いてみた


今日、11月24日に公開した数学の記事の第3問目は

  

になる。

この解は、問題3で求めたとおり

  

である。

(1)の連立常微分方程式は

  

の特殊なものであり、したがって、4次精度ののルンゲ=クッタ法で解くことができる。

連立常微分方程式を解く(4次精度の)ルンゲ=クッタ法のアルゴリズムは次のようなもの。

  

これを上から順次計算し、Step5まで計算しおえたら、xx+hだけ進め、Step1に戻って同じように計算すればよい。

ちなみに、上のアルゴリズムで使われている「=」は、数学の「等号」ではなく、右辺の計算結果を左辺に代入するという代入することをあらわしている。


4次精度のルンゲ=クッタ法はさすがに手計算でこれを行うのは大変なので、2次精度のルンゲ=クッタ法を使って、手計算で最初の値を求めてみることにする。


2次のルンゲ=クッタ法は次のようなアルゴリズム。

  


微分方程式(1)は

  

計算出発点のx=0のとき、y=1z=0

h=0.1とすると、

  


ちなみに、x=0.1のときの正しい値はy≒1.0050z≒0.10017だから、こちらも精度よく計算できていることがわかる。

そして、x=0.2のとき、上の計算で得られたy=1.0025z=0.1を使って同様の計算をする。

次に、4次精度のルンゲ=クッタ法を用いて、h=0.1としたときの計算結果を示す。




誤差は殆どなく、正確に計算できていることがわかると思う。


参考までに2次精度のルンゲ=クッタ法の計算結果を示す。




2次精度のルンゲ=クッタ法であってもこの程度の計算ならば、差は大きくない。

そこで、この両者の計算精度の差を明らかにするために、さらに複雑な次の連立常微分方程式を数値的に解くことにする

  

この微分方程式の解は、

  


4次精度の計算結果は次の通り。




厳密解とよく一致している。

しかし、2次精度のルンゲ=クッタ法では



と、xが大きくなるにつれて、厳密解との差が次第に大きくなってゆく。

h=0.1からh=0.01h1/10にしたとしても、精度はやはりh=0.1のときの4次精度のルンゲ=クッタ法に及ばない。



さらに1/10にしてh=0.001と、h=0.11/100したものが次のもの。


このとき初めて4次精度のルンゲ=クッタ法と2次精度のルンゲ=クッタ法は同程度の誤差になる。

精度が2次違うということはこういうことなんだケロよ。


4次精度のルンゲ=クッタのプログラムは以下の通り。

#include <stdio.h>
#include <math.h>

double f(double x, double y, double z) {
return x*y*z; // 変えていいのはここ
}

double g(double x, double y, double z) {
return x*y/z; //変えていいのはここ
}

// Runge_Kutta本体はいじらない方がいい
// いじらないでもいいように独立化させてある
// 独立化させているので、この部分は他のプログラムにそのまま流用可能!!
void Runge_Kutta(double x, double y, double z, double h ,double n) {
double dely1, dely2, dely3, dely4;
double delz1, delz2, delz3, delz4;
int i;

printf("%f %f %f\n",x,y,z); // 結果の出力(初期値)
for (i = 1; i <= n; i++) {
dely1 = h*f(x,y,z);
delz1 = h*g(x,y,z);
dely2 = h*f(x+h/2, y+dely1/2, z+delz1/2);
delz2 = h*g(x+h/2, y+dely1/2, z+delz1/2);
dely3 = h*f(x+h/2, y+dely2/2, z+delz2/2);
delz3 = h*g(x+h/2, y+dely2/2, z+delz2/2);
dely4 = h*f(x+h, y+dely3, z+delz3);
delz4 = h*g(x+h, y+dely3, z+delz3);
x=x+h;
y=y+(dely1+2*dely2+2*dely3+dely4)/6.;
z=z+(delz1+2*delz2+2*delz3+delz4)/6.;
printf("%f %f %f\n",x,y,z); // 結果の出力
}

}

main() {
int n;
double x, y, z;
double h;

h=0.1; // xの増分の設定 ココは変えてもよい
n=10; // 何回計算するかの設定 最終的にn×h進む
x=1; y=1./3.; z=1; // x,y,zの初期条件の設定
Runge_Kutta(x,y,z,h,n); // Runge-Kutta法を呼び出す

/* dy/dx=f(x,y,z), dz/dx=g(x,y,z)の初期値問題を
4次精度のルンゲ=クッタ法で解くプログラム */
// 注意 このプログラムで解く関数y,zはx=√7で不連続、ココでゼロ割発生!!
}


上のプログラムは、2つ目の微分方程式を解いている。
1つ目の微分方程式を解きたい場合、

double f(double x, double y, double z) {
return z; // 変えていいのはここ
}

double g(double x, double y, double z) {
return y; //変えていいのはここ
}

とし、mainの
x=1; y=1./3.; z=1; // x,y,zの初期条件の設定

 x=0; y=1; z=0;
にすればよい。
nice!(0)  コメント(0)  トラックバック(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

関数方程式関数方程式2 ブログトップ

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。