Difference in calculated and simulated ellipsies

Question

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.

Answer 1

Here's your key error:

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.

Answer 2

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.