数値解シミュレーション

質問

http://www.damp.tottori-u.ac.jp/~ooshida/edu/ode/euler2rk.html
に掲載されていることについてなのですが、そのサイトにも例として載っている物理法則の単純な　dx/dt=v , dv/dt=a=f/m　の近似を、サイトの下の方に載ってる４次Runge-Kutta法を用いてプログラムで表したいのです。最終的には>out.txtにしてエクセルに取り込んでグラフにしたいのですが、その為には多数の近似のプロット点が必要なのです。４次Runge-Kutta近似による多数のプロット点の値を出すようなプログラムを組んでもらえないでしょうか？知識がないもので宜しくお願いします。コンパイラはVBC++です。
以前似たような質問をして、http://www.maths.tcd.ie/~ryan/teaching/harm_osc_rk.c.html
のようなプログラムを組んで頂いたのですが、私の説明不足の為　dx/dt=v , dv/dt=a=f/m　の近似にはなっていないようなので、改めて質問させて頂きました。

回答No.6

tsukasa-12r · ベストアンサー

http://www.damp.tottori-u.ac.jp/~ooshida/edu/ode/euler2rk.html の説明では x[i] の i は時間に相当するものですが、
http://www.maths.tcd.ie/~ryan/teaching/harm_osc_rk.c.html では y[i] の i は
i = 0 : 位置 x
i = 1 : 速度 v
となっています。（私は最初、この点を大きく勘違いしていました。）
また、関数 f() は
i = 0 : 一次導関数 : dx/dt = v
i = 1 : 二次導関数 : d2x/dt2 = a
を表します。
あと、位置 y[0]、速度 y[1] の初期値の設定も忘れずに。

お礼

有難う御座いました。後はなんとか自力でやってみます。

回答No.5

tsukasa-12r

＞要するに、あの作ってもらったプログラムは間違っていないということでしょうか？
ええ、間違っていません。大丈夫です。ただ、関数 f() の中身が dt/dt = a (等加速度運動) の運動になっていないので、f() の中身だけ直さなければなりません。

＞実行してみてもdos画面に何も表示されないのですが・・・。
データは osc.dat というファイルに出力されていると思います。

補足

回答有難う御座います。なるほど、修正をすれば良いわけですね。

＞ただ、関数 f() の中身が dt/dt = a (等加速度運動) の運動になっていないので、f() の中身だけ直さなければなりません。

f()の中身が今はどうなっていて、どのように直せば宜しいのでしょうか？

回答No.4

2718281828

回答者同士の馴れ合いはナニですが、質問者さんにも有意義な情報であると確信するところがありますので、続けます。

> f() の引数 x が使われていないにもかかわらず、

四元四次の Runge-Kutta 法は完成されたアルゴリズムです。あれはイングランドの Trinity College の学生向けサンプルプログラムのようなので、一般のライブラリのソースをそのまま出したものだと思われます。サンプルコードで扱っている系は x が陽に入っていない自律系ですが、一般に i 次元の常微分方程式はx が陽に入って以下の形になりますから。

dy_1/dx = f_1(y_1, y_2, ... , y_i, x)
dy_2/dx = f_2(y_1, y_2, ... , y_i, x)
...
dy_i/dx = f_i(y_1, y_2, ... , y_i, x)

例えば、サンプルの系は線型減衰振子ですが、これに外力 Fsinωt が加わった強制振子になると x が陽の形で表れます。

実は、普通はサンプルコードのような書き方はしません（あれはあくまでもRunge-Kutta 法の学習用でしょう）。関数 f へポインタを runge4() に渡すようにして書きます。runge4() そのものはライブラリにしまっておくものですからね。

お礼

ちょっと自分にはお二方の内容の意味すら理解出来てないのですが、（苦笑）要するに、あの作ってもらったプログラムは間違っていないということでしょうか？ただ、あのプログラムをコンパイルしてみると、エラーは出ませんが、実行してみてもdos画面に何も表示されないのですが・・・。自分の目的は、多数のプロット点の値が欲しいのです。

回答No.3

tsukasa-12r

失礼しました。少しどころか思いっきり勘違いしていました。
次元についての for (i=0;i<N;i++) のループを、時間についての for (j=1; j*dist<=MAX ;j++) のループの意味と勘違いしていました。
2718281828 さん、ご指摘ありがとうございます。

そうすると、dx/dt=v , dv/dt=a になっていないというのは、double f(double x, double y[], int i)
の中身が、それに対応した形になっていないということが原因のようですね。
あと、まぁ、大した問題ではないのかもしれませんが、f() の引数 x が使われていないにもかかわらず、f(x+h, t1, i) のように x+h としているのが、意味がよくわかりません。

回答No.2

2718281828

通りすがりのものですが、気になるところがあるので、アドバイスします（回答ではありません）。失礼ながら、＃１さんは少し勘違いをなされているように思います。

> i が 0、1 以外のとき、返す値がありませんし、

この f は微分方程式を定義している関数です。今、二次元を考えているわけですから（N=2）、それが正常です。百歩譲れば、引数の i が N 以上ならエラーを返すように書くべきかもしれません。ですが、ここは実行時に最も多く呼び出される関数ですので、コードにあるように呼び出し側で制御するのが速度的に◎です（CPUタイムが百時間以内なら気にする必要もありませんが）。

> x[j+1] を求めて、その x[j+1] を使用して v1dt、v2dt、v3dt、v4dt を求めて、x[j+2] を求めて・・・というように求めていくはずなのに、

ご指摘の意図がわかりません。初期値 x[j] があって、そのx[j]を使用して v1dt…を求めて、x[j+1]を求めて…ということ（つまり引用されているアルゴリズム）と等価ではありませんか？あるいは、陰的に解けということでしょうか。陽で解いてまずい理由が特にあるわけでもなし、この考え方で問題ないと思います。

回答No.1

tsukasa-12r

質問っていうか、あなたが書かれているのはお願いですよね。（苦笑）

とりあえず、
http://www.maths.tcd.ie/~ryan/teaching/harm_osc_rk.c.html
見ましたけど、なんかおかしいような気がします。
double f(double x, double y[], int i)
{
　　if (i==0) return(y[1]); /* derivative of first equation */
　　if (i==1) return(-0.2*y[1]-y[0]); /* derivative of second equation */
}
って、i が 0、1 以外のとき、返す値がありませんし、

反復 {
　　v1dt = v(x[j] ) * dt
　　v2dt = v(x[j] + v1dt/2 ) * dt
　　v3dt = v(x[j] + v2dt/2 ) * dt
　　v4dt = v(x[j] + v3dt ) * dt
　　x[j+1] = x[j] + (v1dt + 2*v2dt + 2*v3dt + v4dt)/6
}
っていうのは、x[j+1] を求めて、その x[j+1] を使用して v1dt、v2dt、v3dt、v4dt を求めて、x[j+2] を求めて・・・というように求めていくはずなのに、最初に v1dt、v2dt、v3dt、v4dt をそれぞれ計算して最後に x を計算するようになってしまっていると思います。