Thank you for the source, Pjt33! I’ll try to use it tomorrow…
Concerning Semi-implicit Euler: What is the difference between it and normal Euler? I only had numerical integration mentioned to me in high school, so I can’t really understand the mathy stuff on Wiki… T___T
Update the velocity using the gravitational forces at the current position. Then update the position using the new value of the velocity.
I’d be curious to know whether that makes any noticeable difference to your orbits compared to Euler with the same time step.
Doing a bit more reading, Verlet integration has the same energy-conservation property as semi-implicit Euler, and is of higher order, so that’s another thing to try.
Simon
So semi-implicit Euler is just a fancy name for this code?
double dx = SUN_X - x, dy = SUN_Y - y;
double distSqrd = dx*dx + dy*dy;
double a = GRAVITATIONAL_CONSTANT * SUN_MASS / (distSqrd);
double dist = Math.sqrt(distSqrd);
velX += dx * a / dist * TIMESTEP;
velY += dy * a / dist * TIMESTEP;
x += velX * TIMESTEP;
y += velY * TIMESTEP;
That’s what I’m doing at the moment… Is normal Euler updating the position before the velocity?
Yes, that looks right. If you switch the “velX/Y += …” lines with the “x/y += …” lines then you’d have standard Euler.
Does it make any difference?
Simon
Made some test code to allow me to easily test things faster. So far I’ve implemented Euler, Semi-implicit Euler and Verlet (the Velocity Verlet from Wikipedia).
- Euler: Sucks ballz. Explodes with somewhat high acceleration. Worthless for anything but exactly round orbits.
- Semi-implicit Euler: Much better than standard Euler. Stable orbits, but the whole orbit is rotating. (see my last post)
- Verlet: The best one by a small margin. Marginally more precise orbit compared to Semi-implicit Euler, but the exact same rotation problem. Handles extremes slightly better (only explodes at higher accelerations).
Verlet and Semi-implicit Euler followed each other almost pixel exactly for over 100k iterations, so I assume there is no significant error buildup, except for extremes where Verlet marginally survives and Euler explodes. Not really that relevant as that would just be a collision anyways, but as their computation cost is about the same, there is no reason to not use Verlet.
I still have the problem with rotating orbits. This is the most important problem to solve at the moment. I will try to implement RK4. I’ll be back, either for help or to show my results.
For some reason I can’t seem to wrap my mind around RK4 for my problem. As I don’t have a function depending on time, I don’t see how I should go around implementing it… T___T I realize I should just integrate x and y separately, but I can’t seem to connect f(t) to anything decent for my problem. Gah, so confused.
I did consider refactoring to remove the time parameter from the code, but I decided against. You can just ignore it.
dx/dt = vx
dy/dt = vy
dvx/dt = ax
dvy/dt = ay
You must be using those for the other quadrature implementations, so it’s just a case of putting them in a Vec4 and returning it from the eval method.
Edit: in fact, you’ve posted code further up. So
class MyFn implements Fn {
public Vec4 eval(double t, Vec4 pos) {
double dx = SUN_X - pos.w, dy = SUN_Y - pos.x;
double distSqrd = dx*dx + dy*dy;
double a = GRAVITATIONAL_CONSTANT * SUN_MASS / (distSqrd);
double dist = Math.sqrt(distSqrd);
return new Vec4(pos.y, pos.z, dx * a / dist, dy * a / dist);
}
}
Currently I have a Body class that has a position, velocity and a stored acceleration (needed for Verlet at least). I pass this body to an integrator using the method integrate(Body b, double timestep), which using the function getAcceleration(double x, double y) updates the body’s data. Your eval would be the same to my getAcceleration, and integrate would be takeStep, correct?
What I don’t understand in your code is how you use your Vec4s. For example:
double dx = SUN_X - pos.w, dy = SUN_Y - pos.x;
What’s up with pos.w and pos.x?
I thought it would be educational (for me if no one else) to look at how some different integration methods perform for different step sizes.
So I set up a simple simulation of an asteroid following an elliptical orbit around a star, and ran it for a range of different sizes of time step. The simulations cover the length of time that the asteroid would take to complete one period of the orbit in the exact solution. So, if a simulation were perfectly accurate, the asteroid would end up exactly back at its starting point. I’m taking the distance of the asteroid from its starting point at the end of the simulation as the measure of a simulation’s error.
The asteroid starts at distance 1.0 from the star, and that’s also its maximum distance (for the exact orbit). So an error of 0.01 at the end of the simulation would probably be acceptable, although of course we’d like to do better than that. If the error is greater than 1 then the simulation was a complete mess.
The simulations run for different numbers of time steps. If T is the time simulated (i.e., one period of the exact orbit) and the simulation performs N time steps, then the size of each time step is dt=T/N.
Here are some plots showing the final error when different numbers of time steps are used. The plots cover the Euler method, the semi-implicit Euler method, and the fourth-order Runge-Kutta method. I also tried the basic Verlet and velocity Verlet methods, but the results were so similar to those of the semi-implicit Euler method that they’re not worth showing.
http://www.dishmoth.com/fig1.png
http://www.dishmoth.com/fig2.png
http://www.dishmoth.com/fig3.png
Clearly the Euler method needs a very large number of very small time steps in order for the error to be small. The semi-implicit Euler method gets away with a lot fewer time steps, and the Runge-Kutta method fewer still.
The difference in performance is clear when the error graphs are plotted all together with logarithmic axes.
http://www.dishmoth.com/fig4.png
For this plot it’s interesting to note that, when the number of time steps gets large enough, the slope of the Runge-Kutta curve is roughly -4, the slope of the semi-implicit Euler curve is roughly -2, and the slope of the Euler curve is roughly -1. This is evidence that the methods converge at fourth order, second order and first order respectively. Curiously, the semi-implicit Euler method is, like the Euler method, only supposed to converge at first order (whereas the Verlet methods converge at second order) so for this problem at least it’s performing better than advertised.
In conclusion, the Euler method isn’t worth bothering with for this problem. The semi-implicit Euler method (which is almost identical) is much better, but in terms of accuracy the Runge-Kutta method wins. That said, it takes a few times more work to perform each step of the Runge-Kutta method, and it’s harder to code, so unless you’re chasing very high accuracy it may not be worth the effort. (In effect, the green curve in the logarithmic plot should be shifted a bit to the right to account for the extra computational work involved.)
All this is just for one example orbit of course. The number of time steps needed to get decent accuracy will be different for different orbits, and will probably depend quite a lot on how close the orbit gets to the star.
I do some strange things for fun sometimes…
Simon
The Vec4’s components are called w,x,y,z for generality. I’m using the same convention as in the Trajectory example class that position is (w,x) and velocity is (y,z). If you want to refactor w,x,y,z to x,y,vx,vy, feel free.
eval isn’t quite the same as getAcceleration. It gets the derivatives of all the components, which when the components are speed and velocity would be velocity and acceleration.
Nice work.
As a point of interest, my FunOrb game Torquing uses RK4 on a Vec13. IIRC the 13 components are 3 for position, 3 for velocity, 4 for angular position (using quaternions, and renormalising them regularly), and 3 for angular velocity (and I can’t remember why that isn’t using quaternions).
It’s not a proper plug if you don’t give a link.
It was a plug for RK4 rather than Torquing, but if you insist…
http://www.funorb.com/info.ws?game=torquing
@Dishmoth
Awesome! Very interesting! I will definitely implement RK4 to try it too, it seems like it’s better for smaller time steps but worse for larger (-.-) time steps. It would be interesting to see the difference in precision in my case.
@Pjt33
Thank you!
Am I correct in thinking that the eval() function should return a Vec4 with (vx, vy, ax, ay) where ax, ay is the gravitational force at the position?
Also there seems to be a slight problem with your takeStep(). Currently your test prints:
finalX = (2.0000000000000013, -3.2999999999999656, 1.0, -9.999999999999963)
Shouldn’t Y be close to 0, not -3.3…?
I believe I found the error, k3 should be scaled too:
Vec4 k1 = fn.eval(X);
Vec4 k2 = fn.eval(X.addScaled(k1, h2));
Vec4 k3 = fn.eval(X.addScaled(k2, h2));
//Vec4 k4 = fn.eval(X.add(k3));
Vec4 k4 = fn.eval(X.addScaled(k3, h));
This gives
finalX = (2.0000000000000013, 3.952393967665557E-14, 1.0, -9.999999999999963)
which is a little closer to 0 than -3.3. Hehe. It doesn’t really make sense to ignore the time step there, and h produces a better result compared to h2.
I feel like I understand this a little better now, but I still can’t get it to work with gravity. I’m suspecting a bug in my code at the moment…
I GOT IT WORKING!!! I was sending vx and vy to getAcceleration()! xD
RK4 is very different from Verlet Velocity and Semi-implicit Euler. It does not rotate at all, which is awesome. Instead it loses energy if the time step is too large, eventually crashing into the star and exploding.
However! This is kind of okay! I could just blame it on atmospheric drag or something, as it only occurs when the asteroid is close to the body.
I’ve yet to do performance tests and optimizations, but I believe I will use RK4 from now on!
Thank you so much guys! Physics owns! xD
You’re right about that bug. Good catch. I introduced it when refactoring slightly to use the addScale method, which wasn’t in my previous Vec4 implementation.

