5

I don't think the title is very accurate , sorry for that.

I simulate bodies in space using two timestep:

the TIMESTEP is the Δt wich I use to make the calculation and XTIME is the number of times I make the calculation

I want the each visible step to be 7 days , the Δt will be 86400 (seconds) and XTIME 7 , so it'll be calculate seven days.

I made that because if $Δt = 86400*7$ there is a lot of aberration (especially on moon orbit).

Since I'll use less than 20 planets each time , calculate 7 times Δt is not really a problem , but more could be.

And this is the problem , with $Δt = 86400$ I still have aberation when planet come too close to the massive mass (the sun), the result is the planet going away at the speed of light ! It's because the Δt is too large to calculate the "counter force" .

I heard of leapfrog but I must be stupid , I didn't succeed to adapt it to my code (the result was worst than anything, like the moon going away from earth) and i'm don't sure this will solve my problem.

I precise that I want something that look believable but not realistic , so if there is not "not-eating-cpu" method , I'll certainly going to cheat a little with those cases, but I prefer not to.

I ask this question here because i think it's more computing than physics since I already got the physics

edit: maybe a lead ? something like if (distanceToOldPosition > FIXED_VALUE) calculate more times

edit2: precision on calculations asked on comment: \begin{align} F = \frac{GM_{1} M_{2}}{r} \end{align}

the acceleration \begin{align} a = \frac{F}{M} \end{align}

then the linear vector \begin{align} v_{t+1} = v_{t} + aΔt \end{align}

and finally the new position

\begin{align} x_{t+1} = x_{t} + v_{t+1}Δt \end{align}

J. M.
  • 3,155
  • 28
  • 37
eephyne
  • 151
  • 3
  • welcome to SciComp. Could you be more specific in what algorithm you use to integrate the movement of the planets? – GertVdE Jan 31 '13 at 09:36
  • sorry for the delay , I added the formulas – eephyne Jan 31 '13 at 09:50
  • I am going to guess you are having problems due to the instability of Euler's method with your large time step. If possible, could you switch to backwards Euler to see if that improves your results? – Godric Seer Jan 31 '13 at 16:18
  • After reading docs on backward euler , I don't really get how to adapt it to my formulas – eephyne Jan 31 '13 at 16:50
  • I'm afraid you'll have to try a bit harder. Your question, to me, reads "please write the appropriate code for me". There is an incredible number of books on the forward and backward Euler methods; you should really give them a try, maybe also a fourth order Runge-Kutta method. – Wolfgang Bangerth Feb 01 '13 at 02:13
  • No No @WolfgangBangerth , I don't want the work to be done by another , I'm still working on it since yesterday – eephyne Feb 01 '13 at 06:00
  • By the way , did not the fourth order Runge Kutta is for system with fixed timesteps ? Since I have no "end" in mine this is not adapted, true? – eephyne Feb 01 '13 at 07:42
  • rk4 can be used for this case, just as easily as forward Euler, but I am guessing that it will have the same stability issues that Euler does. As Wolfgang Bangerth states, there is a large amount of literature on implicit methods that should be helpful to you. There also happens to be a recent question: http://scicomp.stackexchange.com/questions/5042/backward-euler-method here on SciComp that deals specifically with Backward Euler in more detail. – Godric Seer Feb 01 '13 at 15:28
  • But I think I missing something , in Euler or in RK4, the $h*f(x_{n}, y_{n})$ is used , $x_{n}$ is the total step since the beginning, right ? But I never use this value in my formula. So is there something wrong at the start with the formulas I use? – eephyne Feb 03 '13 at 07:20
  • Your equations are correct, your $f$ just happens to be independent of $x$. – Godric Seer Feb 08 '13 at 13:32
  • 1
    Related: http://scicomp.stackexchange.com/questions/5018 – Jed Brown Feb 08 '13 at 13:34

3 Answers3

11

To get something that looks realistic for planetary orbits, you shouldn't use the forward or backward Euler methods. These will cause your planets to spiral outward or inward. You should use a symplectic method.

You may also need to adjust the timestep to be smaller when two bodies are very close to each other.

Read Chambers (1999) A hybrid symplectic integrator that permits close encounters between massive bodies to get started with such methods.

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
2

In general I would suggest not writing your own ODE solver. This is a well solved problem, read about how popular ones work and pick one suitable to your problem.

If you let us know what language you are working in I am sure someone can suggest their favorite package. :)

Remember, these are never exact solutions, but some are less wrong that others.

meawoppl
  • 1,918
  • 11
  • 20
1

If don't need physics, you can just cheat and make your planets orbit in a perfect circumference or even better in a ellipse.

$$ pos_{planet} = circumference(angle(time),radius,center)$$

The moon and other satellites would be a bit different:

$$pos_{satellite} = pos_{planet} + circumference(angle(time),radius,center)$$

RSFalcon7
  • 159
  • 8