• ベストアンサー

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

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

質問者が選んだベストアンサー

  • ベストアンサー
  • nitscape
  • ベストアンサー率30% (275/909)
回答No.3

＃１です 宣言 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) を満たすことがないという状態に陥っています。計算方法などを見直してみてはどうでしょうか？

gurizuri4649
質問者

お礼

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

gurizuri4649
質問者

補足

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

その他の回答 (3)

  • momoturbo
  • ベストアンサー率55% (49/88)
回答No.4

皆さんのご指摘があったものを直してますか？ 止まりましたけど・・ #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]); } }

gurizuri4649
質問者

お礼

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

  • neKo_deux
  • ベストアンサー率44% (5541/12319)
回答No.2

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

  • nitscape
  • ベストアンサー率30% (275/909)
回答No.1

ざっとプログラムを眺めただけですが。。。 #define MAX 2 ... double x,y[MAX]; ... if(y[3]>0.0){ というのはまずいのではないでしょうか？（y[3]は宣言されていないので） またループでもMAXを使っているようですのでそれが原因かもしれません。MAX周りを注意してデバッグしてみてはどうでしょうか？

gurizuri4649
質問者

補足

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

  1. カテゴリ
  2. [技術者向] コンピューター
  3. プログラミング・開発
  4. C・C++・C#

