# Difference in calculated and simulated ellipsies

My task here is to determine orbit parameters, using current values:

1. $\mu=GM$ - standard gravitational parameter
2. $r$ - distance to the object with Mass $M$
3. $v$ - speed of the object in the point $r$

I have tried to do this task using Specific orbital energy:

$$\epsilon=\frac{v^2}{2}-\frac{\mu}{r} =-\frac{\mu^2}{2h^2}(1-e^2)=-\frac{\mu}{2a};(1)$$ This gives me $a$ and $e$. $$a = -\frac{\mu}{2\epsilon};$$ $$e = \sqrt{1+\frac{2h^2\epsilon}{\mu^2}}; h=\vec{r}\times \vec{v};$$ In 2d cross-product of vectors are they's manitude;

At the same time, I've launched simulation to draw gravitation only movement with the same $\mu$ as in the given equation.

Simulation is simple: $$\vec g=\frac{\mu}{\vec {r^2}};(2)$$ $$\vec v=\vec gt+\vec {v_0};(3)$$ $$\vec x=\vec gt^2 + \vec vt+\vec{x_0};(4)$$

The code in simulation is written in javascript:

        var rv = Sub(o.p, sun.c); // o.p contains [x,y] coords of point, and sun.c contains coords of sun - gravity object
var rvn = Normal(rv); // position normal
var g =  mu/ (r*r);  // free-fall acceleration magnitude

var gv = Mul(rvn, -g); // frefall vector - negative to position

var dv = Mul(gv, dt); //velocity changes

var dx = Add(Mul(gv , dt*dt) , Mul(o.v, dt)); // position changes

o.v = Add(o.v, dv); // new speed
o.p = Add(o.p, dx); // new position


But, my code gives me wrong results. Ellipse, which is drawn on calculated from (1) parameters ($e$,$a$) is narrower than ellipse, which is drawn by my gravity simulation.

...And in the same way, when I took $r_a$ and $r_p$ parameters from my simulation and drawn an ellipse with $e=\frac{r_a-r_p}{r_a+r_p}$ and $a=\frac{r_a+r_p}{2}$, It was corresponding to calculated (not simulated) ellipse. This brings me to make assumption, that there's something wrong with gravity simulator, but I'm running out of ideas, what's going on with my simulation.

Where could I make a mistake?

PS: Demo code is here. It uses mouse to shoot moving objects. After shooting all ellipses are drawn, but there's several options, when calculated ellipses having wrong orientation - reload page to remove them and start with blank canvas. When shooting objects with eccentricity more than 1, arrow showed in red.

• Before I look too closely at your question, could you clear up two things? 1.) You don't specify the direction of the velocity. Are you given the fact that $\vec{v}$ is perpendicular to $\vec{r}$? 2.) Shouldn't there be a factor of 1/2 in the first term of equation (4)? – garyp Jul 1 '14 at 12:16
• (1)Direction of velocity is specified by ortogonal vector V=[x,y]. equations (2),(3),(4) should be written in vector form of course, I'll edit my question now to point this. But in the source code Mul, Sub, etc. functions - are function over vectors. – Vasiliy Stavenko Jul 1 '14 at 12:23
• (2)I thought, that this factor should be there. But, when I tried to place it, my gravitational traectories become elliptic extending spirals. When this factor removed, I can get elliptical traectories, which remaining solid after 100+ cycles. – Vasiliy Stavenko Jul 1 '14 at 12:26
• Could you clarify, is it necessary to have $\vec v$ perdendecular to $\vec r$ for elliptical orbit? – Vasiliy Stavenko Jul 1 '14 at 13:02
• The problem is your equation (4), $\vec x(t) = \vec g t^2 + \vec v t + \vec x_0$. Remember your high school physics! That should be $\vec x(t) = \frac 1 2 \vec g t^2 + \vec v t + \vec x_0$. See my answer. – David Hammen Jul 1 '14 at 14:28

Even though your integration method seems wrong, like David Hammen pointed out, the results look do not look wrong. The problem rely lies with the way you define your ellipse. Your definition for the semi-major axis, $a$, and eccentricity, $e$, are correct. The way you define the semi-minor axis, $b$, probably is not, which is defined as: $$b = a \sqrt{1 - e^2}$$ A quick measure of your image reveals that your numerical results agree with this and the ellipse you have drawn does not.

However your goal is to determine orbital parameters. You might want to check out this question, since it is related to that.

• drawEllipse() contains the line b=a * Math.sqrt(1 - o.e) that looks suspicious... – Christoph Jul 1 '14 at 21:03
• I took this equation from en.wikipedia.org/wiki/Ellipse#Eccentricity. $e=\sqrt{1-\frac{b^2}{a^2}}$. You could take $b$ from it easely. – Vasiliy Stavenko Jul 1 '14 at 22:18
• @VasiliyStavenko - Look at line 161 of your velocity_based.html. – David Hammen Jul 1 '14 at 22:36

var dx = Add(Mul(gv , dt*dt) , Mul(o.v, dt)); // position changes


In math, that's $\Delta \vec x = \vec g \Delta t^2 + \vec v \Delta t$. That's incorrect. Assuming a constant acceleration $\vec g$, the change in position over some time interval $\Delta t$ is given by $\Delta x = \frac 1 2 \vec g \Delta t^2 + \vec v \Delta t$. That line of code should be

var dx = Add(Mul(gv , 0.5*dt*dt) , Mul(o.v, dt)); // position changes


This will get rid of your huge error, but not all of the error. There's no way to get rid of all of the error with numerical integration (but you can come close with advanced techniques).

You are still going to suffer errors because your position and velocity updates assume constant acceleration. Is there a simple way to reduce those errors? The answer is yes. A simple way to significantly improve both accuracy and stability is to change the velocity calculation

var dv = Mul(gv, dt); //velocity changes


to

var dx = Add(Mul(gv , 0.5*dt*dt) , Mul(o.v, dt)); // position changes
o.p = Add(o.p, dx); // new position

var dv1 = Mul(gv, 0.5*dt); // half-step velocity changes
// Calculate gv using updated position; code not shown (you already have it)
// Note: You can use this updated gv during the next step
var dv2 = Mul(gv, 0.5*dt); // half-step velocity changes
var dv = Add (dv1, dv2); // full step velocity changes


This is a variant of the leapfrog numerical integration technique. You are still going to suffer errors that grow over time, but the errors will be drastically reduced from even the simple correction I offered at first.

Simple techniques such as Euler's method, verlet integration, leapfrog, and Runge Kutta techniques will only take you so far. Eventually accuracy and stability will suffer. High accuracy numerical orbit integrators for solar systems use much more complex numerical integration techniques.

Galactic scale simulators on the other hand oftentimes do use simple techniques such as leapfrog because simulating tens of thousands of stars (or more!) is a very daunting task computationally. One big advantage of the leapfrog technique is that it is "symplectic", which means that it comes close to conserving energy and angular momentum. This is very important for galactic scale simulations. The classical Runge Kutta techniques are not symplectic, so orbits tend to spiral in or spiral out, depending on the technique and the initial conditions.

• That's sad! :). I could realise that myself... But at the same time, how perfect error got out just in case, when you removing $1/2$ factor. Awesome. – Vasiliy Stavenko Jul 1 '14 at 22:21
• I have checked simulation. I takes effect, but the most meaning thing was, that I was incorrectly calculated my minor semi axis. – Vasiliy Stavenko Jul 1 '14 at 23:13