Orbits

Game Mechanics
#21

@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…

#22

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

#23

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.

#24

Ah! I’d played that before. Hadn’t realised it was one of yours! :smiley:

#25

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

#26

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.

#27

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… >___>

#28

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.

#29

I don’t wanna do any hardcore simulations. I just don’t want my subtractions to wobble like jelly. xd