関連するQ&A

  • ルンゲクッタのプロクラム実行できません

    プログラムの実行ができません。 どこを直したらよいのか教えて下さい。 よろしくお願いします。 #include <stdio.h> #include <stdlib.h> #include <math.h> double *dvector(long i,long j); void free_dvector(double *a,long i); double func(double x,double y); double *rk4(double y0,double a,double b, int n,double(*f)()); int main(void) { double *y,h,a=0.0, b=1.0, y0=1.0; int i,n; printf("分割数を入力して下さい------->"); scanf("%d",&n); y=dvector(0,n); y=rk4(y0,a,b,n,func); h=(b-a)/n; for(i=0;i<=n;i++) { printf("x=%f \t y=%f \n",a+i*h,y[i]); } return 0; } double *rk4(double y0,double a, double b,int n,double (*f)()) { double k1,k2,k3,k4,h,*y,x; int i; y=dvector(0,n); h=(b-a)/n; y[0]=y0; x=a; for(i=0;i<n;i++) { k1=f(x,y[i]); k2=f(x+h/2.0,y[i]+h*k1/2.0); k3=f(x+h/2.0,y[i]+h*k2/2.0); k4=f(x+h,y[i]+h*k3); y[i+1]=y[i]+h/6.0*(k1+2.0*k2+2.0*k3+k4); x +=h; } return y; free_dvector(y,0); }

  • 数値解析法

    このHeun法のプログラムをRunge-Kutta法にするにはどうしたらいいですか？ #include <stdio.h> #include <math.h> double f1(double y) { return y; } double f2(double y) { return -4*y; } int main(){ double a=0; double b; int m=10; int n; double h; double x,y; int k; double e; double f; double k1,k2; printf("Heun法計算例：y=e^x, y=1/e^4x\n\n"); // y = e^x b = 1; for(n=100;n<=10000;n*=100){ h = (b-a)/n; printf("y' = y: h(=dx) = %.1e (y=e^x)\n",h); x = a; y = 1; for(k=0;k<=n;k++) { x = k*h; if(k%(n/m)==0) { f = exp(x); e = fabs(y-f); printf("x=%.2f, y=%f, e^x=%f er=%.0e\n",x,y,f,e); } // Heun's method k1 = h*f1(y); k2 = h*f1(y+k1); y += (k1+k2)/2; } }

  • Xcodeでのc言語のプログクムについて。

    Xcodeでのc言語のプログクムについて。 ルンゲークッタ法のプログラムを書いているのですがビルドエラーになってしまい、原因が分かりません。 　ビルド結果は、 "_func",referenced from: _main in main.o "_dvector",referenced from: _main in main.o _rk4 in main.o symbol(s) not found collect2: ld returned 1 exit status 　プログラムは、 #include <stdio.h> #include <stdlib.h> double *dvector(long i, long j); /* ベクトル量領域の確保*/ void free_dvector(double *a, long i); /*　領域の解放*/ double func(double x, double y); /* 関数の定義 */ /* ルンゲークッタ法 */ double *rk4( double y0, double a, double b, int n, double (*f)()); int main (void) { double *y, h, a=0.0, b=1.0, y0=1.0 ; int i, n; printf("分割式を入力してください-->"); scanf("%d",&n); y = dvector( 0, n ); /*領域の確保*/ y = rk4( y0, a, b, n, func ); /*ルンゲクッタ法*/ /*　結果の表示 */ h = (b-a)/n ; /* 刻み幅　*/ for (i = 0 ; i <= n ; i++) { printf("x=%f \t y=%f \n", a+i*h, y[i] ); } return 0; } /* ルンゲクッタ法　*/ double *rk4( double y0, double a, double b, int n, double (*f)() ) { double k1, k2, k3, k4, h, *y, x; int i; y = dvector( 0, n ); /* y[0,1,....,n]の確保　*/ h = (b-a)/n ; /* 初期値の設定　*/ y[0] = y0; x = a; /* ルンゲクッタ法　*/ for ( i = 0 ; i < n ; i++) { k1 = f(x,y[i]); k2 = f(x+h/2.0, y[i]+h*k1/2.0); k3 = f(x+h/2.0, y[i]+h*k2/2.0); k4 = f(x+h, y[i]+h*k3); y[i+1] = y[i] + h/6.0 * ( k1 + 2.0*k2 + 2.0*k3 + k4 ); x += y[i]; } return y; free_dvector(y, 0); /* 領域の解放　*/ }

  • モンテカルロ法を用いた実験

    モンテカルロ法を用いてy=x^2とx軸で囲まれた部分の面積の近似値を求める．これを利用して平均μ，標準偏差σ，を求めたところ下記の写真のようになった．点数nに対して標準偏差は何次であるかを，中心極限定理から求められる次数，実験結果から導かれる次数のそれぞれで求めなさい．次数とはどのように求められるのでしょうか？使用したプログラムを下に書いておきます． #include <stdio.h> #include <time.h> #include <stdlib.h> #include <math.h> double f(double x); int main(void) { int M = 1000; int i,j,n,c,k; int tmp[6] = { 20, 50, 100, 200, 500, 1000 }; double x,y,Y, sumY, sumY2, muY, muY2, sigmaY; srand((unsigned int)time(NULL)); for (k = 0; k < 6; k++){ sumY = 0; sumY2 = 0; for (j = 0; j < M; j++){ c = 0; n = tmp[k]; for (i = 0; i < n; i++){ x = (1.0* (double)rand()) / ((double)RAND_MAX + 1.0); y = (1.0* (double)rand()) / ((double)RAND_MAX + 1.0); if (y <= f(x)){ c++; } } Y = ((double)c) / (double)n; sumY = sumY + Y; sumY2 += Y*Y; } muY = sumY / (double)M; muY2 =sumY2 / (double)M; sigmaY = sqrt(muY2 -(muY*muY)); printf("%d %lf %lf\n", n, muY, sigmaY); } return 0; } double f(double x) { double r = x * x; return r; }

  • 座標をランダムに表示させてx座標順にソートするプログラムを考えています

    座標をランダムに表示させてx座標順にソートするプログラムを考えています とりあえず、以下の様に決まった数の座標でソートすることはできたのですが、ランダムにするとなるとどうすればいいのかわかりません。 ------------------------------------------------- #include <stdio.h> int makepoints(int * pn, double * x, double * y){ double xp,yp; int k; int i,j; int n; n = 7; *pn = n; xp = 1; for(k=0;k<n;k++) { xp = xp/2; yp = xp*xp; x[k] = xp; y[k] = yp; } printf("初期座標列：\n"); for(k=0;k<n;k++) { printf("%f_%f\n",x[k],y[k]); } for(j=1;j<n;j++) { for(i=0;i<j;i++) { if(x[i]>x[j]){ xp=x[i];x[i]=x[j];x[j]=xp; yp=y[i];y[i]=y[j];y[j]=yp; } } } printf("整列後の座標列：\n"); for(k=0;k<n;k++) { printf("%d %f %f\n",k ,x[k],y[k]); } return 0; } ------------------------------------------------- なんとなくrand関数を使えばいいのかな、というのはわかるのですが、プログラミングに弱く困っています。 この後のプログラミング教えてくださる方いればよろしくお願いします。

  • 二次元配列の引数渡し

    二次元配列を関数の引数として渡し、 その配列を戻り値として呼び出し元の変数に返したいです。 具体的にはルンゲクッタ法で微分する関数の引数を行列で扱いたいです。 for (i=0;i<N;i++){ for (j=0;j<N;j++){ y[i][j]=rand(); } } 　　K1=h*ff(t,y); K2=h*ff(t+h/2,y+K1/2); K3=h*ff(t+h/2,y+K2/2); K4=h*ff(t+h,y+K3); Ka=(K1+2*K2+2*K3+K4)/6; double ff(int t,double y[][]){ ここで return y[][]; のような形で二次元配列を戻り値として変数に返したいのです。 } Ｃ言語初心者なのでいまいちよくわかりません、 宜しくお願いいたします。

  • LU分解を利用した逆行列のプログラム（Java）

    LU分解を利用した逆行列のプログラムが作れません… というか、作ったのですが実行するとエラーが出てしまいます(´Д｀;) どこをどう直せばいいか、もしくはこのようにプログラムした方が効率がよい などのアドバイスどなたか下さい double a[][]={{2,5,４}, {2,3,-1}, {6,9,28}}; int N=a.length; double[][] s=new double[N][N]; for(int k=0; k<a[0].length-1; k++){ for(int i=k+1; i<N; i++){ s[i][k]=a[i][k]/a[k][k]; a[i][k]=s[i][k]; for(int j=k+1; j<N; j++){ a[i][j] -= s[i][k] * a[k][j]; } } } double[][] y=new double[N][N]; double[][] X=new double[N][N]; double[][] e=new double[N][N]; for(int i=0;i<N;i++){ for(int j=0;j<N;j++){ if(i==j){ e[i][j]=1; }else{ e[i][j]=0; } } } for(int i=0;i<N;i++){ y[1][i]=e[1][i]; for(int k=2;k<=N;k++){ for(int j=1;j<=N;j++){ y[k][i]=e[k][i]-s[k][j]*y[j][i]; } } X[N][i]=y[N][i]/a[N][N]; for(int k=N-1;k>=1;k--){ for(int j=k+1;j<=N;j++){ X[k][j]=(y[k][j]-s[k][j]*X[j][i])/a[k][k]; } } } for(int i=0;i<N;i++){ for(int j=0;j<N;j++){ System.out.printf(" %6.5f ", X[i][j] ); } System.out.println(""); }

  • ルンゲクッタ法（RK4）で斜方投射の問題を解く

    下のように自由落下する物体の速度の変化を求めるプログラムをC言語で作成しました。これを応用して斜方投射された物体のx、y座標の「位置」の変化を求めるプログラムをルンゲクッタ法（RK4）で作成したいです。x,yについて2つの式を立てないといけないと思うのですがどのような式を立ててプログラムを組めかよいかよく分かりません。どなたか式とできればプログラムを教えて下さい。 k：空気抵抗、m:物体の質量 g：重力加速度としています。 #include <stdio.h> double yd(double t,double y){ double k=0.1; double m=0.1; double g=9.8; double r=g-((k*y)/m); return r; } 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)); r=y+(h/6.0)*((k1+(2.0*k2)+(2.0*k3)+k4)); return r; } void main(){ double y=0.0; double h=0.001; double t=0.0; int i; for(i=0;i<100;i++){ t+=h; y=runge(t,y,h); printf("%f %f\n",t,y); } }

  • 関数宣言

    ３次元で領域を確保するプログラムをmalloc関数を用いて書きました。しかし、プログラムが長いので関数宣言をしなさいといわれたために、以下のプログラムを書きました。しかし、途中でつまづいてしまい、どのように関数を用いたり、関数を定義すれば良いのか混乱しています。初心者ですが、どうかお願いします。 /*ソース*/ #include<stdio.h> #include<stdlib.h> int main(){ double ***C; f3Malloc(C,.,.); //数値を代入（関数の使い方？） f3Free(C,.,.); } /*３次元配列（返し方？）*/ double*** f3Malloc(C,,){ int i,j,x,y,z; x = 2; y = 3; z = 4; C=(double***)malloc(sizeof(double**)*x*y*z); for(i=0;i<y;i++){ C[i]=(double**)malloc(sizeof(double*)*y*z); for(j=0;j<z;j++){ C[i][j]=(double*)malloc(sizeof(double)*z); } } } /*メモリの解放（返し方？）*/ void f3Free(C,.,.){ int i,j,x,y; x = 2; y = 3; for(i=0;i<x;i++){ for(j=0;j<y;j++){ free(C[i][j]); } free(C[i]); } free(C); }

  • 配列のエラーに関して

    java言語を用いて，Householder変換を用いた固有値の数値計算に挑戦してみました．しかし，次のようなエラーが発生し上手くいきません．どなたかこの問題を解決するためにお力をかしていただけないでしょうか． ----------エラー内容-------------------------------------------------------------------------------- Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: Index 0 out of bounds for length 0 at Out.Mhouse(House.java:90) at House.main(House.java:10) ---------------------------------------------------------------------------------------------------- //Householder変換 public class House{ public static void main(String[] args){ double[][] A = new double[3][3]; int n = A.length; Out out = new Out(); for(int i = 0;i < n;i++){ for(int j = 0;j < n;j++){ if(j < n-1){ System.out.print(out.Mhouse(A)[i][j] + " "); }else if (j == n-1) System.out.println(out.Mhouse(A)[i][j]); }; }; }; }; class Out{ double[][] outpro(double[] x){ int n; n = x.length; double[][] A = new double[n][n]; for(int i = 0;i < n;i++ ){ for(int j = 0;j < n;j++){ A[i][j] = x[i] * x[j]; } } return A; }; double[][] Msca(double a,double[][] A){ int n = A.length; for(int i = 0;i < n; i++){ for(int j = 0;j < n;j++){ A[i][j] = a * A[i][j]; } } return A; }; double selfpro(double[] x){ double a = 0; int n = x.length; for(int i = 0;i < n; i++){ a = a + x[i] * x[i]; }; return a; }; double[] minus(double[] x, double[] y){ int n = x.length; double[] z = new double[n]; for(int i = 0;i < n;i++){ z[i] = x[i] - y[i]; }; return z; }; double[][] house_1(double[] x){ int n = x.length; double[][] A = new double[n][n]; for(int i=0;i < n;i++){ for(int j = 0;j < n;j++){ if(i == j){ A[i][j] = 1 - Msca(2/selfpro(x),outpro(x))[i][j]; }else{ A[i][j] = - Msca(2/selfpro(x),outpro(x))[i][j]; }; }; }; return A; }; double[][] house_2(double[] x){ double[][] z = new double[1][1]; z[1][1] = 1 - 2; return z; }; double[][] Mhouse(double[][] A){ int n = A.length; double[][] H = new double[n][n]; for(int i = 0;i < n;i++){ double[] x = new double[n-i]; double[] y = new double[n-i]; double[][][] L = new double[i][n-i][n-i]; for(int j = 0;j < n-i;j++){ x[j] = A[i][i+j]; if(j == 0){ y[j] = 1; }else{ y[j] = 0; }; x[j] = y[j] - x[j]; }; if(i < n-1){ L[i] = house_1(x); for(int k = 0;k < n-i;k++){ for(int l = 0;l < n-i;l++){ H[i+k][i+l] = L[i][k][l]; }; }; }else if(i == n-1){ L[i] = house_2(x); for(int k = 0;k < n-i;k++){ for(int l = 0;l < n-i;l++){ H[i+k][i+l] = L[i][k][l]; }; }; }; }; double[][] B = new double[n][n]; for(int i = 0;i < n;i++){ for(int j = 0;j < n;j++){ for(int k = 0;k < n;k++){ B[i][j] = H[i][k] * A[k][j]; }; }; }; return A; }; };

    • ベストアンサー
    • Java

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

あなたにピッタリな商品が見つかる！ OKWAVEセレクト

質問する