• 締切済み

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

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