Java-Gaming.org Hi !
 Featured games (84) games approved by the League of Dukes Games in Showcase (574) Games in Android Showcase (154) games submitted by our members Games in WIP (620) games currently in development
 News: Read the Java Gaming Resources, or peek at the official Java tutorials
Pages: [1]
 ignore  |  Print
 Orbits  (Read 4363 times) 0 Members and 1 Guest are viewing this topic.
theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Posted 2011-08-04 03:51:20 »

Hello, I'm kinda experimenting a little with gravitational forces and orbits in 2D. I have some objects (asteroids or something) being influenced by n static bodies (planet, stars). 2 questions:
1. What integrator should I use? Euler sucks, it requires too small time steps to be useful, and performance is fairly important. I would really like a code example for gravitation specifically. I've found numerous articles on through Google on integration, but I didn't really understand the actual implementations... T___T
2. How do I calculate the gravitational acceleration? With the standard G*M/R^2, I just get an acceleration, not an 2D acceleration vector. I currently get the angle between the 2 bodies using atan2(dy, dx) and then get the vector using cos and sin, but it is really slow. I suppose I could compute a direction vector by normalizing (dx, dy), but that would require a sqrt and multiple divides... Is there any faster way? It would be awesome if there was an equation that actually produced a vector in the first place.
I may be using some words wrong here, so please don't flame me. xD

Myomyomyo.
ra4king

JGO Kernel

Medals: 374
Projects: 3
Exp: 5 years

I'm the King!

 « Reply #1 - Posted 2011-08-04 05:54:25 »

1. Can't help much, don't know what integration or Euler is :S

2. Faster sin/cos using lookup table: http://www.java-gaming.org/topics/fast-math-sin-cos-lookup-tables/24191/view.html
Faster atan2 also using lookup table: http://www.java-gaming.org/topics/fast-math-atan2-lookup-table/24193/view.html

endolf

JGO Coder

Medals: 7
Exp: 15 years

Current project release date: sometime in 3003

 « Reply #2 - Posted 2011-08-04 06:40:24 »

The acceleration from gravity is the length of the acceleration vector, to get the vector, just subtract the position of the object you are being accelerated towards from the object you are accelerating. Then normalise that vector, that gives you a unit vector describing the direction of acceleration, then you can multiply it by the acceleration from gravity and you get the acceleration vector.

then you can just scale a copy of that by time to get the change in velocity for each frame, then apply that to the current velocity vector for the object, then multiply that by time to get the change in position for the frame. Apply to the transform and job done

Add in angular acceleration + velocity etc and all is nice .

It's what I'm doing in darkvoid which is the full 3D 6dof .

Endolf

pitbuller
 « Reply #3 - Posted 2011-08-04 07:26:43 »

http://www.java-gaming.org/topics/gravitational-slingshot/24085/view.html

There is code for solution with  component calculations if you don't want use vectors. Its quite simple and efficient too. Some cleaning is needed for more general use but you should see how to use it.

 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23 `public final static long staticGravityForce = 1000; //just test what is goodpublic final static long staticX = 400;public final static long staticY = 400;public final static long staticRadius = 50;public void gravity(double x, double y, int radius, double dx,double dy) {          double ddx = (x - staticX);    double ddy = (y - staticY);    double distanceSquared = ddx * ddx + ddy * ddy;    double distance = Math.sqrt(distanceSquared);     if ((distance) > (staticRadius + radius)) {                 ddx /= distance;        ddy /= distance;        dx += (ddx * staticGravityForce / distanceSquared);        dy += (ddy * staticGravityForce / distanceSquared);                 } else {        //COLLAPSE        //JUST make some cool effects here    }    }}`
dishmoth
 « Reply #4 - Posted 2011-08-04 08:38:34 »

Quote
1. What integrator should I use? Euler sucks, it requires too small time steps to be useful
I wouldn't be too quick to dismiss Euler integration.  I'd expect it to be perfectly adequate for games (if that's what you're doing), and certainly for prototyping.  Make the time steps small if you have to, but get everything working before trying to optimize.

But if you're set on using higher-order methods, try the chapter on "Integration of Ordinary Differential Equations" in Numerical Recipes.  Runge-Kutta methods are usually a good starting point.

There presumably are methods that are specific to gravitation, but I'm not familiar with any.  (You could maybe try tweaking the kinetic energy of the asteroids so that kinetic plus potential is conserved during each step.)

