<p>I'm not a Java coder, but it looks like you aren't implementing the velocity Verlet algorithm correctly.</p>
<p>You need to update the velocity vector with the mean of the old acceleration and the new acceleration. The new acceleration must be calculated using the new position vector. But your <code>this.force</code> uses the old <code>this.pos</code>. So <code>this.accel</code> is identical to <code>this.force / this.m</code>.</p>
<hr />
<p>Verlet is a symplectic integrator, which means it conserves energy<sup>1</sup>. So if you're modelling an orbit which ought to be closed, it should close up, not spiral in or out, even with a fairly large time step. So if your orbits are spiralling, that's a sure sign that you've got a bug.</p>
<p>By the way, you don't need those arctan, sin, and cos calls. The acceleration vector is parallel to the radius vector (but of course going in the opposite direction). So you can easily calculate its components from similar triangle ratios, rather than expensive trig function calls.</p>
<p>FWIW, I have code <a href="https://astronomy.stackexchange.com/a/48477/16685">here</a> which plots orbits using synchronised Leapfrog integration. Hopefully, you'll find it readable, even though it's Python, not Java.</p>
<p>The plain synchronised Leapfrog is identical to the velocity Verlet algorithm. In that post I also have a 4th order version of Leapfrog, which is much more accurate.</p>
<ol>
<li>That's a simplification. From <a href="https://en.wikipedia.org/wiki/Symplectic_integrator" rel="nofollow noreferrer">https://en.wikipedia.org/wiki/Symplectic_integrator</a></li>
</ol>
<blockquote>
<p>Symplectic integrators also might possess, as a conserved quantity, a Hamiltonian which is slightly perturbed from the original one (only true for a small class of simple cases).</p>
</blockquote>
<p>It happens to be true in this case. Velocity Verlet / Leapfrog conserves the energy of a system that's approximately equal to the real orbit. We're approximating the true orbit curve with the &quot;fake&quot; Verlet curve. The algorithm preserves the energy of the fake curve. If the time step is small enough, the fake curve is close to the true curve.</p>

