ルンゲクッタ法でRC回路の出力関数を離散的にファイルに出力したい

Question

微分方程式を解くアルゴリズムである
ルンゲクッタ法の課題が出されました。
入力が周期１秒,最大値Eの矩形パルスということで
0.5秒周期でe(初期値10)を0⇒10⇒0⇒10⇒1とクロックさせたいのですが
ファイルに出力してみるとずっと10のままで,０にクロックしていません・・・

何処が原因か教えていただきたいです；；

ほかにも問題があればご指摘していただきたいです＞＜

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

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

Answer 1

計算機の浮動小数点数は内部が二進数表現になっており、またdoubleですと通常仮数が53bit程度ということで当然ながら無限の精度はありません。よって二進数で循環小数になるような数値は正確に表せない（常に誤差を含む）という点を知っておく必要があります。

本件でいえば0.01は循環小数になるので0.01 * 50は0.5に正確には一致しないと考えるべきです。一般に==で等値判定するなら整数を用いるべきです。0.01を50回加算したら0.5になると考えるのではなく1を50回加算して50かどうか判定するように頭を切り替えないといけないんですね。

本件は計算機プログラミングで「数学の世界と矛盾するためにわかっていないと大変混乱する」代表的な例の一つですが、他にもありますよね・・・

int i = 2147483647;
while (i > 0) {
  i++;
}
// なぜか無限ループにならない・・・