Quote
2. I suppose I could compute a direction vector by normalizing (dx, dy), but that would require a sqrt and multiple divides...
Normalizing the vector (sqrt) makes more sense than dealing with angles (sin and cos and atan2).  Also, knowing the (non-squared) distance of the asteroid from the planet could turn out to be useful in other parts of your simulation.

Quote
I have some objects (asteroids or something) being influenced by n static bodies (planet, stars).
Do the asteroids influence each other gravitationally?  If not then the gravitational field of the planets is static, and you can pre-compute it.  That is, compute the force vector (or potential value) at each pixel on the screen (or over some other grid) and just look-up the value for the asteroid's current position (interpolating between grid points if you like).

Simon

theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #5 - Posted 2011-08-04 14:20:32 »

I'm using the same solution as Pitbuller. The only thing I have to do is getting rid of atan2. Sqrt should at least be faster. xd
Euler integration is fine with circular orbits, but the whole orbit starts to rotate when the orbit is elliptic. Look at the first test in this article: http://codeflow.org/entries/2010/aug/28/integration-by-example-euler-vs-verlet-vs-runge-kutta/ (run it by holding your mouse over it).

The main goal is to get stable orbits with good performance. Asteroids do NOT influence each other. I may however have multiple stars/planets affecting the body. From a performance perspective Euler integration is perfect, but the problem with high acceleration needs a solution. I'm going to experiment with using a smaller time step when the gravitational force gets higher. That would only increase the computation cost for a fraction of the time. Will report back!

EDIT: Adapting the time step to the force was a bad idea. The orbit gets unstable...
Does anybody have an implementation of RK4 or the means to teach me how to implement it myself? I would like to test it. xd

Myomyomyo.
dishmoth
 « Reply #6 - Posted 2011-08-04 18:40:54 »

It might be worth giving the semi-implicit Euler method a try.  It conserves (kinetic plus potential) energy, which might help with the stability of your orbits.
Simon

pjt33

« JGO Spiffy Duke »

Medals: 40
Projects: 4
Exp: 7 years

 « Reply #7 - Posted 2011-08-04 20:25:57 »

