Exploring the Reliability of Euler Method for Orbital Gravity Equations

[Timothy S]

Hello everyone,

I am curious as to if it is possible to use the Euler Method to solve the gravity differential equations? Would the approximations quickly diverge to inaccurate solutions, or would it stay relatively reliable?

Thanks

 

[HallsofIvy]

I presume you mean the basic "Euler method" for solving differential equations. (Given a differential equation dx/dt= f(x, t) with initial value $x(t_0)= x_0$, calculate $f(t_0,x_0)$ and construct the line through $(t_0, x_0)$ with slope $f(t_0,x_0)$. Follow that line for some short distance, $\Delta t$ to find $x_1= x_0+ f(t_0, x_0)\Delta t$ and $t_1= t_0+ \Delta t$.)

Such an approximation method can be made as precise as you want by taking very short $\Delta t$ but generally, no, the step required to make orbit calcultations would have to be so small as to make the method unreasonable, even on a computer. I would recommend a "Runge-Kutta" method of order 3 or higher.

 

[pbuk]

https://en.wikipedia.org/wiki/N-body_problem

ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/ExplSupplChap8.pdf explains how NASA's JPL does it using a variable-step, variable-order multistep method: Euler's method is the most simple method of this class being a single-step method of order one and is only useful for demonstrating how quickly errors propagate! An alternative class of methods is the Runge-Kutta methods: the fourth order method referred to as "RK4" is available in a number of libraries for C, Java etc.

 

[Timothy S]

Thank you, good point about the error propagation.

 

[Drakkith]

Hmm. Maybe the info in this thread can help me figure out why my orbital simulation I created Saturday likes to make the Earth spiral into the Sun and then go flying off into interstellar space.

 

[Timothy S]

Your simulation is probably not using a high enough order Runge - Kutta method, and is giving you increasingly inaccurate results as time goes on.

 

[Drakkith]

My simulation is using a couple of laws from my physics book I just learned, lol. I can't even tell you what method it is!

 

[Timothy S]

can you show what equations you are using?

 

[Drakkith]

I'm using the basic equations for motion along with Newton's 2nd law and law of gravitation:

X = X_o +Vt+1/2At²
V = V_o + At
F=MA
F = GM₁M₂/r²

I break up the equations into X and Y components with some trigonometry and then calculate each one using a time step. Then I add the resulting values to their respective variables and repeat. I didn't look up anything on this prior to creating the simulation. I just wanted to see if I could do it using what I've learned in my physics class.

 

[Timothy S]

Oh I see. You cannot do that with orbits as gravity is a differential equation and it is impossible to solve the orbital trajectory explicitly as a function of time.

You must first write the equation of motion of gravity in polar coordinates like so, \ddot{r}-r \dot{ \theta }^{2}=-\frac{GM}{ r^{2} }.
Because theta-dot however complicates things, we replace it with, \dot{ \theta } = \frac{h}{ r^{2} }. h is angular momentum. After this, to find r as a function of time, use a numerical method such as Runge-Kutta to find the trajectory as a function of time.

BTW, I am a high school sophomore too, good luck.

 

[Drakkith]

BTW, I am a high school sophomore too, good luck.

Gah! I've been outdone by a kid!
I'm 31 and in a basic college physics class, lol.

 

[Drakkith]

You must first write the equation of motion of gravity in polar coordinates like so, \ddot{r}-r \dot{ \theta }^{2}=-\frac{GM}{ r^{2} }.
Because theta-dot however complicates things, we replace it with, \dot{ \theta } = \frac{h}{ r^{2} }. h is angular momentum. After this, to find r as a function of time, use a numerical method such as Runge-Kutta to find the trajectory as a function of time.

Doesn't look like the forum liked your math.

 

[pbuk]

The OP's method is perfectly valid, and all you have done is replaced a numerical solution in cartesian coordinates with a numerical solution in polar coordinates.

However your assumption that angular momentum of a body is constant is only valid in general for a two body problem, where an analytical solution is available negating the need for numerical integration.

 

[Drakkith]

