1

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

ultraking x2
  • 43

3 Answers3

8

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).

fibonatic
  • 5,916
6

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.

Winther
  • 688
2

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}}$$

Alfred Centauri
  • 60,633