Does anybody have an implementation of RK4 or the means to teach me how to implement it myself? I would like to test it. xd
Here you go. I've included a simplified version of one calculation I was using it for as a test case.
 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99  100  101  102  103  104  105  106  107  108  109  110  111  112  113  114  115  116  117  118 `/**  * Represents a point in R^4. Doubles are used. Instances of this class are immutable. * @author Peter Taylor */public class Vec4{    private final double w, x, y, z;        /** Constructs a point in R^4 from its projections onto the four axes. */    public Vec4(double _w, double _x, double _y, double _z) { w = _w; x = _x; y = _y; z = _z; }        /** Vector addition. Adds components in each axis. */    public Vec4 add(Vec4 q) { return new Vec4(w + q.w, x + q.x, y + q.y, z + q.z); }    /** Scalar multiplication. */    public Vec4 scale(double scalar) { return new Vec4(w * scalar, x * scalar, y * scalar, z * scalar); }    /** Combines scalar multiplication and addition. */    public Vec4 addScaled(Vec4 q, double scalar) {        return new Vec4(w + q.w * scalar, x + q.x * scalar, y + q.y * scalar, z + q.z * scalar);    }        /** Gets projection onto w axis; this is the first component. */    public double getW() { return w; }    /** Gets projection onto x axis; this is the second component. */    public double getX() { return x; }    /** Gets projection onto y axis; this is the third component. */    public double getY() { return y; }    /** Gets projection onto z axis; this is the fourth component. */    public double getZ() { return z; }        /** Output as "(w, x, y, z)". */    public String toString() { return "(" + w + ", " + x + ", " + y + ", " + z + ")"; }    /** Function from (time, Vec4) to Vec4 */    public static interface Fn {        Vec4 eval(double t, Vec4 arg);    }    /** The most familiar of the fourth order Runge-Kutta methods for first order differential equations. */    public static class RK4 {        private final Fn fn;        public RK4(Fn _fn) { fn = _fn; }        /** Takes a step in the integration. This is defined by
* k1 = h f(t, X)         * k2 = h f(t + h/2, X + k1/2)         * k3 = h f(t + h/2, X + k2/2)         * k4 = h f(t + h, X + k3)         * X' = X + (k1 + 2k2 + 2k3 + k4)/6
* @param X The current position in R^4         * @param h The step size */        private Vec4 takeStep(Vec4 X, double t, double h) {            double h2 = h/2;            Vec4 k1 = fn.eval(t, X);            Vec4 k2 = fn.eval(t + h2, X.addScaled(k1, h2));            Vec4 k3 = fn.eval(t + h2, X.addScaled(k2, h2));            Vec4 k4 = fn.eval(t + h, X.add(k3));            Vec4 sum = k2.add(k3);            sum = k1.addScaled(sum, 2).add(k4);            return X.addScaled(sum, h / 6);        }        /** Solves the equation from t1 to t2 in steps of delta.         * t2 must be at least as great as t1.         * @param X The value of X(t1)         * @param t1 The start value of t         * @param t2 The end value of t         * @param delta The step size         * @return The computed value of X(t2), using dX/dt = fn(X). */         public Vec4 integrate(Vec4 X, double t1, double t2, double delta) {             double range = t2 - t1;             if (range < 0) throw new IllegalArgumentException(t2 + " is smaller than " + t1 + " - negative range");             // Do as much as can be done in steps of delta             double t;             for (t = t1; t < t2 ; t += delta) X = takeStep(X, t, delta);             if (t2 - t > 0.0) X = takeStep(X, t, t2 - t);             return X;         }    }    /** The differential equation corresponding to a trajectory with drag proportional to the square of the speed.     * We use the fields of the Vec4 as (w, x) for position and (y, z) for velocity. */    public static class Trajectory implements Fn {        private final static double g = 10.0;        private final double mass, surfaceArea, dragCoeff;        /** Generate a differential equation giving the trajectory of an object with         * the specified mass (from which we calculate a surface area for purposes         * of evaluating drag) and drag coefficient. */        public Trajectory(double _mass, double _dragCoeff) {            mass = _mass;            surfaceArea = Math.pow(_mass, 2.0 / 3);            dragCoeff = _dragCoeff;        }        /** We take the drag force to be proportional to the surface area (and hence         * to the mass raised to 2/3) times the velocity times the drag coefficient.         * This gives us: dw/dt = y; dx/dt = z; dy/dt = -surfaceArea * y * dragCoeff / mass;         * dz/dt = -surfaceArea * z * dragCoeff / mass - g */        public Vec4 eval(double t, Vec4 X) {           double w = X.getY();           double x = X.getZ();           double y = -surfaceArea * X.getY() * dragCoeff / mass;           double z = -surfaceArea * X.getZ() * dragCoeff / mass - g;           return new Vec4(w, x, y, z);        }    }    /** Test RK4 on Trajectory by checking a known trajectory. */    public static void main(String[] args) {       Trajectory traj = new Trajectory(1, 0.0);       Vec4 initialX = new Vec4(0.0, 0.0, 1.0, 10.0);       // Should take two seconds to hit ground again       Vec4.RK4 rk4 = new Vec4.RK4(traj);       Vec4 finalX = rk4.integrate(initialX, 0.0, 2.0, 0.01);       System.out.println("finalX = " + finalX);    }}`
theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #8 - Posted 2011-08-04 22:32:44 »

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

Myomyomyo.
dishmoth
 « Reply #9 - Posted 2011-08-05 08:42:49 »

Concerning Semi-implicit Euler: What is the difference between it and normal Euler?

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

theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #10 - Posted 2011-08-05 16:17:29 »

So semi-implicit Euler is just a fancy name for this code?
 1  2  3  4  5  6  7  8  9  10 `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?

Myomyomyo.
dishmoth
 « Reply #11 - Posted 2011-08-05 17:29:58 »

So semi-implicit Euler is just a fancy name for this code?

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

theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #12 - Posted 2011-08-05 20:08:24 »

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.

Myomyomyo.
theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #13 - Posted 2011-08-05 20:47:17 »

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.

Myomyomyo.
pjt33

« JGO Spiffy Duke »

Medals: 40
Projects: 4
Exp: 7 years

 « Reply #14 - Posted 2011-08-05 21:25:04 »

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
 1  2  3  4  5  6  7  8  9 `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);    }}`
theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #15 - Posted 2011-08-05 21:46:21 »

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:
 1 `double dx = SUN_X - pos.w, dy = SUN_Y - pos.x;`

What's up with pos.w and pos.x?

Myomyomyo.
dishmoth
 « Reply #16 - Posted 2011-08-06 14:39:39 »

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.

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.

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

pjt33

« JGO Spiffy Duke »

Medals: 40
Projects: 4
Exp: 7 years

 « Reply #17 - Posted 2011-08-06 15:05:20 »

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:
 1 `double dx = SUN_X - pos.w, dy = SUN_Y - pos.x;`

What's up with pos.w and pos.x?
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.

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.

...

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.
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).
dishmoth
 « Reply #18 - Posted 2011-08-06 18:25:42 »

As a point of interest, my FunOrb game Torquing uses RK4 on a Vec13.

It's not a proper plug if you don't give a link.

pjt33

« JGO Spiffy Duke »

Medals: 40
Projects: 4
Exp: 7 years

 « Reply #19 - Posted 2011-08-06 20:29:46 »

As a point of interest, my FunOrb game Torquing uses RK4 on a Vec13.

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
theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #20 - Posted 2011-08-07 06:15:13 »

@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:
 1  2  3  4  5 `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...

