# Energy anomolies in many body simulation

I am trying to simulate the gravitational interaction between many bodies. I am using a direct PP force calculation and a 4th order symplectic integrator with a variable step size.

The energy of the simulation is calculated by this function:

real energy(body* system,int n){
int i,k,j;
real energy_var=0;
for(i=0;i<n;i++){
energy_var += 0.5*system[i].mass*pow(mag(system[i].vel),2);
for(j=i+1;j<n;j++){
if(i==j){continue;}
real position_vect[3]={0,0,0};

for(k=0;k<3;k++){
position_vect[k] = system[i].pos[k] - system[j].pos[k];}

real distance = mag(position_vect);

energy_var += -Gravitational_constant*system[i].mass*system[j].mass/distance;
}
}
return energy_var;}


In a test run with 100 bodies the energy( calculated every 100000's ) appears to be constant apart from one anomaly.

Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -4.24807e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41


If it was a numerical or coding error i.e. the timestep being too high then I would not expect the energy to return to the same constant value. Is there a known phenomenon that might cause such an error?

• If this was an inherent anomaly of the model, your code should still work for normal parameters. Start debugging your code with small examples that give results that are easy to interprete.
– Jan
Jul 31 '13 at 15:41
• I have run the simulation for 3 bodies with stable circular orbits and the energy stays constant to 5dp for a long time. They also seem to orbit correctly from looking at the data. Jul 31 '13 at 16:00
• Following on from @Jan, as you approach the anomaly take the energy and small time intervals and see if behaves in a physical way (e.g. smooth changes). How do you calculate the internal energy? Could this be error from numerical integration? Jul 31 '13 at 18:55
• There is something wrong in your example... you declare energy_var as real, but the energy is $\approx 10^{41}$ which would cause an overflow of a 32 bit IEEE floating point variable... Jul 31 '13 at 19:54
• real is a typdef from elsewhere in the program, in this run it is actually a long double. Jul 31 '13 at 19:58