ルンゲクッタ法(RK4)で斜方投射を解く

このQ&Aのポイント
  • C言語で作成した自由落下する物体の速度の変化を求めるプログラムを応用し、ルンゲクッタ法(RK4)を使用して斜方投射の位置の変化を求めるプログラムを作成したいです。
  • ルンゲクッタ法(RK4)を使い、斜方投射された物体のx座標とy座標の位置の変化を計算する方法を教えてください。
  • 斜方投射の問題を解くために、ルンゲクッタ法(RK4)を使用して物体の位置の変化を計算するプログラムをC言語で作成したいです。どのような式を立ててプログラムを組むべきか教えてください。
回答を見る
  • ベストアンサー

ルンゲクッタ法(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); } }

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

  • ベストアンサー
  • Tacosan
  • ベストアンサー率23% (3656/15482)
回答No.1

「プログラムを作成したい」ということだから, どのような式をたてればよいか自分で調べてみてはどうでしょうか. ちなみに式そのものは C や C++ とは無関係だね.

関連するQ&A

  • 微分方程式(ルンゲ=クッタ)がうまくいきません

    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); } }

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

    プログラムの実行ができません。 どこを直したらよいのか教えて下さい。 よろしくお願いします。 #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); }

  • Rubyのルンゲクッタ法がうまくいきません

    Rubyでルンゲクッタ法でy=exp(x**2)と比較しようという問題で、下のようなプログラムを組むとyの値がかなり大きくなってしまいます。どこが間違っているのでしょうか。教えてください。 #! ruby -Ks #ルンゲックッタ def df(a,b) z=2*a*b return z end n=10 x=Array.new(n) y=Array.new(n) x[0]=0.0 y[0]=1.0 h=1.0/n #任意で変更 fw=File.open("output2.txt","w") for i in 0..n k1=df(x[i],y[i])*h k2=df(x[i]+h/2,y[i]+k1/2)*h k3=df(x[i]+h/2,y[i]+k2/2)*h k4=df(x[i]+h,y[i]+k3) k=(k1+2*k2+2*k3+k4)/6 x[i+1]=x[i]+h y[i+1]=y[i]+k fw.print x[i]," ",y[i]," ","\n" end fw.close()

    • ベストアンサー
    • Ruby
  • 常微分方程式の数値解法: 陰解法のルンゲクッタ法の公式について

    陽的なRunge-Kutta法は, y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4) ただし,k1=hf(xi,yi), k2=hf(xi+h/2, yi+k1/2) k3=hf(xi+h/2, yi+k2/2) k4=hf(xi+h,yi+k3) と表すことは x=[x0,xn]の範囲でルンゲ・クッタ法により数値的にy(x)を解きたいのですが, 解こうとしている問題の初期条件がx=xnの時,y=0となっており,陰的(xnからx0に向かって)に解かなければならないのだろう,と考えています. 上記の公式でhの代わりに-hを入れて,プログラムを走らせても求めたい結果に大きな差異が生じてしまい困っています. そこで,陰的なRunge-kutta法の公式には,陽的な解き方と比較してどのような修正をすればよいか,教えてください. ちなみに解きたい微分方程式は, d2y/dx2 = x と仮定します. よろしくお願いいたします!!

  • 数値解析法

    この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; } }

  • 標的への斜方投射

    小球を打ち出し(l、h)にある的に当てる。θ度で投げ上げる I=V0cosθt h=V0sinθt-(1/2)gt^2 ここから初速度V0は? これを上の式からt=の形にして、下の式に代入 するのはわかるのですが、そっから全く計算があいません ちなみに答えはI√{g/2cosθ(Isinθ-hcosθ)} です、これになるための途中式を教えていただきたいです

  • ルンゲクッタ法で数値解析(C言語)

    コンピュータ(プログラミング)のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。 単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。 書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました(C言語を使って)。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。 ---- #include <stdio.h> double f1(double t,double x,double v); double f2(double t,double x,double v); int main() { double t,x,v,dt,tmax; double k1[2],k2[2],k3[2],k4[2]; x=1.0; //位置の初期値 v=0.0; //速度の初期値 dt=0.01; //刻み幅 tmax=500.0; //繰り返し最大回数 FILE *output; output=fopen("output.dat","w"); for(t=0;t<tmax;t+=dt) { k1[0]=dt*f1(t,x,v); k1[1]=dt*f2(t,x,v); k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k3[0]=dt*f1(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k3[1]=dt*f2(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k4[0]=dt*f1(t+dt,x+k3[0],v+k3[1]); k4[1]=dt*f2(t+dt,x+k3[0],v+k3[1]); x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0; v=v+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0; fprintf(output,"%f %f %f\n",t,x,v); } fclose(output); return 0; } double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); } ---- このソースは単振動のもので、減衰振動のときは、最後の方の ---- double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); ---- の部分が変わるのだと思うのですが、よく分かりません。 減衰振動は、mx"=-kx-rx'で表されるので、この式を変形(?)したものが、入るのかな、という予想程度にしか分かりません。 ソースをどのように変えれば減衰振動を解析できるのでしょうか。 どなたか詳しい方、教えてください。お願いします。

  • 微分方程式を解くアルゴリズムである ルンゲクッタ法

    微分方程式を解くアルゴリズムである ルンゲクッタ法の課題が出されました。 入力が周期1秒,最大値Eの矩形パルスということで 0.5秒周期でe(初期値10)を0⇒10⇒0⇒10⇒1とクロックさせたいのですが ファイルに出力してみるとずっと10のままで,0にクロックしていません・・・ 何処が原因か教えていただきたいです;; ほかにも問題があればご指摘していただきたいです>< ![イメージ説明](7309564d21bbda3453db49c0ec6147d9.jpeg) #include <stdio.h> #include <stdlib.h> #include <math.h> double clock(double E){ if(E==10) return 0; else return 10; } double dxdt(double x,int E){ double c = 0.001; double r = 50000; return (E-x)/(r*c); } // ルンゲクッタ法(初期条件x0, 区間[t0, tn]) double runge(double x0, double t0, double tn, int n) { FILE *output; FILE *output1;FILE *output2; output=fopen("output.txt","w"); output1=fopen("output1.txt","w"); output2=fopen("output2.txt","w"); double i; double x, t, h, d1, d2, d3, d4,cnt,e; x = x0; t = t0; e = 10; cnt=0; h = 0.01; // 漸化式を計算 for ( i=1; i <= n ; i++){ if(cnt==0.500000){e=clock(e); cnt=0;} t = t0 + i*h; d1 = dxdt(x,e); d2 = dxdt(x + d1*h*0.5,e); d3 = dxdt(x + d2*h*0.5,e); d4 = dxdt(x + d3*h,e); cnt= cnt + h; x += (d1 + 2 * d2 + 2 * d3 + d4)*(h/6.0); fprintf(output,"%f\n",t); fprintf(output1,"%f\n",x); fprintf(output2,"%f\n",e); } return x; } int main(void) { runge(0, 0, 1000, 50000); return 0; } ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

  • (C言語)プログラムが動きません。

    下のC言語のプログラムを書いたのですが、動きません。上手く作動していれば画面に「結果をrunge1.csvに書き込みました」と表示されてテキストファイルが作成されるはずですが、プログラムを実行しても何のメッセージも表示されずファイルも作成されません。ですが、bcpadのエラーメッセージの欄には何も表示されず、どこを直せばいいのか分かりません。アドバイスを下さい。お願いします。 /*ルンゲ・クッタ法による微分方程式y'+y*y=xの解*/ #include<stdio.h> #include<math.h> #define N 100 int main(void) { char *file="runge1.csv"; double x0=0.0,x1=5.0,y0=0.0,x,y,h,k1,k2,k3,k4; int i; FILE *fp; h=(x0-x1)/N; x=x0; y=y0; fp=fopen(file,"w"); fprintf(fp,"x,Runge-y\n"); fprintf(fp,"%3.2f,%3.5f\n",x,y); for(i=1;i<=N;i++){ x+=h; k1=(x-y*y)*h; k2=(x-y*y+k1/2.0)*h; k3=(x-y*y+k2/2.0)*h; k4=(x-y*y+k3)*h; y=(x-y*y)+(k1+2*k2+2*k3+k4)/6.0; fprintf(fp,"%3.2f,%3.5f\n",x,y); } fclose(fp); printf("結果を%sに書き込みました。\n",file); return(0); }

  • 斜方投射

    ボールを初速度V0で斜め上方に角度αで投げ上げたとき、ボールの到達距離と飛行時間がそれぞれSとTで、SとTよりV0およびαを求めたいのですが、下の式の(1)式と(3)式を使って代入を繰り返しているのですがうまく求まりません。どうか教えて下さい。よろしくお願いします。m(_ _)m 初速度V0の垂直成分をVyとすると、 Vy=V0sinα 飛行時間をt、垂直方向の変位をyとすると、 y=V0sinα・t-gt^2 変位y=0のとき、t=T(T≠0)となるので、 0=V0sinα・T-gT^2 ∴T=2V0sinα/g・・・(1) また、初速度V0の水平成分をVxとすると、 Vx=V0cosα ボールの到達距離Sを求めると、 S=T・Vx  =(2V0sinα/g)・V0cosα(∵(1)、(2)より)  =2V0^2/gsinαcosα =V0^2/gsin2α・・・(3) (1)、(3)式より・・・

専門家に質問してみよう