No idea why 0.83+0.01=0.839999

[gracinda710]

In the main program, I want to print out the value of t and the kinetic energy every 10*step only but it never prints. I found out the reason of this is because t never equals 10*step: t is initialized to 0 at the beginning and add 0.01(which is dt) in each iteration and it does the calculation correctly until it reaches t=0.83. In this iteration, 0.01 is added to t gives 0.839999 which mess up the calculation. Can somebody tell me how to fix it and why it happens?? Thanks!

[ code ]
#include<iostream>
#include<cstdlib>
#include"dislin5"
#include<cmath>

using namespace std;
using namespace dislin;

float X_max=10.0;
float Y_max=6.0;
float dt=0.01;
int Nparticles=48;

template<class T>
class Particle
{
T x;
T y;
T vx;
T vy;
T fx;
T fy;
float col[3];
static int ntot;

void wall(T& s, T& vs, T max)
{
float rad=0.1;
if (s>max-rad)
{
s=2.0*(max-rad)-s;
vs=-vs;
}
else if (s<rad)
{
s=2.0*rad-s;
vs=-vs;
}
}

public:
Particle(){
x=0.5*X_max;
y=0.5*Y_max;
vx=(1.0*rand()/RAND_MAX)-0.5;
vy=(1.0*rand()/RAND_MAX)-0.5;

for(int i=0; i<3; i++)
col[i]=0.8*rand()/RAND_MAX+0.2;

int nr=ntot++;
int Nx=static_cast<int>(sqrt(Nparticles/0.6));
int Ny=static_cast<int>(Nparticles/Nx);

if(Nx>1)
x*=(0.2+1.6*(nr%Nx)/(Nx-1));

if(Ny>1)
y*=(0.2+1.6*(nr/Nx)/(Ny-1));
}

static T KE;

void setforce(const T& fx_set, const T& fy_set)
{fx=fx_set; fy=fy_set;}

void addforce(const T& fx_add, const T& fy_add)
{fx+=fx_add; fy+=fy_add;}

void addvelocity(const T& vx_add, const T& vy_add)
{vx+=vx_add; vy+=vy_add;}

void addposition(const T& x_add, const T& y_add)
{x+=x_add; y+=y_add;}

void step()
{
addvelocity(fx*dt, fy*dt);
addposition(vx*dt, vy*dt);
wall(x, vx, X_max);
wall(y, vy, Y_max);
KE+=0.5*(vx*vx+vy*vy);
}

void show()
{DislinParticle(x,y,col);}

friend void calculate_forces(Particle* p, int N)
{
for(int i=0; i<N; i++)
p[i].setforce(0.0, -1.0);

T x_dist, y_dist, r2, r;

for(int i=0; i<N; i++)
for(int j=i+1; j<N; j++)
{
x_dist=p[i].x-p[j].x;
y_dist=p[i].y-p[j].y;
r2=(x_dist*x_dist)+(y_dist*y_dist);

if((r=sqrt(r2))<1.0)
{
T f_r=0.4/(r2*r);
p[i].addforce(f_r*x_dist, f_r*y_dist);
p[j].addforce(-f_r*x_dist, -f_r*y_dist);
}
}
}
};

template<class T> T Particle<T>::KE=0.0;
template<class T> int Particle<T>::ntot=0;

int main()
{
Particle<float> p[Nparticles];
float t=0.0;

DislinInit(10,"","Molecular Dynamics");

float step=0.1;

while(t<50)
{
t+=dt;
DislinFrame(X_max,Y_max);
Particle<float>::KE=0.0;
calculate_forces(p,Nparticles);

for(int i=0; i<Nparticles; i++)
{
p[i].step();
p[i].show();
}
endgrf();

if(t==10*step)
{
cout<< "The value of t is " << t << " and the kinetic energy is " << Particle<float>::KE << "\n";
step+=0.1;
}
}

DislinFinish();

return 0;
}

[ / code]

[Athar]

Floating point arithmetic is inaccurate.
See
http://en.wikipedia.org/wiki/Floating_point and
http://docs.sun.com/source/806-3568/ncg_goldberg.html

The bottom line is that you can't compare floating point values using ==. You need to test whether the value is within a small range of the number you want to compare against, i.e. if (abs(t-10*step)<epsilon)...

[shroomdoom]

Couldnt you use double?

[jimc]

Doubles are obviously more accurate, but they still suffer from the same issues; you can not trust comparison operations to function as they would for integers. This has to do with the fact that most floating point values are not exact, they are an approximation unless they happen to be something that comes out evenly in base 2 and fit within the precision of the value. Just like 1/3 is infinitely repeating in decimal (0.33333...), most real numbers do not come out evenly in base 2. So 1/10 comes out evenly in base 10 as 0.1, but not in base 2. Fractions that would be non-repeating in base 2 are those with powers of 2 in the denominator; 1/2, 1/4, 1/8, 1/16, etc.