The OP's method is perfectly valid, and all you have done is replaced a numerical solution in cartesian coordinates with a numerical solution in polar coordinates.

That was the OP.

 

[pbuk]

That was the OP.

Oh yes, Drakkith, your method is perfectly valid!

 

[Drakkith]

Oh yes, Drakkith, your method is perfectly valid!

Dang! I thought I had a solution to why the Earth spiraled into the Sun in my simulation, but I guess not.

 

[tony873004]

You can compare Euler vs. RK4 here. http://orbitsimulator.com/BA/compare.html . It actually uses Euler-Cromer, the difference being that velocities are updated before positions are updated. It's javascript, so you can view the code and look at functions RK4() and function Euler(). Euler's method is only a few lines long. Runge-Kutta 4 is slower, but you can take larger time steps which makes up for the difference. But Euler's method still works pretty well.

I didn't look up anything on this prior to creating the simulation. I just wanted to see if I could do it using what I've learned in my physics class.

That's what I did about 12 years after graduating high school. I remembered my physics, but that capital letter G in all the gravitation formulas still confused me. I decided to give it a go anyway. My units were pixels per iteration. I was so happy when I saw a trajectory get bent by a mass. I spent all night zeroing in on the perfect value for a circular orbit. Months later, I wanted to use meters / second instead of pixels per iteration. So I introduced a "fudge factor" in my equations which made "Earth" take 365.24 days to go around "Sun". Later, I learned that what I came up with was Euler-Crommer method, and my fudge factor was G.

 

[pbuk]

Dang! I thought I had a solution to why the Earth spiraled into the Sun in my simulation, but I guess not.

One problem can be too SMALL a time step: you need δt to be large enough that δx/x >> ∈ where ∈ measures the precision of the arithmetical engine (or rather that |δx|/|x| >> ∈ where δx and x are vectors).

 

[Drakkith]

That's what I did about 12 years after graduating high school. I remembered my physics, but that capital letter G in all the gravitation formulas still confused me. I decided to give it a go anyway. My units were pixels per iteration. I was so happy when I saw a trajectory get bent by a mass. I spent all night zeroing in on the perfect value for a circular orbit. Months later, I wanted to use meters / second instead of pixels per iteration. So I introduced a "fudge factor" in my equations which made "Earth" take 365.24 days to go around "Sun". Later, I learned that what I came up with was Euler-Crommer method, and my fudge factor was G.

I used meters and seconds in all of the math and just divided the distances by about 2 million or something when I updated their position on screen.

 

[Drakkith]

One problem can be too SMALL a time step: you need δt to be large enough that δx/x >> ∈ where ∈ measures the precision of the arithmetical engine (or rather that |δx|/|x| >> ∈ where δx and x are vectors).

I'm not familiar with any of that notation. For my timestep I used 12 hours. Since my program runs 60 timesteps per second max, making it much smaller would make the orbit take forever. However, I think I can get around this 60 fps limit by using a 'repeat X number of times' in my program, so it will run each calculation 60X times. Haven't tried it yet though.

 

[pbuk]

Oh 12 hours is MUCH too large, you won't get anything stable with Euler's method until you are around a minute* and with double precision arithmetic you can get down as low as a second* and still have sensible rounding errors that won't become a problem for tens of thousands of steps. Having said that, working in polar coordinates is less subject to loss of precision (because in Cartesian coordinates δx/x is < ∈ when x is small, which it is twice each orbit, whereas in polar coordinates both δr/r and δθ/θ are bounded).

You should be able to calculate thousands of steps per second, it is displaying them that takes the time.

* these are guesses off the top of my head

 

[Drakkith]

You should be able to calculate thousands of steps per second, it is displaying them that takes the time.

Right now I'm limited by the software I'm using, which is video game development software (Construct 2). I think I can get the needed calculations per second, but I'll have to try some things.

 

[tony873004]

Oh 12 hours is MUCH too large, you won't get anything stable with Euler's method until you are around a minute*
* these are guesses off the top of my head

