風の抵抗を考慮した放物運動

質問

風の抵抗を考慮した放物運動のプログラムを作ったのですが、無限ループが終了しません。どうすればよいでしょうか？以下、プログラムです。

#include<stdio.h>
#include<math.h>
#define MAX 2

double g=9.8; /*重力加速度*/
double k=0.0; /*抵抗係数*/

void main()
{

int rhs(double,double*,double*);

double x,y[MAX];
double h=0.01;
double x_old,y_old[MAX],f[MAX];
double k1[MAX],k2[MAX],k3[MAX],k4[MAX],k_work[MAX];
double z;
int i,j;

x=0;y[0]=28.3;y[1]=0.0;y[2]=28.3;y[3]=0.0;
printf("%f %f %f\n",x,y[1],y[3]);

for(i=1; ;i++){ /*ここは無限ループに*/
if(y[3]>0.0){
break;
} /*地面に落ちたら計算終了*/

x_old=x;
for(j=0;j<MAX;j++)y_old[j]=y[j];

x=i*h;

rhs(x_old,y_old,f);
for(j=0;j<MAX;j++){k1[j]=h*f[j];k_work[j]=y_old[j]+k1[j]/2.;}

rhs(x_old+h/2,k_work,f);
for(j=0;j<MAX;j++){k2[j]=h*f[j];k_work[j]=y_old[j]+k2[j]/2.;}

rhs(x_old+h/2,k_work,f);
for(j=0;j<MAX;j++){k3[j]=h*f[j];k_work[j]=y_old[j]+k3[j];}

rhs(x_old+h,k_work,f);
for(j=0;j<MAX;j++) k4[j]=h*f[j];

for(j=0;j<MAX;j++)
y[j]=y_old[j]+(k1[j]+2*k2[j]+2*k3[j]+k4[j])/6;

if(i%1==0)printf("%f %f %f\n",x,y[1],y[3]);
}
}

int rhs(double x,double y[],double f[])
{
f[0]=y[1];
f[1]=y[3];
f[2]=-k*pow(f[0]*f[0]+f[1]*f[1],0.5)*f[0];
f[3]=-k*pow(f[0]*f[0]+f[1]*f[1],0.5)*f[1]-g;
return 0;
}

回答No.3

nitscape · ベストアンサー

＃１です

宣言
int rhs(double,double*,double*);

実装
int rhs(double x,double y[],double f[])
{}

宣言と実装が違うというのもまずいと思います。

また実際に実行してみました。
if(i%1==0)printf("%f %f %f\n",x,y[1],y[3]);
によって終了条件であるy[3]が表示されていますが、始まりの値が0.0で、それからどんどん負の値になっています。そのため終了条件である
if(y[3]>0.0)
を満たすことがないという状態に陥っています。計算方法などを見直してみてはどうでしょうか？

お礼

おかげさまでちゃんと止まりました！
アドバイスありがとうございました！

補足

あ、間違っていました。
if(y[3]>0.0)はif(y[3]<0.0)です。
でも、終了しない。なんでだろう・・・
もう一度、根本から見直してみます。

回答No.4

momoturbo

皆さんのご指摘があったものを直してますか？
止まりましたけど・・

#include<stdio.h>
#include<math.h>
#define MAX 4

double g=9.8; /*重力加速度*/
double k=0.0; /*抵抗係数*/

int rhs(double x,double *y,double *f)
{
f[0]=y[1];
f[1]=y[3];
f[2]=-k*pow(f[0]*f[0]+f[1]*f[1],0.5)*f[0];
f[3]=-k*pow(f[0]*f[0]+f[1]*f[1],0.5)*f[1]-g;
return 0;
}

void main()
{

double x,y[MAX];
double h=0.01;
double x_old,y_old[MAX],f[MAX];
double k1[MAX],k2[MAX],k3[MAX],k4[MAX],k_work[MAX];
double z;
int i,j;

x=0;y[0]=28.3;y[1]=0.0;y[2]=28.3;y[3]=0.0;
printf("%f %f %f\n",x,y[1],y[3]);

for(i=1; ;i++){ /*ここは無限ループに*/
if(y[3]<0.0){
break;
} /*地面に落ちたら計算終了*/

x_old=x;
for(j=0;j<MAX;j++)y_old[j]=y[j];

x=i*h;

rhs(x_old,y_old,f);
for(j=0;j<MAX;j++){k1[j]=h*f[j];k_work[j]=y_old[j]+k1[j]/2.;}

rhs(x_old+h/2,k_work,f);
for(j=0;j<MAX;j++){k2[j]=h*f[j];k_work[j]=y_old[j]+k2[j]/2.;}

rhs(x_old+h/2,k_work,f);
for(j=0;j<MAX;j++){k3[j]=h*f[j];k_work[j]=y_old[j]+k3[j];}

rhs(x_old+h,k_work,f);
for(j=0;j<MAX;j++) k4[j]=h*f[j];

for(j=0;j<MAX;j++)
y[j]=y_old[j]+(k1[j]+2*k2[j]+2*k3[j]+k4[j])/6;

if(i%1==0)printf("%f %f %f\n",x,y[1],y[3]);
}
}

お礼

はい、ちゃんと止まりしました！
わざわざすみません！

回答No.2

neKo_deux

if(y[3]>0.0){

の手前にprintfでiとy[3]を表示する処理を入れると、手がかりになるのでは？

回答No.1

nitscape

ざっとプログラムを眺めただけですが。。。

#define MAX 2
...
double x,y[MAX];
...
if(y[3]>0.0){

というのはまずいのではないでしょうか？（y[3]は宣言されていないので）
またループでもMAXを使っているようですのでそれが原因かもしれません。MAX周りを注意してデバッグしてみてはどうでしょうか？

補足

MAX 2をMAX　4に修正しましたが、やはり改善はされません・・・