微分方程式（ルンゲ＝クッタ）がうまくいきません

質問

C言語で物体の自由落下運動（空気抵抗を考慮する）を微分方程式を用いて解かなければならないのですが、どうしてもグラフの形が曲線になりません。どこが間違っているのか指摘お願いします。

空気抵抗を考慮しているので、落下速度はどんどんある一定の値に近づいていくはずなのですが、グラフは直線になってしまいます。いろいろと方法を試しましたがどれも外れました。
写真は赤が今の状態（間違っています）、黒が理想形です。
ここではルンゲ＝クッタ法を使用しています。
#include <stdio.h>
#include <math.h>

double yd(double t,double y){
double m = 50.0; 質量
double g = 9.8; 重力加速度
double k = 0.01; 空気抵抗の定数係数
double a= g - ((k *y)/m); 運動方程式
return(a);
}

double runge(double t,double y,double h){
double k1,k2,k3,k4,r;
double h2 = h / 2.0;
k1 = yd(t,y);
k2 = yd(t+h2,y+(h2*k1));
k3 = yd(t+h2,y+(h2*k2));
k4 = yd(t+h,y+(h*k3));
a = y + (h/6.0) * ((k1 + (2.0 * k2) + (2.0 * k3) + k4));
return(a);
}

void main(){
double y = 0.0;
double h = 0.001;
double t = 0.0;
int i;
for(i = 0; i < 1000; i++){
t += h;
y = runge(t,y,h);
printf("%f %f\n",t,y);
}
}

回答No.5

maru_yoshi_ · ベストアンサー

どこが間違っているかといえば、微分方程式の係数でしょうね。
double m = 50.0; 質量
double k = 0.01; 空気抵抗の定数係数
50Kgの質量に対して空気抵抗係数が0.01では終端速度が49Km/s（a = g - kv/m で a = 0となる v が終端速度）。
曲線になるには、かなりの時間がかかります。
係数をいろいろ変えてみれば、それなりのカーブになります。

> #4
> 空気抵抗力は，速度に比例するの？それとも速度の2乗に比例するの？それともそれ以外？
わかっていると思いますが「速度に比例する」モデルで計算しているようですね。

> #3
> double yd()の式は位置と加速度の関係を求めているようですが意味不明です。
違いますよね。速度と加速度でしょ。
まあ、仮引数名を y にしているのがなんだかなぁ。って感じですけど。
微分方程式を　y' = f(t, y); としているからかなぁ？

ところで、プログラムが「ルンゲクッタ法（RK4）で斜方投射の問題を解く」と全く同じで係数だけが異なっているように見えるのだが、質問者は同じ人？

回答No.4

f272

第1点
空気抵抗力は，速度に比例するの？それとも速度の2乗に比例するの？それともそれ以外？
第2点
どういう単位系の数値を使っているのかはっきりしないが，今の値で大丈夫なの？

とりあえず，今のままでもt=10000になるくらいまで計算させると，曲線が見えてくるよ。

回答No.3

ninoue

>double runge(double t,double y,double h){
,,,,
>a = y + (h/6.0) * ((k1 + (2.0 * k2) + (2.0 * k3) + k4));

>double yd(double t,double y)
....
>double a= g - ((k *y)/m); 運動方程式

上記の２個の式の対応が取れていません。
加速度、速度、位置の関係を再度検討して下さい。

double yd()の式は位置と加速度の関係を求めているようですが意味不明です。
速度と時間間隔(と空気抵抗、重力)をパラメータとして加速度を求め、結果として速度を返すのではないでしょうか。

runge()の式では速度と時間間隔をパラメータとして、その間の平均速度を返すようにするように出来ませんか。
平均速度や加速度には現在位置yの値は関係ないのでパラメータには不要と思われます。

平均速度と時間間隔を基に位置を求める計算は main()の中で行うようにすれば良いのではないでしょうか。

以上を元に再検討した結果を補足で記して下さい。
不明な点があれば追記下さい。

回答No.2

Tacosan

もともとの微分方程式.

回答No.1

asuncion

内容についてはよくわかりませんが、
コードを見て気になった点をいくつか。

＞double yd(double t,double y){

引数tを全く使っていませんが、正しいでしょうか。

＞double runge(double t,double y,double h){
＞a = y + (h/6.0) * ((k1 + (2.0 * k2) + (2.0 * k3) + k4));
＞return(a);

変数aの定義がありません。コンパイルできないはずです。