Runge-Kutta法
をテンプレートにして作成
home
>
サイトマップ
開始行:
RIGHT:寄稿:東條遼平
* 4次近似の常微分方程式の数値解 [#zfcee7dd]
Heun法によって
随分と精度の高い常微分方程式の数値解を得ることができるよ...
もっと精度の良い4次のRunge-Kutta法を作成してみます。
Runge-Kutta法は5次以上のものもあるようなのですが、4次のも...
4次のRunge-Kutta法に手を加える形になるようです。
ところで自分は4次のRunge-Kutta法を導出することはできませ...
以下の図等で直感的に精度が良さそうだと思って頂ければ幸い...
#ref(runge1.png,nolink)
まずEuler法のときのような感じで
幅を半分で計算し傾きk1を算出します。幅dt/2の間、yの関数の...
あったと仮定しy(dt/2)の値を仮定します。その後y(dt/2)を初...
dt/2からdtまでの間のyの関数の傾きk2を算出します。Euler法...
これを解とし、どんどん先に進んでいくわけですが、Runge-Kut...
#ref(runge2.png,nolink)
0からdt/2までの範囲でyの関数の傾きが先程求まった傾きk2で...
もう一度y(dt/2)を求めてみます。その上で新たに求まったy(dt...
傾きk3を求めて再びy(dt)を算出します。
#ref(runge3.png,nolink)
今度は先程求まった傾きk3が0からdtの範囲の傾きであると仮定...
y(dt)を求めます。そしてその値より傾きk4を算出します。以上...
k1、k2、k3、k4を用いて、
#ref(img191.png,nolink)
のようにし、0からdtまでの傾きを算出し、最終的なy(dt)を決...
式で書くと次のようになります。
#ref(img192.png,nolink)
#ref(img193.png,nolink)
#ref(img194.png,nolink)
#ref(img195.png,nolink)
#ref(img196.png,nolink)
ここでk2とk3は共にdt/2の部分での傾きです。
そう考えると数値積分であるSimpson則の式、
#ref(img96.png,nolink)
と非常に似ている事が判ります。
ところで微分方程式を解くということは実は積分しているとい...
傾きを出して幅を掛け、高さを求めるというのも正に積分です。
とすると、数値積分のSimpson則の式と似ている事もうなずけ、
いきなりk2やk3を2倍したりしている事にもある程度納得できま...
しっかりと導出しようとすると大変な4次のRunge-Kutta法です...
Simpson則と絡める事で、Runge-Kutta法を使うには困らない程...
理解はできたと思います。ちなみにEuler法は積分の矩形公式、...
それぞれ対応します。
* ソースコード [#vb1b423b]
それではソースコードを見てみます。
h = (tn - t0);
if(div) h /= div;
t = t0;
for(i=0; i<div; i++){
for(j=0; j<num; j++){
k1[j] = (*f[j])(t, x);
temp[j] = x[j] + h*k1[j]/2;
}
for(j=0; j<num; j++){
k2[j] = (*f[j])(t+h/2, temp);
}
for(j=0; j<num; j++){
temp[j] = x[j] + h*k2[j]/2;
}
for(j=0; j<num; j++){
k3[j] = (*f[j])(t+h/2, temp);
}
for(j=0; j<num; j++){
temp[j] = x[j] + h*k3[j];
}
for(j=0; j<num; j++){
k4[j] = (*f[j])(t+h, temp);
x[j] += (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])*h/6;
}
t += h;
}
関数の引数は全く同じです。ただ、k1、k2、k3、k4と
確保しておく傾きの数が増えましたのでメモリの確保やループ...
連立の数を自由に指定することができるように、関数ポインタ...
初期値を配列で渡したりしているためにごちゃごちゃしてしま...
基本的に公式をそのまま計算しています。
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
終了行:
RIGHT:寄稿:東條遼平
* 4次近似の常微分方程式の数値解 [#zfcee7dd]
Heun法によって
随分と精度の高い常微分方程式の数値解を得ることができるよ...
もっと精度の良い4次のRunge-Kutta法を作成してみます。
Runge-Kutta法は5次以上のものもあるようなのですが、4次のも...
4次のRunge-Kutta法に手を加える形になるようです。
ところで自分は4次のRunge-Kutta法を導出することはできませ...
以下の図等で直感的に精度が良さそうだと思って頂ければ幸い...
#ref(runge1.png,nolink)
まずEuler法のときのような感じで
幅を半分で計算し傾きk1を算出します。幅dt/2の間、yの関数の...
あったと仮定しy(dt/2)の値を仮定します。その後y(dt/2)を初...
dt/2からdtまでの間のyの関数の傾きk2を算出します。Euler法...
これを解とし、どんどん先に進んでいくわけですが、Runge-Kut...
#ref(runge2.png,nolink)
0からdt/2までの範囲でyの関数の傾きが先程求まった傾きk2で...
もう一度y(dt/2)を求めてみます。その上で新たに求まったy(dt...
傾きk3を求めて再びy(dt)を算出します。
#ref(runge3.png,nolink)
今度は先程求まった傾きk3が0からdtの範囲の傾きであると仮定...
y(dt)を求めます。そしてその値より傾きk4を算出します。以上...
k1、k2、k3、k4を用いて、
#ref(img191.png,nolink)
のようにし、0からdtまでの傾きを算出し、最終的なy(dt)を決...
式で書くと次のようになります。
#ref(img192.png,nolink)
#ref(img193.png,nolink)
#ref(img194.png,nolink)
#ref(img195.png,nolink)
#ref(img196.png,nolink)
ここでk2とk3は共にdt/2の部分での傾きです。
そう考えると数値積分であるSimpson則の式、
#ref(img96.png,nolink)
と非常に似ている事が判ります。
ところで微分方程式を解くということは実は積分しているとい...
傾きを出して幅を掛け、高さを求めるというのも正に積分です。
とすると、数値積分のSimpson則の式と似ている事もうなずけ、
いきなりk2やk3を2倍したりしている事にもある程度納得できま...
しっかりと導出しようとすると大変な4次のRunge-Kutta法です...
Simpson則と絡める事で、Runge-Kutta法を使うには困らない程...
理解はできたと思います。ちなみにEuler法は積分の矩形公式、...
それぞれ対応します。
* ソースコード [#vb1b423b]
それではソースコードを見てみます。
h = (tn - t0);
if(div) h /= div;
t = t0;
for(i=0; i<div; i++){
for(j=0; j<num; j++){
k1[j] = (*f[j])(t, x);
temp[j] = x[j] + h*k1[j]/2;
}
for(j=0; j<num; j++){
k2[j] = (*f[j])(t+h/2, temp);
}
for(j=0; j<num; j++){
temp[j] = x[j] + h*k2[j]/2;
}
for(j=0; j<num; j++){
k3[j] = (*f[j])(t+h/2, temp);
}
for(j=0; j<num; j++){
temp[j] = x[j] + h*k3[j];
}
for(j=0; j<num; j++){
k4[j] = (*f[j])(t+h, temp);
x[j] += (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])*h/6;
}
t += h;
}
関数の引数は全く同じです。ただ、k1、k2、k3、k4と
確保しておく傾きの数が増えましたのでメモリの確保やループ...
連立の数を自由に指定することができるように、関数ポインタ...
初期値を配列で渡したりしているためにごちゃごちゃしてしま...
基本的に公式をそのまま計算しています。
- &ref(main.c);
- &ref(calculation.c);
- &ref(calculation.h);
ページ名:
home
>
Modified by
物理のかぎプロジェクト
PukiWiki 1.4.5_1
Copyright © 2001-2005
PukiWiki Developers Team
. License is
GPL
.
Based on "PukiWiki" 1.3 by
yu-ji
Powered by PHP 5.3.29HTML convert time to 0.002 sec.