Myomyomyo.
theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #21 - Posted 2011-08-07 06:29:25 »

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

Myomyomyo.
pjt33

« JGO Spiffy Duke »

Medals: 40
Projects: 4
Exp: 7 years

 « Reply #22 - Posted 2011-08-07 06:44:50 »

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.
dishmoth
 « Reply #23 - Posted 2011-08-07 08:44:45 »

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!

theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #24 - Posted 2011-08-07 20:39:17 »

I optimized away all object creation in my code. I then did a small performance test on Semi-explicit Euler, Verlet and RK4. Results:
 1  2  3  4  5  6  7  8  9  10  11  12 `100.000 objects, 1 attractorSemi-implicit Euler:    4.5 msVerlet:         4.6 msRK4:         10.2 ms100.000 objects, 10 attractorsSemi-implicit Euler:    17.0 msVerlet:         17.0 msRK4:         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?
 1  2  3  4  5  6  7  8  9  10  11  12  13 `/*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

Myomyomyo.
delt0r

JGO Knight

Medals: 33
Exp: 18 years

Computers can do that?

 « Reply #25 - Posted 2011-08-10 11:05:13 »

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.

I have no special talents. I am only passionately curious.--Albert Einstein
theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #26 - Posted 2011-08-14 13:07:24 »

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

Myomyomyo.
delt0r

JGO Knight

Medals: 33
Exp: 18 years

Computers can do that?

 « Reply #27 - Posted 2011-08-15 12:09:35 »

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.

I have no special talents. I am only passionately curious.--Albert Einstein
theagentd

« JGO Bitwise Duke »

Medals: 429
Projects: 2
Exp: 8 years

 « Reply #28 - Posted 2011-08-15 15:59:06 »

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

Myomyomyo.
Pages: [1]
 ignore  |  Print

You cannot reply to this message, because it is very, very old.

 Riven (32 views) 2015-04-16 10:48:47 Duke0200 (43 views) 2015-04-16 01:59:01 Fairy Tailz (33 views) 2015-04-14 20:13:12 Riven (35 views) 2015-04-12 21:36:37 bus hotdog (50 views) 2015-04-10 02:39:32 CopyableCougar4 (52 views) 2015-04-10 00:51:04 BurntPizza (52 views) 2015-04-06 22:06:58 ags1 (54 views) 2015-04-02 10:58:48 Riven (53 views) 2015-04-01 18:27:05 ags1 (70 views) 2015-03-31 10:55:12
 theagentd 27x BurntPizza 16x wessles 15x 65K 11x kingroka123 11x Rayvolution 11x alwex 11x KevinWorkman 9x phu004 8x kevglass 8x Olo 7x ra4king 7x Ecumene 7x Roquen 7x chrislo27 7x Hanksha 7x
 How to: JGO Wikiby Mac702015-02-17 20:56:162D Dynamic Lighting2015-01-01 20:25:42How do I start Java Game Development?by gouessej2014-12-27 19:41:21Resources for WIP gamesby kpars2014-12-18 10:26:14Understanding relations between setOrigin, setScale and setPosition in libGdx2014-10-09 22:35:00Definite guide to supporting multiple device resolutions on Android (2014)2014-10-02 22:36:02List of Learning Resources2014-08-16 10:40:00List of Learning Resources2014-08-05 19:33:27
 java-gaming.org is not responsible for the content posted by its members, including references to external websites, and other references that may or may not have a relation with our primarily gaming and game production oriented community. inquiries and complaints can be sent via email to the info‑account of the company managing the website of java‑gaming.org