That's a very good guess. Here's another copy of the above simulation. In this version, The planets and Moon are set in their positions today, using JPL Horizons data. When you press play, they are propagated to August 2045, when a total solar eclipse passes through Redding CA. The time step automatically slows as you arrive at the eclipse. The simulation begins with the integrator set to RK4 and a time step of 2048 seconds. You will notice the eclipse happens on schedule. You can refresh and try again at different time steps. RK4 gives you pretty good results with time steps up to about 4 hours. Switch to Euler. You need time steps of about 2 minutes or less to get the eclipse to happen on schedule.

You should be able to calculate thousands of steps per second, it is displaying them that takes the time.

In the following sim, the number "Do Events", which is set to 500, means do 500 iterations for every 1 graphic update.
Simulation:
http://orbitsimulator.com/BA/compareEclipse.html

Oh 12 hours is MUCH too large, you won't get anything stable with Euler's method until you are around a minute* and with double precision arithmetic you can get down as low as a second* and still have sensible rounding errors that won't become a problem for tens of thousands of steps. Having said that, working in polar coordinates is less subject to loss of precision (because in Cartesian coordinates δx/x is < ∈ when x is small, which it is twice each orbit, whereas in polar coordinates both δr/r and δθ/θ are bounded).

You should be able to calculate thousands of steps per second, it is displaying them that takes the time.

* these are guesses off the top of my head

 

[Drakkith]

Reducing the time step from 43,200 seconds (12 hours) to 432 seconds (7.2 minutes) seems to have a drastic effect. The Earth no longer spirals into the Sun after a few orbits. Now it just drifts slowly outwards!

 

[pbuk]

If you are programming a game, do not waste clock cycles on calculating planetary orbits, just use elliptical formulae. Save your processing power for what is going on in the region of the player(s).

 

[Drakkith]

If you are programming a game, do not waste clock cycles on calculating planetary orbits, just use elliptical formulae. Save your processing power for what is going on in the region of the player(s).

I'm not making a game. I'm just using the engine since it's the only one I'm familiar with.

 

[tfr000]

For what it's worth... I did up the whole solar system in MS Excel. I used a fancy Runge-Kutta-Fehlberg method, even did the relativity corrections to gravity. It worked "ok"... the relativity thing either had problems, or the method just wasn't as accurate as it should have been, because Mercury tended to drift off from where it was supposed to be. So there you go... it's harder than you'd think to get good results with this stuff, even with high-order methods.

 

[jerromyjon]

http://orbitsimulator.com/BA/compareEclipse.html

That is spectacular! I know where I want to be in 30 years.

 

[stedwards]

Several years ago, I used a time step method for the multi-body problem, but the orbits kept growing in diameter. In trying to refine the code, I made an error and transposed two statements. The orbital growth became unnoticable. I could increase the time step by a significant factor, yielding the same amount of orbital error as before.

 

[tfr000]

Several years ago, I used a time step method for the multi-body problem, but the orbits kept growing in diameter. In trying to refine the code, I made an error and transposed two statements. The orbital growth became unnoticable. I could increase the time step by a significant factor, yielding the same amount of orbital error as before.

Yah, I played around with everything - time steps, various terms in the equations of motions, various methods of doing the relativity correction - all to no avail. I was starting to program an "interpolation method" of doing the integration, but ran out of time to work on it...

 

[Kodogo]

I've been tooling around with this, and found that the method known in some quarters as "Leapfrog" gives pretty good results even with modest computing power. The method is the same as described in Feynman Lectures volume 1 chapter 9. It seems to be the same as what some people call Verlet. There's an excellent discussion here:
http://young.physics.ucsc.edu/115/leapfrog.pdf

For a free iPad app to show the results of the method see here:
https://itunes.apple.com/nz/app/orbit-simulator/id1048345753?mt=8

The main problem I see is that highly elliptical orbits do start to precess, depending on the time step you use. Also when the speed gets very high the orbits develop sharp corners. To get round this you's need smaller steps or a method with adaptive step size, but Leapfrog is well good enough just to play around with.