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;
にすればよい。

関数方程式 [ネコ騙し数学]

関数方程式


問題1 abの任意の値に対してつねに、次の関係が成り立つように微分可能な関数f(x)を求めよ。

  

【解】

任意の定数aに対して

  

両辺をxで微分すると

  

f'(0)=Aとすると、

  

したがって、

  

また、

  

よって、B=0

したがって、

  

(解答終了)

 


問題2 すべての実数xに対して0でない関数f(x)があり、任意のx₁x₂に対して

  

を満足するものとする。

(1) f(0)=1であることを示せ。

(2) f'(0)=2であるとき、f'(x)の定義より、f'(x)=2f(x)であることを示せ。

(3) 任意の実数xに対してf(x)>0であることを示せ。

(4) (2)で与えられた微分方程式からf(x)を求めよ。

【解】

(1) x₁=x₂=0とすると

  

f(0)≠0だから、両辺をf(0)で割って、f(0)=1である。

(2)

  


(3)

  hk-eq-02.png

(4) y=f(x)とおくと、微分方程式f'(x)=2f(x)

  hk-eq-03.png

f(0)=1だから、C=1

よって、

  

(解答終了)


問題3 f(x)はすべての実数に対して定義され、正の値をとる連続関数で、

  

とおくとき、

  

である。

(1) F(x)=f(x)+g(x)G(x)=f(x)−g(x)とするとき、

  

であることを証明せよ。

(2) f(x)を求めよ。

【解】

(1)

  

xで微分すると、
  

F(x)=f(x)+g(x)xで微分すると

  

G(x)=f(x)−g(x)xで微分すると、

  

(2) y=F(x)とおくと、

  hk-eq-08.png

条件

  

x=0を代入すると、

  

すべての実数xに対してf(x)>0だから、f(0)=1

よって、

  

同様に微分方程式

  

を解くと

  

x=0のとき

  

したがって、

  

これをf(x)g(x)について解くと

  

(解答終了)

 


問題では、f(x)だけを求めよとあるけれど、g(x)も簡単に求められるので、求めた。

  

が成り立っていることが分かると思う。

実は、問題3で求めたf(x)g(x)は双曲線関数と呼ばれるもので

  

三角関数と類似した性質を持っている。

たとえば、

  


  

などなど。

そして、

  


タグ:微分積分

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