• 締切済み

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

微分方程式を解くアルゴリズムである ルンゲクッタ法の課題が出されました。 入力が周期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; } ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

みんなの回答

  • wormhole
  • ベストアンサー率28% (1619/5653)
回答No.1

0.01を浮動小数点で表現すると0.01の近似値であって、0.01そのものではありません。 そのため0.01を50回足しても0.5にはなりません。 例えば double d = 0.0; for (int i = 0; i < 50; i++) { d += 0.01; } printf("%s\n", (d == 0.5?"d equal 0.5":"d not equal 0.5")); を試してみてください。"d equal 0.5"とは出力されませんから。

関連するQ&A

  • 偏微分方程式、差分法

    Fitzhugh-Nagumo方程式 dU/dt = d/dx(dU/dx) + (a-U)(U-1)U-V (d/dx(dU/dx)はUのxに対する二回微分) dV/dt =e(bU-gV) a=0.1, b=0.5, g=1.0, e=0.01 初期条件 (U,V)=(1,0) if x=0 (U,V)=(0,0) if x>0 境界条件 dU/dx=0 at x=0 dV/dx=0 at x=0 を差分して陽解法で解くと添付の答えになるらしいのですが、自分で解いたところ添付の結果のようになりました。答えと一致しないため、プログラム上で何を間違えているのかご指摘頂けると助かります。 #include <stdio.h> #include <stdlib.h> //exsit()で必要 #include <math.h> int main(){ double a,b,g,e; double dt,dx; int x,t; double **u,**v; int i,j,k; a=0.1;b=0.5;g=1.0;e=0.01; x=100;t=5000.; dt=0.001;dx=0.1; FILE *fp; if((fp = fopen("ResultNagumo.txt","w"))==NULL){ printf("Can't open file\n"); exit(2); } u = (double**)malloc(sizeof(double*)*t); v = (double**)malloc(sizeof(double*)*t); for(i=0;i<t;i++){ u[i]=(double*)malloc(sizeof(double)*x); v[i]=(double*)malloc(sizeof(double)*x); } //初期値 u[0][0]=1.0; v[0][0]=0.0; for(i=1;i<x;i++){ u[0][i]=0.0; v[0][i]=0.0; }    //差分計算 for(i=0;i<t-1;i++){ for(j=1;j<x-1;j++){ u[i+1][j] = u[i][j] + dt*((u[i][j+1]-2*u[i][j]+u[i][j-1])/pow(dx,2)+(a-u[i][j])*(u[i][j]-1)*u[i][j]-v[i][j]); v[i+1][j] = v[i][j] + dt*(e*(b*u[i][j]-g*v[i][j])); } //境界条件 u[i+1][0]=u[i+1][1]; v[i+1][0]=v[i+1][1]; u[i+1][x-1]=u[i+1][x-2]; v[i+1][x-1]=v[i+1][x-2]; }   //結果の出力 for(i=0;i<t;i++){ if((i%100)==0){ fprintf(fp,"%d\n",i); for(j=0;j<x;j++){ fprintf(fp,"%2.4e,",u[i][j]); } } } for(i=0;i<t;i++){ free(u[i]); free(v[i]); } free(u);free(v); fclose(fp); getchar(); return 0; }

  • ルンゲクッタ法によるマクスウェルモデルの微分方程式

    大学の課題で dγ/dt=1/G(dσ/dt)+σ/η (但しt=時刻、Gは弾性要素の定数、ηは粘性要素の定数) に γ=1(t≧0) γ=0(t<0) のような階段状の歪みを加えた時 応力σ(t)の時間変化をルンゲクッタ法を用いて時間の刻み幅0.1 時刻20まで求めるプログラムを作成しろ というのがあります。作ってはみたものの ex.c: In function `main': ex.c:20: called object is not a function ex.c:22: called object is not a function ex.c:24: called object is not a function ex.c:26: called object is not a function と表示されて動きません。だれかどこをどう修正すれば良いか教えてください。 #include <stdio.h> double g(double t,double sigma,double z); double f(double t); int main(void){ double t,sigma,dt,count,tstart,z; double x1,x2,x3,x4; sigma=2.0; /*σの初期値を2とした*/ dt=0.1; /*刻み幅0.1*/ count=20.0; /*20まで表示*/ tstart=0.0; /*0から表示*/ for(t=tstart;t<count;t+=dt){ z=(f(t+dt)-f(t))/dt; /* z=dγ/dt */ x1=dt*g(t,sigma,z,); x2=dt*g(t+dt/2,sigma+x1,z); x3=dt*g(t+dt/2,sigma+x2,z); x4=dt*g(t+dt/2,sigma+x3,z); sigma=sigma+(x1+2.0*x2+2.0*x3+x4)/6.0; t=t+dt; printf("%f %f\n", t, sigma); /*tを表示  σを表示*/ } } double f(double t){ double ganma; if(t>=0)ganma=1; /*  γ=1(t≧0)  */ else ganma=0; /*  γ=0(t<0)  */ return(ganma); } double g(double t,double sigma,double z){ double ita,g,uhen; ita=10.0; g=2.0; /*与式を変形をし dσ/dt=(dγ/dt-σ/η)*G  右辺の関数がg()である*/ uhen=(z-sigma/ita)*g; /*η=10 G=2とした*/ return(uhen); } いろいろ調べてみたのですが、dy/dx=x+y などのような微分方程式の解き方は分かるのですが、上記のような微分方程式の解き方がわかりません。だから無理やりdγ/dtを求めて自分が解ける形にもっていったのですが、もっときれいな求め方があるならそれも教えてください。お願いします。

  • n自由度の振動のルンゲクッタ法

    n自由度の振動をルンゲクッタ法で計算したいと思い、プログラムの練習をしています。 現在2自由度まで計算するプログラムは組めているのですが、 nを入力し、そのnに対応した計算を行うルンゲクッタのプログラムが分かりません。 どなたか詳しい方に教えて頂きたいです。 以下に2自由度のプログラムを表示します。 #include <stdio.h> #include <math.h> #define NDIV 500; static double m1=1.0; static double m2=1.00; static double k1=1.0; static double k2=1.0; static double c1=0.30; static double c2=0.30; // dx/dt = kx double f(double x1,double v1){ return v1; } // dv/dt = kv double g(double x1,double v1,double x2, double v2){ return -((c1+c2)*v1-c2*v2+(k1+k2)*x1-k2*x2)/m1; } double j(double x2, double v2){ return v2; } double p(double x1, double v1, double x2, double v2){ return (c2*v1-c2*v2+k2*x1-k2*x2)/m2; } int main(void){ double x1; // 変位 double x2; double v1; // 速度 double v2; double c1; double c2; // 減衰 double t; // 時刻 double h; // =dt double hh; // =dt/2 double kx11,kv11,kx21,kv21; double kx12,kv12,kx22,kv22; double kx13,kv13,kx23,kv23; double kx14,kv14,kx24,kv24; int i,n=NDIV; char a[20],*pa; FILE *fp1,*fp2,*fp3,*fp4,*fp5; h = 0.1; hh = h/2; x1 = 1; v1 = 0; x2 = 0; v2 = 0; fp1 = fopen("ju-t.txt","w"); fp2 = fopen("ju-x1.txt","w"); fp3 = fopen("ju-v1.txt","w"); fp4 = fopen("ju-x2.txt","w"); fp5 = fopen("ju-v2.txt","w"); if(fp1 == NULL){ printf("can't open/n"); } else{ printf("open/n"); } printf(" # t x1 v1 x2 v2\n"); for(i=0;i<=n;i++){ t = h*i; printf("i,t,x1,v1,x2,v2"); fprintf(fp1, "%.4f\n",t); fprintf(fp2, "%.4f\n",x1); fprintf(fp3, "%.4f\n",v1); fprintf(fp4, "%.4f\n",x2); fprintf(fp5, "%.4f\n",v2); kx11 = f(x1,v1); kv11 = g(x1,v1,x2,v2); kx21 = j(x2,v2); kv21 = p(x1,v1,x2,v2); kx12 = f(x1+hh*kx11,v1+hh*kv11); kv12 = g(x1+hh*kx11,v1+hh*kv11,x2+hh*kx21,v2+hh*kv21); kx22 = j(x2+hh*kx21,v2+hh*kv21); kv22 = p(x1+hh*kx11,v1+hh*kv11,x2+hh*kx21,v2+hh*kv21); kx13 = f(x1+hh*kx12,v1+hh*kv12); kv13 = g(x1+hh*kx12,v1+hh*kv12,x2+hh*kx22,v2+hh*kv22); kx23 = j(x2+hh*kx22,v2+hh*kv22); kv23 = p(x1+hh*kx12,v1+hh*kv12,x2+hh*kx22,v2+hh*kv22); kx14 = f(x1+h*kx13,v1+h*kv13); kv14 = g(x1+h*kx13,v1+h*kv13,x2+h*kx23,v2+h*kv23); kx24 = j(x2+h*kx23,v2+h*kv23); kv24 = p(x1+h*kx13,v1+h*kv13,x2+h*kx23,v2+h*kv23); x1 = x1 + h*(kx11+2*kx12+2*kx13+kx14)/6; v1 = v1 + h*(kv11+2*kv12+2*kv13+kv14)/6; x2 = x2 + h*(kx21+2*kx22+2*kx23+kx24)/6; v2 = v2 + h*(kv21+2*kv22+2*kv23+kv24)/6; } fclose(fp1); fclose(fp2); fclose(fp3); fclose(fp4); fclose(fp5); return 0; }

  • 2自由度減衰振動のルンゲクッタ法

    2自由度の減衰振動をルンゲクッタ法で計算したいと思い、プログラムを組んでみました。 ただし、実行結果がどうも正しいものになっていないようです。 どなたか詳しい方、教えて下さらないでしょうか? 以下が組んでみたプログラムです。 #include <stdio.h> #include <math.h> #define NDIV 500; static double m1=1.0; static double m2=1.0; static double k1=1.0; static double k2=1.0; static double c1=0.00; static double c2=0.00; // dx/dt = kx double f(double x1,double v1){ return v1; } // dv/dt = kv double g(double x1,double v1,double x2, double v2){ return -((c1+c2)*v1-c2*v2+(k1+k2)*x1-k2*x2)/m1; } double j(double x2, double v2){ return v2; } double p(double x1, double v1, double x2, double v2){ return (c2*v1-c2*v2+k2*x1-k2*x2)/m2; } int main(void){ double x1; // 変位 double x2; double v1; // 速度 double v2; double c1; double c2; // 減衰 double t; // 時刻 double h; // =dt double hh; // =dt/2 double kx11,kv11,kx21,kv21; double kx12,kv12,kx22,kv22; double kx13,kv13,kx23,kv23; double kx14,kv14,kx24,kv24; int i,n=NDIV; char a[20],*pa; FILE *fp1,*fp2,*fp3,*fp4,*fp5; h = 0.1; hh = h/2; x1 = 1; v1 = 0; x2 = 0; v2 = 0; fp1 = fopen("ju-t.txt","w"); fp2 = fopen("ju-x1.txt","w"); fp3 = fopen("ju-v1.txt","w"); fp4 = fopen("ju-x2.txt","w"); fp5 = fopen("ju-v2.txt","w"); if(fp1 == NULL){ printf("can't open/n"); } else{ printf("open/n"); } printf(" # t x1 v1 x2 v2\n"); for(i=0;i<=n;i++){ t = h*i; printf("i,t,x1,v1,x2,v2"); fprintf(fp1, "%.4f\n",t); fprintf(fp2, "%.4f\n",x1); fprintf(fp3, "%.4f\n",v1); fprintf(fp4, "%.4f\n",x2); fprintf(fp5, "%.4f\n",v2); kx11 = f(x1,v1); kv11 = g(x1,v1,x2,v2); kx21 = j(x2,v2); kv21 = p(x1,v1,x2,v2); kx12 = f(x1+hh*kx11,v1+hh*kv11); kv12 = g(x1+hh*kx11,v1+hh*kv11,x2+hh*kx21,v2+hh*kv21); kx22 = j(x2+hh*kx21,v2+hh*kv21); kv22 = p(x1+hh*kx11,v1+hh*kv11,x2+hh*kx21,v2+hh*kv21); kx13 = f(x1+hh*kx12,v1+hh*kv12); kv13 = g(x1+hh*kx12,v1+hh*kv12,x2+hh*kx22,v2+hh*kv22); kx23 = j(x2+hh*kx22,v2+hh*kv22); kv23 = p(x1+hh*kx12,v1+hh*kv12,x2+hh*kx22,v2+hh*kv22); kx14 = f(x1+h*kx13,v1+h*kv13); kv14 = g(x1+h*kx13,v1+h*kv13,x2+h*kx23,v2+h*kv23); kx24 = j(x2+h*kx23,v2+h*kv23); kv24 = p(x1+h*kx13,v1+h*kv13,x2+h*kx23,v2+h*kv23); x1 = x1 + h*(kx11+2*kx12+2*kx13+kx14)/6; v1 = v1 + h*(kv11+2*kv12+2*kv13+kv14)/6; x2 = x2 + h*(kx21+2*kx22+2*kx23+kx24)/6; v2 = v2 + h*(kv21+2*kv22+2*kv23+kv24)/6; } fclose(fp1); fclose(fp2); fclose(fp3); fclose(fp4); fclose(fp5); return 0; }

  • 乱数について

    乱数の分布を見るために以下のようなプログラムを書きました。 #include <stdio.h> #include <stdlib.h> #include <math.h> int main() { int i,imax, S[RAND_MAX], r; double x,y; FILE *output1; output1=fopen("random2.data","w"); imax=100000; for(i=0;i<=imax;i++){ r = rand(); S[r] += 1; } for(i=0;i<=RAND_MAX;i++){ fprintf(output1,"%d %d \n",i,S[i]); } return 0; } するとコンパイルできて実行もできるのですが、なぜか乱数が30000を 超えるくらいのところでおかしな値になりました。 原因がわからないのでどなたか教えてください。

  • ルンゲクッタ法で数値解析(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'で表されるので、この式を変形(?)したものが、入るのかな、という予想程度にしか分かりません。 ソースをどのように変えれば減衰振動を解析できるのでしょうか。 どなたか詳しい方、教えてください。お願いします。

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

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

    プログラムの実行ができません。 どこを直したらよいのか教えて下さい。 よろしくお願いします。 #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
  • C言語 テキストの場合分け

    C言語プログラミングの質問です。 x1, y1, r1, distance1 x2, y2, r2, distance2 ・・・ と続くテキストを読み込み、distanceを10刻みで場合分けをするプログラムを作成しています。初心者で苦戦しているところですが、以下のサンプルの間違いをご指摘いただけませんでしょうか。 よろしくお願いいたします。 テキストの行数とdistanceの最大値は不明です。 ーーーーーーーー #include <stdio.h> #include <stdlib.h> #include <math.h> #include <string.h> #include <locale.h> int main(int argc, char *argv[]) { int a,b,c,distance; int i=0; int num[200]; FILE* inifile = fopen("file.txt", "r"); if(inifile==NULL) { printf("err__file.txt is nothing!!"); return 0; } FILE* outfile = fopen("outfile.txt", "w"); if(outfile==NULL) { printf("err__outfile.txt is nothing!!"); return 0; } while (fscanf(inifile,"%d\t %d\t %d\t %d\n", &a, &b, &c, &distance) != EOF) fclose(inifile); for(int i =1; i<101; i++) { if((i-1)*10 <= distance && distance < i*10) { num[i]+=1; } } for(int i =1; i<101; i++) { fprintf(outfile, "%d\t", num[i]); } fclose(inifile); fclose(outfile); printf("Normal END\n"); return 0; }