3

I am trying to create a custom solar system for a research project, to test different numerical methods and see their difference.

We are making our own solar system instead of trying to simulate earth/sun since that seemed harder.

For now we have a sun M1 and earth M2.

  • M1 = 3500
  • M2 = 500
  • distance from M1 to M2 = 200
  • We use a different G = 0.2

When we calculate the gravitional force it should be, F = 8.75 or -8.75 (F = GM1M2/r^2)

When we calculate the constant speed the earth should have to orbit around earth we find v = 2 (v = sqrt(G(M1+M2)/r))

To calculate the vec2 of the force we use the follow code

    static double GravitationalConstant = 0.2f;//6.674e-11;
    static public Vector2f GravitationalForce(SolarObject onObj, SolarObject fromObj)
    {
        Vector2f result = fromObj.OldPosition - onObj.OldPosition;

        float distance = result.Lenght();
        Console.WriteLine(distance);
        result = new Vector2f(result.X / distance, result.Y / distance); // Normalized, but I've the lenght al ready so

        result *= (float)((GravitationalConstant * onObj.Mass * fromObj.Mass) / (distance * distance));

        return result;
    }

And to update the position/velocity we use this

public void Update(float dt, List<SolarObject> objects)
    {
        foreach(SolarObject s in objects.Skip(1))
        {
            Vector2f f = Utility.GravitationalForce(s, objects[0]);

            Vector2f a = f / s.Mass;

            s.OldPosition  = s.NewPosition;
            s.NewPosition += dt * s.Velocity
            s.Velocity += dt * a;
        }
    }

The object does fly around the sun, but it's not an orbit at all. The distance between M1/M2 is not constant <-> the force is also not always equal to 8.75f. We know euler has an 'error', but this seems to big, since without even orbiting one circle 25% of the distance is al ready increased. So there has to be a bug somewhere.

Ruubb
  • 39
  • 2
  • Does it change if you are using doubles instead of floats? They are more precise. – Daniel Hilgarth Mar 28 '18 at 08:50
  • decimal has most precision – Ramunas Mar 28 '18 at 08:52
  • Also did you step through the code and verify the calculated values at each step? If there is an issue with the formulas you used, this should be the easiest way to detect it. – Daniel Hilgarth Mar 28 '18 at 08:56
  • Okay so I created an vector2d, everything is in doubles now. Still the same thing happens. At the start of the simluation the Force is 8.75 as expected, but things goes downhill really fast. – Ruubb Mar 28 '18 at 09:02
  • Ruubb use decimal not double. double d1 = 0.2 / (5.123456789); //0.039036144586872989 decimal d2 = 0.2m / (5.123456789m); //0.0390361445868729859214588957 And add breakpoints to code and check calculus. Double can be tricky, for example declared double d1 = 4.5 can become 4.49999999999 – Ramunas Mar 28 '18 at 09:59
  • I switched to doubles. At t = 1, I've a F = (0, -8.75), the position.X +2. And the velocity is making a circle too the sun. The frame after the F = (0, -8.75) and position seems correct. But is it not weird that x +=2, while there is no y movement. – Ruubb Mar 28 '18 at 10:30
  • Do you _know_ the error is too big or are you just assuming? The Euler method can get wrong results very quickly. – Turksarama Mar 28 '18 at 22:50

1 Answers1

5

Unfortunately, this behavior is intrinsic to Euler integration - you will always overshoot a curved path by some degree. This effect can be suppressed by using a smaller timestep, which works even without doubles:

  • dt = 0.01:

    enter image description here

  • dt = 0.001:

    enter image description here

  • dt = 0.0001:

    enter image description here

  • dt = 0.00001:

    enter image description here

As you can see, the Euler method's accuracy improves with decreasing timestep. The inner planets (smaller orbital radius = greater curvature = more overshoot) start to follow their projected orbits (green) more consistently. The overshoot becomes only just visible for the innermost planet at dt = 0.0001, and not at all for dt = 0.00001.

To improve on the Euler method without having to resort to ridiculously small timesteps, one could use e.g. Runge-Kutta integration (the 4-th order variant is popular).

Also, the orbital velocity should be v = sqrt(G*Msun/r)) rather than v = sqrt(G(M1+M2)/r)), although for the large star limit this should not cause too much of a problem.

(If you would like my test code please let me know - although it is very badly written, and the core functions are identical to yours.)

meowgoesthedog
  • 14,670
  • 4
  • 27
  • 40
  • 1
    Nice expalnation +1 for that btw the precision of integration on huge numbers can be boosted if needed see [Is it possible to make realistic n-body solar system simulation in matter of size and mass?](https://stackoverflow.com/a/28020934/2521214) – Spektre Apr 04 '18 at 14:41
  • @spektre while that is true, one can still normalize all relevant quantities to scale of unity in order to negate (to a limited extent) the effect of imprecision; what I tried to illustrate is that using Euler integration is impractical for any kind of simulation requiring stable results, regardless of precision - and thus the *algorithm*, not the data type, is the most important aspect. – meowgoesthedog Apr 04 '18 at 15:35