CC BY-SA 4.0
2022-07-27T19:37:23.093Z
i tried that but now the orbits are getting smaller, is there still a mistake ( i added the code into the question)
Ght007
CC BY-SA 4.0
2022-07-27T19:51:15.480Z
@Ght007 You need to use the mean of the acceleration calculated using the old position and the acceleration calculated using the new position. Your old code used the mean of the old & the old, which is just the old. Your new code uses the mean of the new and the new.
PM 2Ring
CC BY-SA 4.0
2022-07-27T20:46:41.030Z
answer(s) to my question [What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?](https://scicomp.stackexchange.com/q/29149/17869) are helpful here and reminded me that we shouldn't say things like "symplectic integrator... means it conserves energy".
uhoh
CC BY-SA 4.0
2022-07-27T21:02:22.753Z
thank you for bearing with me! last question i promise: the orbit doesnt get interrupted or changed now but is very slightly deviating from the original one. Im assuming here all the comments about deltaTime are relevant? The problem that would arise then is that i cant have the simulation speed up for the user (but im guessing finding a solution to that isnt physics related and ratehr something for me to figure out myself(
Ght007
CC BY-SA 4.0
2022-07-27T21:40:12.617Z
@Uhoh Ok, I was being a bit lazy, and trying to keep things simple. A symplectic integrator preserves a Hamiltonian that's closely related to the energy.
PM 2Ring
CC BY-SA 4.0
2022-07-27T21:42:27.943Z
@Ght007 Yes, try smaller time steps. Higher eccentricity requires small steps. Try a perfectly circular orbit & see how your code behaves.
PM 2Ring
CC BY-SA 4.0
2022-07-27T21:49:56.727Z
@Ght007 There are various things you can do to speed up your code. Firstly, replace that trig stuff. Your sim doesn't need to plot every point it calculates, so if plotting is relatively slow you can speed things up by (say) plotting every 2nd point. Finally, you can use a larger step size if you use a higher order integrator, but that makes the code more complex. And for highly eccentric orbits, you have to use a fairly small time step.
PM 2Ring
CC BY-SA 4.0
2022-07-27T17:52:01.823Z
2022-07-27T23:33:57.443Z
2022-07-27T23:33:57.443Z
PM 2Ring
<p>I am using the velocity Verlet integration to make an N-Body simulator but the results are odd. If I use simple Newtonian physics I get a closed orbit (btw im testing with 2 planets at first) but if I implement the algorithm the planet leaves the orbit after reaching the closest point to the other planet.
Here is the code:
EDIT: Upon PMRing2's answer I edited the code but for some reason, now the orbit keeps getting smaller.</p>
<pre><code>                float deltaT = Gdx.graphics.getDeltaTime();
        float deltaT2 = deltaT * deltaT;


        float x = this.pos.x;
        float y = this.pos.y;


        x += this.vel.x * deltaT + 0.5f * this.accel.x * deltaT2;
        y += this.vel.y * deltaT + 0.5f * this.accel.y * deltaT2;

        float deltaX = planet.getPos().x - x;
        float deltaY = planet.getPos().y - y;
        float alpha = (float) Math.toDegrees(Math.atan2(deltaY, deltaX));

        float distance = (float) Math.sqrt(Math.pow(deltaX, 2) + Math.pow(deltaY, 2));

        float F = MainScreen.G * this.m * planet.getM() / (distance * distance);
        this.force.x = F * MathUtils.cosDeg(alpha);
        this.force.y = F * MathUtils.sinDeg(alpha);

        this.accel.x = this.force.x / this.m;
        this.accel.y = this.force.y / this.m;

        this.pos.x = x;
        this.pos.y = y;

        this.vel.x += 0.5f * (this.accel.x + (this.force.x / this.m)) * deltaT;
        this.vel.y += 0.5f * (this.accel.y + (this.force.y / this.m)) * deltaT;
</code></pre>
<p>Most implementations or explanatory videos I've watched don't use the force so I fear that might play a role. I am also not sure if the implementation is even right (I am quite new to this).</p>

CC BY-SA 4.0
2022-07-27T12:37:04.477Z
Default advice: try lowering your `deltaT`
Barry Carter
CC BY-SA 4.0
2022-07-27T13:30:45.097Z
it did something but not a solution for the whole thing as the orbit is still not fixed. the orbit got bigger and bigger. theres a different problem i think
Ght007
CC BY-SA 4.0
2022-07-27T13:35:32.657Z
Experiment with different `deltaT`'s. I think you'll find smaller values yield more accurate results. You're effectively using a numerical method to solve a differential equation, and there's a lot of things to consider
Barry Carter
CC BY-SA 4.0
2022-07-27T14:09:41.323Z
Can you provide a working example for your code, ideally something minimalised for testing, and which we can run in an online compiler or similar. Creating a minimal working example might take a bit of time, but it is time well spent.  Often the process of mimimising shows up where the bugs are.
James K
CC BY-SA 4.0
2022-07-27T15:08:02.913Z
I am not sure if I am able to do that as the code is using a game engine in Java and I'm not fluent in JavaScript to have it on the web etc. But it would be helpful to know whether my implementation is correct as if it is, I can simply debug it myself and have this question deleted due to possible uselessness for future visitors
Ght007
CC BY-SA 4.0
2022-07-27T15:08:41.233Z
Not sooo surprising that a numerical method/situation is "unstable", in exactly the sense you've discovered. That could be an inescapable aspect of the physical situation, so that taking really tiny increments is the only improvement, but, even then, it will only be a matter of degree. Or, it might just reflect a mathematical issue, not necessarily physical. This type of problem does illustrate why "minimization/extremal" approaches (like Lagrangian mechanics) are more stable mathematically.
paul garrett
CC BY-SA 4.0
2022-07-27T18:14:23.427Z
@paulgarrett Verlet / Leapfrog is quite stable for plotting simple orbits. Sure, it helps to use a small time step, but it's a [symplectic integrator](https://en.wikipedia.org/wiki/Symplectic_integrator), so it doesn't blow up if the time step isn't small enough.
PM 2Ring
CC BY-SA 4.0
2022-07-27T18:16:48.757Z
@PM2Ring, ah, thanks for the info! I'm only an amateur in such things. :)
paul garrett
CC BY-SA 4.0
2022-07-27T18:29:09.157Z
As others have said, the step size per iteration is a common cause of errors.  Too large and you miss important forces, too small and rounding becomes an issue.  You don't specify the orbit being modeled, but if there's a huge different between apoapsis and periapsis, it's going to be a lot harder to get right.  So try with a mostly circular orbit first.
Greg Miller
<orbit><orbital-mechanics><eccentric-orbit>
Velocity Verlet algorithm creating odd orbits
Ght007
# Velocity Verlet algorithm creating odd orbits

(Asked by: username_0 on 2022-07-27)

I am using the velocity Verlet integration to make an N-Body simulator but the results are odd. If I use simple Newtonian physics I get a closed orbit (btw im testing with 2 planets at first) but if I implement the algorithm the planet leaves the orbit after reaching the closest point to the other planet.
Here is the code:
EDIT: Upon PMRing2's answer I edited the code but for some reason, now the orbit keeps getting smaller.
                float deltaT = Gdx.graphics.getDeltaTime();
        float deltaT2 = deltaT * deltaT;


        float x = this.pos.x;
        float y = this.pos.y;


        x += this.vel.x * deltaT + 0.5f * this.accel.x * deltaT2;
        y += this.vel.y * deltaT + 0.5f * this.accel.y * deltaT2;

        float deltaX = planet.getPos().x - x;
        float deltaY = planet.getPos().y - y;
        float alpha = (float) Math.toDegrees(Math.atan2(deltaY, deltaX));

        float distance = (float) Math.sqrt(Math.pow(deltaX, 2) + Math.pow(deltaY, 2));

        float F = MainScreen.G * this.m * planet.getM() / (distance * distance);
        this.force.x = F * MathUtils.cosDeg(alpha);
        this.force.y = F * MathUtils.sinDeg(alpha);

        this.accel.x = this.force.x / this.m;
        this.accel.y = this.force.y / this.m;

        this.pos.x = x;
        this.pos.y = y;

        this.vel.x += 0.5f * (this.accel.x + (this.force.x / this.m)) * deltaT;
        this.vel.y += 0.5f * (this.accel.y + (this.force.y / this.m)) * deltaT;

Most implementations or explanatory videos I've watched don't use the force so I fear that might play a role. I am also not sure if the implementation is even right (I am quite new to this).

**Tags:** orbit, orbital-mechanics, eccentric-orbit
> **Comments:**
> **username_1**: Default advice: try lowering your `deltaT`  
> **username_0**: it did something but not a solution for the whole thing as the orbit is still not fixed. the orbit got bigger and bigger. theres a different problem i think  
> **username_1**: Experiment with different `deltaT`'s. I think you'll find smaller values yield more accurate results. You're effectively using a numerical method to solve a differential equation, and there's a lot of things to consider  
> **username_2**: Can you provide a working example for your code, ideally something minimalised for testing, and which we can run in an online compiler or similar. Creating a minimal working example might take a bit of time, but it is time well spent.  Often the process of mimimising shows up where the bugs are.  
> **username_0**: I am not sure if I am able to do that as the code is using a game engine in Java and I'm not fluent in JavaScript to have it on the web etc. But it would be helpful to know whether my implementation is correct as if it is, I can simply debug it myself and have this question deleted due to possible uselessness for future visitors  
> **username_3**: Not sooo surprising that a numerical method/situation is "unstable", in exactly the sense you've discovered. That could be an inescapable aspect of the physical situation, so that taking really tiny increments is the only improvement, but, even then, it will only be a matter of degree. Or, it might just reflect a mathematical issue, not necessarily physical. This type of problem does illustrate why "minimization/extremal" approaches (like Lagrangian mechanics) are more stable mathematically.  
> **username_4**: @paulgarrett Verlet / Leapfrog is quite stable for plotting simple orbits. Sure, it helps to use a small time step, but it's a [symplectic integrator](https://en.wikipedia.org/wiki/Symplectic_integrator), so it doesn't blow up if the time step isn't small enough.  
> **username_3**: @PM2Ring, ah, thanks for the info! I'm only an amateur in such things. :)  
> **username_5**: As others have said, the step size per iteration is a common cause of errors.  Too large and you miss important forces, too small and rounding becomes an issue.  You don't specify the orbit being modeled, but if there's a huge different between apoapsis and periapsis, it's going to be a lot harder to get right.  So try with a mostly circular orbit first.  

---

## 1 answers

### Answer by username_4 on 2022-07-27. Score: 5 | Accepted answer

I'm not a Java coder, but it looks like you aren't implementing the velocity Verlet algorithm correctly.

You need to update the velocity vector with the mean of the old acceleration and the new acceleration. The new acceleration must be calculated using the new position vector. But your `this.force` uses the old `this.pos`. So `this.accel` is identical to `this.force / this.m`.

---

Verlet is a symplectic integrator, which means it conserves energy1. So if you're modelling an orbit which ought to be closed, it should close up, not spiral in or out, even with a fairly large time step. So if your orbits are spiralling, that's a sure sign that you've got a bug.

By the way, you don't need those arctan, sin, and cos calls. The acceleration vector is parallel to the radius vector (but of course going in the opposite direction). So you can easily calculate its components from similar triangle ratios, rather than expensive trig function calls.

FWIW, I have code [here](https://astronomy.stackexchange.com/a/48477/16685) which plots orbits using synchronised Leapfrog integration. Hopefully, you'll find it readable, even though it's Python, not Java.

The plain synchronised Leapfrog is identical to the velocity Verlet algorithm. In that post I also have a 4th order version of Leapfrog, which is much more accurate.

1. That's a simplification. From <https://en.wikipedia.org/wiki/Symplectic_integrator>

> Symplectic integrators also might possess, as a conserved quantity, a Hamiltonian which is slightly perturbed from the original one (only true for a small class of simple cases).

It happens to be true in this case. Velocity Verlet / Leapfrog conserves the energy of a system that's approximately equal to the real orbit. We're approximating the true orbit curve with the "fake" Verlet curve. The algorithm preserves the energy of the fake curve. If the time step is small enough, the fake curve is close to the true curve.

> **Comments:**
> **username_0**: i tried that but now the orbits are getting smaller, is there still a mistake ( i added the code into the question)  
> **username_4**: @Ght007 You need to use the mean of the acceleration calculated using the old position and the acceleration calculated using the new position. Your old code used the mean of the old & the old, which is just the old. Your new code uses the mean of the new and the new.  
> **username_6**: answer(s) to my question [What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?](https://scicomp.stackexchange.com/q/29149/17869) are helpful here and reminded me that we shouldn't say things like "symplectic integrator... means it conserves energy".  
> **username_0**: thank you for bearing with me! last question i promise: the orbit doesnt get interrupted or changed now but is very slightly deviating from the original one. Im assuming here all the comments about deltaTime are relevant? The problem that would arise then is that i cant have the simulation speed up for the user (but im guessing finding a solution to that isnt physics related and ratehr something for me to figure out myself(  
> **username_4**: @Uhoh Ok, I was being a bit lazy, and trying to keep things simple. A symplectic integrator preserves a Hamiltonian that's closely related to the energy.  
> **username_4**: @Ght007 Yes, try smaller time steps. Higher eccentricity requires small steps. Try a perfectly circular orbit & see how your code behaves.  
> **username_4**: @Ght007 There are various things you can do to speed up your code. Firstly, replace that trig stuff. Your sim doesn't need to plot every point it calculates, so if plotting is relatively slow you can speed things up by (say) plotting every 2nd point. Finally, you can use a larger step size if you use a higher order integrator, but that makes the code more complex. And for highly eccentric orbits, you have to use a fairly small time step.  

---

