# Why does my gravity simulation do this? [closed]

For a school project i created a simple 2D gravity sim in Matlab using the simplest possible method. There are 2 nested loops so that the total force and acceleration of every object can be calculated. I found that conservation of momentum in x and y works perfectly. I also use a damping so if two object are very close, their acceleration doesn't explode.

Then i place two objects of equal mass on coördinates (-100,0) and (100,0) with initial speeds 0. If i start the sim, the two objects will move towards eachother, pass though eachother and move away again. However, they move away much farther than where they originally started. Would that mean that kynetic energy + pot energy is violated, or is this normal?

In the real world this obviously never happens because objects can't pass through eachother.

It is not very important to solve the problem, the project is from a few months ago. I just wonder what went wrong.

Here is my code. The circle.m must be in the same folder as the gravity.m Some names of variables are in Dutch, but i translated the comments.

http://speedy.sh/4wVqm/Matlab-gravity2.rar

-

## closed as off-topic by tpg2114, kleingordon, Brandon Enright, Kyle Kanos, DavePhDJun 21 '14 at 13:42

• This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

The problem the @kleingordon describes is common. It is practically impossible to conserve energy with a discrete approximation to a continuous time evolution. There are techniques to reduce the effect by updating velocity and position in a certain sequence, but you probably won't want to bother with that. Make your time step smaller. This has problems, too (e.g. slows computation), but you should be able to improve your results. On the other hand, if you have to wait several periods to see the effect, you might not want to bother at all. – garyp Jun 20 '14 at 21:12
The adjective that describes algorithms that do not have this problem is "symplectic". Note that in general adjusting the time-step does not "fix" non-symplectic algorithms unless you only care about a fairly small number of time-steps. – dmckee Jun 20 '14 at 21:33
This question appears to be off-topic because it is about numerical methods and not physics. SciComp.SE would be a better place. – tpg2114 Jun 20 '14 at 21:53
You're almost certainly introducing error during the passing event. Try giving one of them a slight initial velocity so that they don't exactly collide. – Jerry Schirmer Jun 20 '14 at 22:00
Also, momentum conservation does not necessarily imply energy conservation. Probably the easiest way to debug this is to compute the energy at each timestep, and you can see whether your nonconservation is a gradual leak, or happens during the collision. Either way, there is almost certainly an error in your calculation, and what it is is unknowable without explicitly showing what you've done. – Jerry Schirmer Jun 20 '14 at 22:03

You are using Euler method to integrate this differential equation. This method has a smaller region of stability compared to other methods, such as Verlet integration (which is used in this web application) or one of the methods from the Runge-Kutta family.

Euler method uses the position at the beginning of the time step to compute the acceleration, which if the result would be 100% correct would mean that this value would be equal to the time average of the acceleration during this time step. However if there is some kind of exponentially like growth with will not be the case. This is also demonstrated on the Wikipedia page. In your simulation this is also the case since when the two object close in on each other and thus increasing the errors of the results. Analytically the velocities of the two objects should also go to infinity when $\vec{r}_i=\vec{r}_j$ if using energy conservation. Therefore it would be better to test this with (highly) elliptical orbits.

Since you where using Matlab you could even use ode45, which uses a 4th and 5th order Runge-Kutta method to also have variable time step sizes (so that it would take bigger times steps when the accelerations are small).

-
Unless I'm misremembering, standard Runge-Kutta, as ubiquitous as it is, isn't symplectic. So it's not as unconditionally unstable as forward-Euler, but it's still not very good for long-term N-body simulations. – Chris White Jun 21 '14 at 8:30
Also, Mathematica is more suitable for numeric integration in my opinion, because it supports multiple algorithms for integration, including symplectic ones, and is nowadays approximately as fast as matlab, provided that you know how to program it effectively. – auxsvr Jun 21 '14 at 9:19

As you point to yourself this problem happens because of the particles passing through each other. At the point when this happen $r\approx 0$ and the force $1/r^2$ blows up and numerically anything can happen here.

This is a generic problem in gravity simulations (so called N-body simulations). A numerical trick that is usually employed to avoid this is called force-softening. Instead of the inverse square law one instead use

$$\vec{F}_{12} = \frac{GMm(\vec{r}_1-\vec{r}_2)}{(|\vec{r}_1-\vec{r}_2|^2 + \epsilon^2)^{3/2}}$$

where $\epsilon$ is a small (but not too small) constant. This makes sure that the accelleration never blows up. The best value of $\epsilon$ depends on the system you want to solve and is usually found by experimenting.

I would also suggest you use the LeapFrog algorithm to solve the system (http://en.wikipedia.org/wiki/Leapfrog_integration). This is as simple to implement - and computationally just as cheap - as a standard forward Euler, but is much more stable for particles in orbit.

-
Yes! Leapfrog! +1. I plug this every chance I get. (Of course that fixes problems the OP hasn't encountered yet; force softening is more important here.) – Chris White Jun 21 '14 at 8:37

I found that conservation of momentum in x and y works perfectly.

Then i place two objects of equal mass on coördinates (-100,0) and (100,0) with initial speeds 0

So, I assume, the total momentum at each time step sums to zero to whatever precision you're working to.

But this tells you nothing about energy conservation. You can calculate the initial PE (the initial KE is zero) and the PE + KE at each step and, I'll bet, that the sum at each step does not equal the initial PE.

Momentum conservation only requires that the velocities are equal and opposite (in the case of equal masses) at all times since the initial momentum is zero.

$$\vec v_1[t] = - \vec v_2[t] = \vec v[t]$$

However, energy conservation, requires that

$$PE[0] = PE[t] + \frac{m}{2}(v^2_1[t] + v^2_2[t]) = PE[t] + mv^2[t]$$

So, the fact that the velocities are equal and opposite does not guarantee energy conservation; the velocities must be related to the initial PE and instantaneous PE as follows:

$$|\vec v_1[t]| = |\vec v_2[t]| = \sqrt{\frac{PE[0] - PE[t]}{m}}$$

-