Why is f1 not equal to zero in this C code for a six spring mass simulation?

[Beeza]

I'll try to make this brief and more of a programming question than a physics question-- since I'm nearly positive my physics is correct.

I've written a simulation/code for a mass attached to six springs oscillating in 3-dimensions. The springs are each attached on the standard cartesian x-y-z axis. Take their anchor points to be + or - a on the x-axis. + or - b on the y axis, and + or - c on the z-axis. Now take a = b = c. Now, we should expect that given some initial displacement ONLY in the x-direction, there should be some oscillations only in the x-direction. This is true in my code, but I'm getting extremely small variations in the y and z directions.. oscillations on the order of 10^(-20).

One of the functions is as follows :
double result = 0.0,
f1 = 0.0,
f2 = 0.0,
f3 = 0.0,
f4 = 0.0;

f1 = b * ((y + b) / sqrt(x*x + y*y + z*z + 2.0*b*y + b*b) + (y - b) / sqrt(x*x + y*y + z*z - 2.0*b*y + b*b));
f2 = a * y * (1.0 / sqrt(x*x + y*y + z*z + 2.0*x*a + a*a) + 1.0 / sqrt(x*x + y*y + z*z - 2.0*x*a + a*a));
f3 = c * y * (1.0 / sqrt(x*x + y*y + z*z + 2.0*z*c + c*c) + 1.0 / sqrt(x*x + y*y + z*z -2.0*z*c + c*c));
f4 = 6.0 * y - (f1 + f2 + f3);
// printf ("%e %e %e %e\n", f1, f2, f3, f4);
result = - wo * wo * f4;
return (result);

Now, since given the initial conditions the first time through this function, as expected since y = 0, then f2 = 0 and f3 = 0. BUT, f1 does not equal zero! If we evaluate f1 for y = 0, then analytically equals zero, but the commented out print statement tells me something on the order of approximately 10^(-20). Is this just C being stupid and storing values wayyyy out in the decimals places or am I going crazy here? Have I just been staring at it too long and missed something stupid?

 

[Dr Transport]

You're compiler is using double precision, which machine are you using? When I was working with a CRAY, single precision was 16 decimals, double 32. If I remember correctly, on a PC single is 8 and double is 16 depending on the compiler. Now I would chalk up you $$10^{-20}$$ to round-off error and your compiler.

 

[Beeza]

Thanks Dr. Transport!

I'm just using a PC windows machine with JGrasp as my compiler.

The results are interesting. I expected the trajectory to be much more bizarre.

 

[dimensionless]

I'm getting weird errors too. I'm using GCC. When I add the compilation tag -O3, I get a different error.

 

[Beeza]

I'm getting weird errors too. I'm using GCC. When I add the compilation tag -O3, I get a different error.

huh? I don't understand.

 

[nmtim]

Floating point math is non-associative. More aggressive optimization levels permit the compiler to reorder operations, thus leading to different results.

You're counting on two numbers being identical except for sign and thus perfectly cancelling. That's not happening, indicating that perhaps the two numbers are arrived at with a different order of operations, or different operations are used. Without digging into the assembly, you probably won't know exactly what the compiler and the machine are doing.

On a Pentium-4 type PC (or Athlon, pretty much anything decent sold as a "PC" in the last five years) and with any reasonable compiler, a double means an IEEE-754 double. That's 64 bits. That doesn't map into an exact number of decimal digits, but typically I wouldn't worry too much about a relative error less than 10^{-14}.

 

[Beeza]

One would expect a plot of the total energy, call it E, to be constant right? Am I going crazy?

$$E = KE(\dot{x},\dot{y},\dot{z}) + U(x,y,z)$$

where U(x,y,z) is the potential energy function of the springs based on how far the mass is displaced from the origin (0,0,0). I'm getting a strange oscillating total energy. It seems as if not all of the potential energy is not being converted to kinetic energy.

 

[AlephZero]

If there are no external forces than energy should be conserved.

If your time-marching solution was not conserving energy, most likely the amplitude of the motion would increase with time, or die away to zero.

So I would guess your energy calculations are wrong.

Try making the model simpler, so you know what the "correct" answers are. For example set the springs in the Y and Z directions to zero and simulate

1. Oscillation along the X axis (simple harmonic motion)
2. Oscillation along the Y axis with springs in the X direction only. This is not simple harmonic motion, but you can check that the strain energy in the springs at the maximum displacement is correct.

 

[Beeza]

It's fixed. I re-derived all three Lagrangian equations -- What a pain in the neck! Here is the fixed energy plot.. but the title is wrong, it's 10 seconds.

Thanks Alephzero.

 

[dimensionless]

huh? I don't understand.

I tried compiling some C code at the Linux command line. Adding " -O3" to the command tells the compiler to optimize the output executable for faster execution.