It was a plug for RK4 rather than Torquing, but if you insist…
http://www.funorb.com/info.ws?game=torquing
Ah! I’d played that before. Hadn’t realised it was one of yours!
I optimized away all object creation in my code. I then did a small performance test on Semi-explicit Euler, Verlet and RK4. Results:
100.000 objects, 1 attractor
Semi-implicit Euler: 4.5 ms
Verlet: 4.6 ms
RK4: 10.2 ms
100.000 objects, 10 attractors
Semi-implicit Euler: 17.0 ms
Verlet: 17.0 ms
RK4: 68.5 ms
Verlet is insignificantly slower than Semi-implicit Euler, but the precision difference is probably less than significant too. xD When the number of attractors increases RK4 falls behind. The actual operations for each integrator is not very costly, even for RK4. The bottleneck seems to be my gravitational force calculations. RK4 gets really slow because it does 4 of these each step.
This would’ve been fine if RK4 had 4 times the precision of Verlet, but unfortunately RK4’s precision gets better with smaller time steps. For small time steps it is awesome, almost perfect, but for larger time steps it is equally or even more inaccurate than Verlet.
They also behave very differently. If your time step is too large and you have an elliptical orbit, the whole orbit will spin around the attractor, making a cute flower pattern. RK4 however just loses energy eventually crashing into the attractor and then exploding. It seems like that would be an important point when deciding on what integrator to use for a game. A space game could use RK4 and justify the energy loss with atmospheric friction (the reason satellites eventually enter the atmosphere), while a puzzle game or a particle engine using magnets or “black holes” or something would be better of with Verlet, as it performs better and the spinning might not be noticeable if objects do not orbit. Space games also usually have a lot larger time steps, you probably don’t want to simulate an entire orbital period in 1/60 second time steps. xD
Precision loss occurs when objects get close to attractors, so the shape of the orbit definitely plays a part in the decision. If you just have circular orbits, RK4 is waaay overkill.
Does anybody have any ideas on how to improve my gravitational acceleration function?
/*Acceleration is just an object containing the x- and y-accelerations (ax, ay). Supplied by the integrator to avoid having to create a new one each run. (Verlet and Euler keeps 1 instance, RK4 has 4.)*/
public void getAcceleration(Acceleration a, double x, double y) {
a.ax = 0;
a.ay = 0;
for(int i = 0; i < gravityBodies.size(); i++){
GravityBody p = gravityBodies.get(i);
double dx = p.x - x, dy = p.y - y;
double distSqrd = dx*dx + dy*dy;
double scale = GRAVITATIONAL_CONSTANT * p.mass / (distSqrd * Math.sqrt(distSqrd)); //Calculates the acceleration (G*M/r^2) and divides it by the distance.
a.ax += dx*scale;
a.ay += dy*scale;
}
}
And yes, GravityBody is a really bad name. xD
Stable Newtonian dynamics integration is in fact a big area. Why? Because anything that is not 2 body’s is not stable long term in reality. So in Astronomy questions like how long is earth orbit stable is an interesting, important and hard to answer question. Since they want to predict the orbits thousands of years into the future, these methods are probably not really applicable here ;).
But as others have kinda of stated but not really made explicit, is that you want a method that preserves time reversibility (momentum and energy kinda with roundoff error). The classic example is the leap frog method. This is perfect in the sense that you can always go back a step as well as forward. Note this only works for nice gravity problems where the acceleration is a function of positions only. Once acceleration is dependent on velocity there are no “really good methods” but its not just a matter of the order of the method. How it bleeds energy and momentum via round off are important.
It should be noted that RK4 is not time reversible, and is often not used for these types of simulations.
Personally I would start with a leap frog or a Velocity Verlet method. I would never bother with RK4. Engineers love to bash that method to death, but it is overrated IMO.
There is another solution if you are only considering 2 body. Use the well developed analytical orbital equations. No good for many body case however.
ps we are pretty sure that earths orbit is going to be stable for at least a millions years or so… but even a 10 meter error in current measurements means you can’t predict further out than about 10 millions years IIRC.
Hello! Sorry to resurrect this thread again, but I completely missed Delt0r’s post!
Further testing has led me to believe that RK4 is impractical because it loses energy. Verlet just spins, so the satellite will at least probably never crash into whatever it’s orbiting. I’ll just have to keep my orbits as round as possible.
After Delt0r’s post I made a test with simplified gravity. Only the smallest body gets pulled (the sun pulls the Earth, but not the other way around). I put the sun fixed at (0,0). I then put the Earth in an orbit around the sun, then the moon around the Earth and finally a satellite around the moon. All positions and velocities were based on scientific real world values except my moon satellite, which I based on a newspaper article. I got “stable” non-exploding orbits with a time step of 12 seconds using Verlet Velocity.
Verlet Velocity has so much better performance than RK4 and seems to work just fine. I will probably use it if I decide to make a game out of this. Implementing RK4 was fun though, and I think it’s great to know how it works.
I’m a little worried about floating point imprecision. The addition on position and velocity doesn’t seem to be that bad, but the distance calculation for gravity is really bouncy for my satellite. It’s because of the subtractions (dx = moon.x - satellite.x;). I don’t know how much this will affect stuff, but if I’m gonna have multiple solar systems, I definitely need a local coordinate system for each system, probably with the sun at (0, 0).
EDIT: Yep, stable for 817 years… >___>
The problem with solar system dynamics and roundoff is how nonlinear roundoff is with floats. That is round off is a function of the number. Switching to fixed point (aka integer math) does help quite a bit IME (In My Experience). For example with longs at 1mm per unit, 2^63 is still about a light year. So its pretty manageable.
For real hard core simulations they use baselines and local coordinates a lot IIRC. Bear in mind that hard core simulations need to take many things into account that you don’t care about, like the Yarkovsky effect, photon pressure, the large asteroids so on and so forth.