Velocity Verlet leading to faster simulation than Euler in an n-Body simulation?

Question

I have all the constants set to the same values for each set of code, G, the timestep, the masses of the planets etc. But using Velocity Verlet doesn't work unless I lower the gravitational constant because otherwise the gravitational pull is too high. The code is the same for both simulations up until the actual Velocity Verlet integration. Is this normal? If it's not, I'll edit the code into the question as I am unable to do so right now. The simulation uses one massive planet and one less massive planet with a starting speed, which works for symplectic and forward Euler. Here is the code:

/*// Euler Integration
    float deltaT = Extras.TIME_STEP * simSpeed;

    float deltaT2 = deltaT * deltaT;

    float deltaX = planet.getPos().x - this.pos.x;
    float deltaY = planet.getPos().y - this.pos.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 = Extras.G * ((float) this.m * planet.getM()) / (distance * distance);
    this.force.x = F * MathUtils.cosDeg(alpha);
    this.force.y = F * MathUtils.sinDeg(alpha);

    this.accel.x = (float) (this.force.x / this.m);
    this.accel.y = (float) (this.force.y / this.m);

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

    this.pos.x += this.vel.x * deltaT;
    this.pos.y += this.vel.y * deltaT;*/

    // Velocity Verlet Integration
    float deltaT = Extras.TIME_STEP * simSpeed;
    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 = Extras.G * (float) this.m * planet.getM() / (distance * distance);
    this.force.x = F * MathUtils.cosDeg(alpha);
    this.force.y = F * MathUtils.sinDeg(alpha);

    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;

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

    this.accel.x = this.force.x / (float) this.m;
    this.accel.y = this.force.y / (float) this.m;

This is what the orbit looks like in Verlet vs Euler:

Edit: I've found out why, but still not if this is normal: The velocity is updated similarily in both but the position is updated WITH the acceleration in Verlet and not in Euler, which is the correct implementation however.

You can also implement velocity Verlet as v += 0.5*h*a; x += h*v; a = acc(x); v += 0.5*h*a, or use leapfrog Verlet as the main propagation and store the mid-values of the velocities for the "integer indices". — Lutz Lehmann, Feb 26 '23 at 08:53

score 3 · Accepted Answer · answered Feb 24 '23 at 06:57

3

Surprisingly, it is your symplectic Euler code that is wrong. There should be no deltaT2 in any of the Euler updates. Practically this additional factor can be seen as a modification of the force, instead of G you are using G*deltaT, which changes what initial velocities give a numerically stable almost circular orbit.

It is a bad idea to change a numerically sensible time step by the demands of the simulation speed. It would be better to do a number of simspeed steps of the original step size and only use the last step in the animation or plot. Making the step size smaller is usually not a problem. So if you want simspeed=2.5, use simsteps=ceil(simspeed)=3 numerical steps with the slightly smaller step size TIMESTEP*simspeed/simsteps=0.82*TIMESTEP.

Remember that for a circular orbit the initial tangential speed has to be V=sqrt(G*M/R), M the central mass, R the radius of the orbit.

answered Feb 24 '23 at 06:57

Lutz Lehmann

6,044
1
17
25

Good idea with the simspeed as well, thanks! (The orbit isnt necessarily circular its random, the eccentricity turned out to be low in the example i showed) – ght007 Feb 24 '23 at 12:56
Good. Without looking it up, the difference between circular velocity and escape velocity is a factor of 2 or the root of it, so getting the magnitudes wrong can still have drastic consequences. – Lutz Lehmann Feb 24 '23 at 13:22
I edited the code so now it is indeed not using the square. However, now the Euler simulation simply drifts off without completing the first orbit, whereas the Velocity Verlet system achieves stable orbit – ght007 Feb 24 '23 at 15:11
You still have a factor 0.5 in the velocity update that does not belong to the method, thus gets combined into the factor G*M. – Lutz Lehmann Feb 24 '23 at 15:17
Ah yeah, the factor comes from my little knowledge of physics (Newtons equations of motion) but would i be correct in saying it isnt there because the velocity uses a taylor expansion? – ght007 Feb 24 '23 at 15:21
It is much more simple, you have $dv/dt=F/M$ and thus also $\Delta v=(F/m)\Delta t$ in the approximation of the Euler method. Of course this is also the start of the Taylor expansion (or its inversion). But it is only to the linear terms, the factor $1/2$ belongs to the quadratic term. – Lutz Lehmann Feb 24 '23 at 15:32
Let us continue this discussion in chat. – ght007 Feb 24 '23 at 20:35

Velocity Verlet leading to faster simulation than Euler in an n-Body simulation?

1 Answers1