0
$\begingroup$

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?

Improve this question
$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – Jan
    Jul 31, 2013 at 15:41
  • $\begingroup$ 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. $\endgroup$
    – user1750289
    Jul 31, 2013 at 16:00
  • $\begingroup$ 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? $\endgroup$
    – boyfarrell
    Jul 31, 2013 at 18:55
  • $\begingroup$ 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... $\endgroup$
    – Stefano M
    Jul 31, 2013 at 19:54
  • $\begingroup$ real is a typdef from elsewhere in the program, in this run it is actually a long double. $\endgroup$
    – user1750289
    Jul 31, 2013 at 19:58

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.