Approximate change in position and velocity over timestep with varying acceleration due to gravitational interactions of many bodies

Question

I'm writing a program that models the gravitational interactions of a number of objects. After the number of bodies exceeds two, the motion can no longer be solved analytically, so the behavior of the chaotic system must be approximated using small timesteps.

At the moment, my code loops through each object hundreds of times a second, sums up the forces from the other objects using Newton's law of universal gravitation, and determines the resulting change in position and velocity over the timestep using $\Delta v = at$ and displacement $= v_it + \frac{1}{2}at^2$. (My actual code breaks this down into $x$ and $y$ components.)

However, the basic formulas I am using assume acceleration is constant throughout the timestep, when in reality it can can be quite a bit different at the end of a larger timestep. I implemented variable timesteps to make it a bit smoother, where the objects are iterated on more frequently when they are experiencing large acceleration and less frequently when the forces acting on them are small. However, there is still quite a bit of energy drift in order to make it run at a reasonable pace.

Is there a better way to more accurately predict change in velocity and position over each timestep?

Answer 1

You might try to decrease the number of loops so that the program is a bit faster though it becomes less accurate also .Also you may consider to program each body with respect to other and then integrate all the other possibilities(This might not work ,I am not confident).

Answer 2

there is still quite a bit of energy drift in order to make it run at a reasonable pace

This is a characteristic of poor numerical methods like Euler's method. You will need to use a better method. Unfortunately, it will require some re-coding, but these methods can run faster because you can use larger time steps. Euler's method is only an instructional tool.

What you want to look for is a symplectic integrator method. One of the more common and my personal recommendation would be the Verlet method. With these methods energy will be conserved. You may find numerical deviations, but they will not be of the kind that violates the conservation of energy.

See a related thread here.