Building a symplectic numerical integrator

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
In a recent post, 'meson800' asked:

How does one go about building the numerical integrator? I understand how to numerically integrate using Simpson's rule and friends, so I imagine it would go something like the following:
Calculate total acceleration due to gravity using gravity equations and locations of ship+planets
Use Simpson's rule to calculate dV for a certain dT
Using that dV, calculate a dX(and dY and dZ) for that dT
Do it again!
Looking at the Runge-Kutta methods, it looks an awful lot like Simpson's rule applied to differential equations. Does the higher-order Runge-Kutta simply integrate using a cubic/quintic/etc function analogue instead of the quadratic Simpson's uses?"

In response I offered to write up some more detailed notes on numerical integrators and how to build them. This post constitutes these notes. In these notes, I'm going to concentrate on symplectic integrators

A second-order symplectic integrator
First, let's focus on a basic second-order symplectic integrator. It is worthwhile spending a little time on this because, so long as the time-step of the integration is small enough, it actually works pretty well. It is also an integrator from which it is possible to build higher-order syplectic integrators.

A couple of introductory points:
  1. The integrator is called 'symplectic' because (up to second order at least) it preserves the symplectic character of the physical system. This means that the long-run properties of the integrator tend to be much more stable than other kinds of integrators and, largely because of this, they have become very popular amongst physicists in the last couple of decades.
  2. The integrator is second-order because it has been constructed so as to 'kill' error terms up to the second-order only.

So, what is this integrator? If we just focus on a one-dimensional system for a moment, then the integrator maps a pair of points [MATH]\left\{Q_0,P_0\right\}[/MATH] to a new pair of points [MATH]\left\{Q_2,P_2\right\}[/MATH] at a time [MATH]\delta t [/MATH] later. We can think of this as taking a pair of numbers, [MATH]\left\{Q_0,P_0\right\}[/MATH], that describe the state of some object at some start time [MATH]t = t_0[/MATH] and updating this to a new pair of numbers, [MATH]\left\{Q_2,P_2\right\}[/MATH], that describes the new state of the same object at time [MATH]t = t_0 + \delta t[/MATH].

There are three steps to this updating rule:
[MATH]Q_{1} \leftarrow Q_{0}+\frac{1}{2}\,\delta t\, P_{0}[/MATH][MATH]P_{2} \leftarrow P_{0}+\delta t\, F\left(Q_{1}\right)[/MATH][MATH]Q_{2} \leftarrow Q_{1}+\frac{1}{2}\,\delta t\, P_{2}[/MATH]
OK, so what do all of these symbols mean? Let's start with 'Q' and 'P'. In physics, largely because of a longstanding convention in nomenclature, 'Q' is often used to denote a spatial coordinate of something. Here, in our one-dimensional example, you can think of 'Q' as representing the x-coordinate of some object moving in some gravitational field. In the same convention, 'P' is used to to denote the momentum of the same object. (Formally, it is the generalised momentum conjugate to 'Q' but that's a nuance that we don't need to worry about here.). Here, we can think of 'P' as representing the x-coordinate of the momentum of a particle. But since momentum is just mass * velocity (in this coordinate system), we can think of 'P' as the x-coordinate of velocity (multiplied by the mass of the object). So, the pair of numbers [MATH]\left\{Q_0,P_0\right\}[/MATH] just represents the position and velocity of the object at some initial time. And this is just the state-vector of the object written in cartesian coordinates. In other words, the symplectic integrator is an updating rule that takes an object's state vector at some initial time, and returns a new state vector at some time [MATH]\delta t [/MATH] later. ( Now, you can extend this idea to three dimensions but for the time being we'll just stick with one-dimension.)

Next, let's focus on the 'F' term. Again by convention, 'F' is used to denote the force acting on an object. Specifically, [MATH]F\left(Q_{1}\right)[/MATH] represents the force on the object being modelled when it is at position [MATH]Q_1[/MATH]. So, to carry out the three steps of the integration updating rule, we need to provide some force function, [MATH]F(Q)[/MATH] which allows us to calculate the force on an object at point in our one-dimensional space (i.e., at any 'Q'). For object moving subject to a gravitational force from a single body, we know that:

[MATH] F(Q) = -\frac{\mu\, m}{(Q-Q^*)^{2}}[/MATH]
(provided that [MATH]Q > Q^*[/MATH])

where [MATH]\mu[/MATH] is the gravitational constant for the body in question; 'm' is the mass of the object that we are modelling, and [MATH]Q^*[/MATH] is the location of the source of the gravitational field (e.g., the centre of the Sun or the Earth). If work in Gaussian units where we measure distances in AU (Astronomical Units); velocity in AU/day; and mass in units of the mass of the Sun, then for an object moving in the gravitational potential of the Sun, [MATH]\mu = 0.00029591220828559115[/MATH] and [MATH]m = 1[/MATH]. For an object moving in the gravitational field of the Earth, [MATH]m = 1/354710[/MATH]. To convert from AU and days to metres and seconds, one needs to know that 1 AU = 149597870700 metres; and that 1 day = 86400 seconds. Unless the mass of the object that we are modelling is changing, it is convenient to set [MATH]m = 1[/MATH].

In this one-dimensional system, the object moving in a gravitational potential is a bit limited in terms of directional options. It can either go up, or it can go down. Clearly, as the object moves towards the location of the gravitational source, the force acting on the object is going to increase without limit and the integration is going to going hay-wire, but so long as we are a reasonable distance away from [MATH]Q^*[/MATH] then the integration scheme given above will provide a reasonable description of the object's motion.

In three dimensions
In three dimensions, the integration scheme looks much the same. But now we have to apply it to three spatial coordinates and three velocity components - i.e., the integration scheme becomes one that updates one set of 6 numbers [MATH]\left\{Q_{x,0},Q_{y,0},Q_{z,0},P_{x,0},P_{y,0},P_{z,0}\right\}[/MATH] to a new set of six numbers [MATH]\left\{Q_{x,2},Q_{y,2},Q_{z,2},P_{x,2},P_{y,2},P_{z,2}\right\}[/MATH]. And the second order integration scheme that does this is as follows:

[MATH]Q_{x,1} \leftarrow Q_{x,0}+\frac{1}{2}\,\delta t\, P_{x,0}[/MATH][MATH]Q_{y,1} \leftarrow Q_{y,0}+\frac{1}{2}\,\delta t\, P_{y,0}[/MATH][MATH]Q_{z,1} \leftarrow Q_{z,0}+\frac{1}{2}\,\delta t\, P_{z,0}[/MATH]
[MATH]P_{x,2} \leftarrow P_{x,0}+\delta t\, F_x\left(Q_{x,1},Q_{y,1},Q_{z,1}\right)[/MATH][MATH]P_{y,2} \leftarrow P_{y,0}+\delta t\, F_y\left(Q_{x,1},Q_{y,1},Q_{z,1}\right)[/MATH][MATH]P_{z,2} \leftarrow P_{z,0}+\delta t\, F_z\left(Q_{x,1},Q_{y,1},Q_{z,1}\right)[/MATH]
[MATH]Q_{x,2} \leftarrow Q_{x,1}+\frac{1}{2}\,\delta t\, P_{x,2}[/MATH][MATH]Q_{y,2} \leftarrow Q_{y,1}+\frac{1}{2}\,\delta t\, P_{y,2}[/MATH][MATH]Q_{z,2} \leftarrow Q_{z,1}+\frac{1}{2}\,\delta t\, P_{z,2}[/MATH]
Each step of the integrator is the same - except now we apply it three times - once for each spatial coordinate.

You should not that we now have three force functions rather than just one. This is because force is really a 'vector' rather than 'scalar' quantity. [MATH]F_x\left(Q_{x},Q_{y},Q_{z}\right)[/MATH] is the force acting in the x-direction on the object located at position [MATH]{Q_x, Q_y, Q_z}[/MATH]; [MATH]F_y\left(Q_{x},Q_{y},Q_{z}\right)[/MATH] is the force acting in the y-direction on the object located at position [MATH]{Q_x, Q_y, Q_z}[/MATH]; and [MATH]F_z\left(Q_{x},Q_{y},Q_{z}\right)[/MATH] is the force in the x-direction acting on the object located at position [MATH]{Q_x, Q_y, Q_z}[/MATH]. And if we are in the gravitational field of a single body, then these functions become:

[MATH]F_x\left(Q_{x},Q_{y},Q_{z}\right) = - \frac{\mu\,m\,\left(Q_{x}-Q_x^*\right)}{\left(\left(Q_{x}-Q_x^*\right)^2+\left(Q_{y}-Q_y^*\right)^2+\left(Q_{z}-Q_z^*\right)^2\right)^{3/2}}[/MATH]
[MATH]F_y\left(Q_{x},Q_{y},Q_{z}\right) = - \frac{\mu\,m\,\left(Q_{y}-Q_y^*\right)}{\left(\left(Q_{x}-Q_x^*\right)^2+\left(Q_{y}-Q_y^*\right)^2+\left(Q_{z}-Q_z^*\right)^2\right)^{3/2}}[/MATH]
[MATH]F_z\left(Q_{x},Q_{y},Q_{z}\right) = - \frac{\mu\,m\,\left(Q_{z}-Q_z^*\right)}{\left(\left(Q_{x}-Q_x^*\right)^2+\left(Q_{y}-Q_y^*\right)^2+\left(Q_{z}-Q_z^*\right)^2\right)^{3/2}}[/MATH]
where [MATH]\left(Q_x^*,Q_y^*,Q_z^*\right)[/MATH] are the coordinates of the gravitating body (when the object whose motion you are modelling is at [MATH]\left(Q_x,Q_y,Q_z\right)[/MATH])


And that's about all there is to this this second-order symplectic integrator. To integrate, one chooses a suitable small step-size (so that the overall errors are as small as some tolerance that you require for your calculations); one specifies the initial conditions of the object - its position and velocity - and then one repeatedly applies the integration step for as long as you wish.

Now, for those interested, it is worthwhile setting up this integrator in, say, an spreadsheet and seeing how it performs under various sizes of time-steps. and initial conditions. Of course, if there is more than one gravitating body, then you have additional terms in the force functions, but the basic scheme of the updating rule remains the same.

As a sequel to this post, I will sketch how you can quickly build fourth and sixth order symplectic integrators from this simple second-order integrator.

And as a sequel to that I will talk about the distinction between 'implicit' and 'explicit' integrators.
 
Last edited:

meson800

Addon Developer
Addon Developer
Donator
Joined
Aug 6, 2011
Messages
405
Reaction score
2
Points
18
Awesome post!

Here's some questions I have, to see if I understand.


The integrator is second-order because it has been constructed so as to 'kill' error terms up to the second-order only.

Ok, so now I see that, if I just do a "naive" integrator, errors will accumulate in position.

Assuming one-dimensional, with the simulated point accelerating towards some gravity source.

For example, If the momentum is updated as in the second-order system
[MATH]P_{1} \leftarrow P_{0}+\delta t\, F\left(Q_{0}\right)[/MATH]
and then do a naive update based on the initial momentum
[MATH]Q_1 \leftarrow Q_0 + \delta t P_0[/math]then, over time, the position value will be underestimated.

On the other hand, if the position is based on the final momentum
[math]Q_1 \leftarrow Q_0 + \delta t P_1[/math]then the position value will be overestimated.​

Therefore, the second-order updater, by doing the position update in two steps,
[MATH]Q_{1} \leftarrow Q_{0}+\frac{1}{2}\,\delta t\, P_{0}[/MATH]
[MATH]Q_{2} \leftarrow Q_{1}+\frac{1}{2}\,\delta t\, P_{2}[/MATH]
has the overestimate and underestimate cancel out. Neat!

However, doesn't the momentum of the object get under/overestimated over time with this scheme, like the position in the naive approach? After a while, wouldn't that accumulate as position, and therefore energy, errors?


As a sequel to this post, I will sketch how you can quickly build fourth and sixth order symplectic integrators from this simple second-order integrator

But I assume that this ^^ is what fixes the momentum error.
 
Last edited:

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
Hi, 'meson800'

Although it may be possible to view symplectic integrators in terms of 'predictor-corrector' methods, I think the reason why symplectic integration often proves to be so stable is due to a somewhat different line of reasoning.

To explain this reasoning, let's just focus on the standard Kepler problem of a (small) satellite in orbit around some other massive body (e.g., the Earth). The stability of symplectic integration applied to the Kepler problem is actually quite stunning: even after integrating for thousands of orbits, or even hundreds of thousands of orbits, there is very little increase or decrease in the orbital energy. Although the energy doesn't remain absolutely constant, the energy of the system is bound so that it tends to oscillate up and down with each orbit but never steadily increases or decreases. This stability of symplectic integration is due entirely to the fact that these integrators preserve the symplectic character of the underlying system. But why this is important for long-run stability takes a bit of explaining

First, let's state that the standard Kepler problem is an example of an 'integrable system'. Basically, to be 'integrable' means you can always find coordinates in which orbital motion is just uniform motion around a circle. These coordinates are generally known as 'action-angle' coordinates. For the Kepler problem the relevant 'action-angle' coordinates are total angular momentum (the 'action' variable), and the mean anomaly (the 'angle' variable). For the Kepler problem the full set of action-angle variables (of which there are six) are the Delaunay variables. No need to worry too much about the full set here - just focusing on the two (angular momentum and mean anomaly) will suffice here. So, imagine a sheet of paper with a whole series of concentric circles inscribed on it. Simplifying the problem a bit, any closed orbit can be represented as uniform motion around one of these circles. And no matter where you start on the piece of paper, for an integrable system, you will always find yourself following uniform motion on one of these concentric circles. The generic name for these circles is 'tori' - and each one is a 'torus'.

OK, so far this has nothing to do with symplectic integrators. But here is the chain of logic that makes the connection. (I warn you, though, what follows is mind bending stuff.)

There are two ways of viewing a symplectic integrator. There is the conventional way - that is, to see it as an approximate time evolution operator for the true physical system - in this case the standard Kepler problem. The less conventional way is to see the symplectic integrator as the exact time evolution operator of some perturbed Kepler problem. Now, we may not know exactly what that perturbation is, but we know it exists. And if we could integrate that system exactly, then we would get the exactly the right prediction for our perturbed system using our symplectic integrator. Now, if the integrator is not symplectic, then you can't view the integrator in this way so the rest of the logic doesn't apply. What comes next is only valid if you are using a symplectic integrator.

This alternative view of symplectic integrators - namely that they are the exact time evolution operator of a perturbed system is important because of a famous piece of theory called 'KAM theory'. 'KAM' stands for 'Komogorov, Arnol'd and Moser' - three people who did some ground-breaking work in this area. What KAM theory says if you take an 'integrable' physical system and then add a perturbation to it, not all of the closed circles (the tori) survive - but an appreciable number do. They are bent out of a shape a bit by the perturbation, but it is still true that if you start on one of those tori, then you will always move around in a circle on it. As the perturbation gets smaller - e.g., if the 'delta_t' of our symplectic integration gets smaller - then more of the tori survive following the addition of the perturbation. But what happens between the surviving tori? Well, basically motion becomes chaotic. So, KAM theory says that after you add a perturbation to an integrable system, extreme order will be replaced with chaos - interspersed with islands of stability - the surviving tori.

So, for the Kepler problem, a symplectic integrator effectively adds a (hopefully small) perturbation to the system. Generally, for an arbitrary starting position one is going to find oneself on one of the chaotic orbits and so, at first blush, you might think that long-run stability might be compromised. But if you get go back to the picture of the concentric circles on a piece of paper, we know that some of the concentric circles survive and that for a small perturbation (e.g., small delta_t of our integrator) those concentric circles, the tori, that survive the addition of the perturbation are usually very close together. Let's suppose, then, that you start of one of the seas of chaos between the surviving tori. Then motion is bound above and below. The dynamics of the system do not allow an orbit to jump across a surviving tori so that even though though the long-run evolution is chaotic, the orbit that a symplectic integrator would predict is going to be tightly restricted. And it is this tight constraint on chaotic motion in physical systems that gives rise to the long-run stability of symplectic integrators. Curious, but true.

Now, this picture of concentric circles is a pretty reasonable one for simple systems like the Kepler problem. In this two-dimensional system the principle of confinement of seas of chaos is pretty easy to visualise. But for systems with more degrees of freedom, this simple picture breaks down. It is no longer necessarily true that seas of chaos are bound by the surviving tori. Generally, there are 'bolt holes' that connect chaotic regions. However, so long as the perturbation (i.e., delta_t) is small enough, the chance that an integrated trajectory will go through one of these 'bolt-holes' tends to be exponentially small - so that still, for very long times, the trajectory calculated using a symplectic integrator remains tightly bound.

This way of thinking about things also has implications for why fourth and higher order symplectic integrators tend to be better than, say, a second-order integrator: with a fourth order integrator, the perturbation of the integrable system tends to be smaller, which means that more of the tori survive, which means the chaotic spaces between the surviving tori are even closer, which means that the integrated trajectory is even more closely bound.

Having said all of this, the plain fact of the matter is that symplectic integrators are just more stable than non-symplectic integrators. And so, unless you have a particularly good reason for not using one, you should use a symplectic integrator.
 
Last edited:

MikeB

Member
Joined
Feb 25, 2009
Messages
185
Reaction score
0
Points
16
Location
Seattle
This is very interesting, but will take some work to digest. Can you recommend a reference for further reading? Maybe one with pictures? ;)
 

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
I will see what I can find.
 
Last edited:

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
On 14 April, 'MikeB' wrote:

This is very interesting, but will take some work to digest. Can you recommend a reference for further reading? Maybe one with pictures?

Hi, 'MikeB'. I've scoured the web for some for of introductory reference on this. Frankly, I'm struggling a bit to find something that isn't ridiculously technical. The trouble is, I think, that KAM theory in general - and its specific application to symplectic integrators - isn't something that has been well captured by popular science media. Or at least if it is, I haven't found it so far.

There are, of course, reams of technical papers on KAM theory which, if you suffer from insomnia and have an acutely mathematical bent, I am happy to forward. But an introductory article (with pictures) - nada!

Perhaps the best non-technical text I've found is a book called "The KAM Story" (see http://www.amazon.com/The-KAM-Story-Introduction-Kolmogorov-Arnold-Moser/dp/9814556580) - but even in Kindle form, it isn't cheap. Having said that, the picture on the front cover pretty much captures the essence of the story.

Given the shortage of readily available stuff, it may be best if I put together a few words focusing on the simple pendulum showing how tori are destroyed if one uses a symplectic integrator - but showing how even though tori are destroyed, integrated paths are forced to remain forever bound between neighbouring surviving tori.
 
Last edited:

MikeB

Member
Joined
Feb 25, 2009
Messages
185
Reaction score
0
Points
16
Location
Seattle
Perhaps the best non-technical text I've found is a book called "The KAM Story" (see http://www.amazon.com/The-KAM-Story-Introduction-Kolmogorov-Arnold-Moser/dp/9814556580) - but even in Kindle form, it isn't cheap. Having said that, the picture on the front cover pretty much captures the essence of the story.

The price is a bit steep for my level of interest, but it's interesting to use Amazon's "look inside". And the cover picture is suggestive of the stability/chaos regions you mentioned. Thanks for the tip, and I look forward to your further exploration of this topic for us amateurs.

That reminds me that this topic might fit in with the idea of Rocket Science For Amateurs. Would you consider developing an article for the OrbiterWiki site?
 

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
Hi, MikeB

Happy to write something for OrbiterWiki. However, I'm a bit busy at the moment and don't have much 'bandwidth' left over.

If one were to write something on symplectic integrators, one should probably say something about the uses to which they can be put as well as how to build them.
 
Last edited:

MikeB

Member
Joined
Feb 25, 2009
Messages
185
Reaction score
0
Points
16
Location
Seattle
Happy to write something for OrbiterWiki. ... If one were to write something on symplectic integrators, one should probably say something about the uses to which they can be put as well as how to build them.

The great thing about the wiki approach is that an article can start small, and be expanded as its authors/editors find the time and energy. If you don't mind, I can take what you've provided here, and post it as an article. Then as your bandwidth permits, you might go in and edit to correct my mistakes, or add info about uses, etc. Keep in mind we're not writing a textbook, just trying to give amateur rocket scientists some background to inspire their own exploration.
 

Topper

Addon Developer
Addon Developer
Donator
Joined
Mar 28, 2008
Messages
666
Reaction score
20
Points
33
For those who are intrested:
Here are some code shreds for some different integrators in Java:
Comments are as far as I understood without guarentee :)

These are the different methods I investigated:
1. Euler Method: Simple, not precise but fast (Not Symplectic?)
But it seems there is another "Euler" - Method which is: http://en.wikipedia.org/wiki/Semi-implicit_Euler_method (Sorry I had no time to study this)
2. Runge Kutta 4: Good precision (Not Symplectic?)
3. Verlet: Not as precise as RK4, but realistic "Look and feel") ==> Keeps alsways the total energy of the system (Symplectic?)

4. Physics: In the physics, you can see that it's optimized a bit, because it calculates the acceleration for both body in one step, so the distance etc. has only to be calculated once. For sure, there are more thing to optimize...

5. Collition detection: Off-Topic, bit I'm so happy with it ;-)

I used var names as they are used in the Orbiter API so I hope VECTOR3 and GBody are self explanatory.
But one important point: The VECTOR3 members of GBody have 4 fields for r=pos and v=velocity for x,y,z just for RK4 and Verlet.

Sorry I didn't got the differences regarding what is "Symplectic" and "not Symplectic" but I will study this point soon...
1:
Code:
package gl_sceneProject.src.sceneProject;

import java.util.ArrayList;

public class CalculatorEuler 
{
    ArrayList<GBody> gbodies =  new ArrayList<GBody>();
    ArrayList<GBody> ignoredGbodies = new ArrayList<GBody>();
    double updateEuler(double h1, int dt, ArrayList<GBody> gbodies2)
    {
        
        this.gbodies.clear();
        this.gbodies.addAll(gbodies2);
        this.ignoredGbodies.clear();
        for (GBody gb : this.gbodies)
            if (gb.mergedObject != null)
                ignoredGbodies.add(gb);
        
        this.gbodies.removeAll(ignoredGbodies);
        
        for (int step=0; step<dt; step++)
        {        
            TrackPlotter.storeTracksDeleted(ignoredGbodies);
            for (GBody star : gbodies)
            {
                star.r[0].add(star.v[0].mul(h1));
                TrackPlotter.storeTrack(star);
            }
            for (int i = 0; i < (gbodies.size()); i++)
            {                
                GBody body1 = gbodies.get(i);                
                for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
                {
                    GBody body2 = gbodies.get(i2);                
                    if (!body1.equals(body2))
                    {                    
                        updateEulerStep(body1, body2, h1);    
                    }                
                }            
            }
            ArrayList<GBody> collidet = CollisionDetection.check(gbodies);
            ignoredGbodies.addAll(collidet);
            gbodies.removeAll(collidet);
        }
        return h1*dt;
    }    
    
    private void updateEulerStep(GBody body1, GBody body2, double dt )
    {        
        double distance = body1.r[0].minus(body2.r[0]).length();
        VECTOR3 forceVector = Physic.getForceVector(body1.r[0], body2.r[0]);
        VECTOR3 force = Physic.getForce(forceVector, distance, body1.m, body2.m);                
        body1.a[0] = (Physic.getAcceleration(force.neg(), body1.m).mul(dt));        
        body2.a[0] = (Physic.getAcceleration(force, body2.m).mul(dt));   
        body1.v[0].add(body1.a[0]);
        body2.v[0].add(body2.a[0]);
    }

}
2: RK4
Code:
package gl_sceneProject.src.sceneProject;

import java.util.ArrayList;

public class CalculatorRungeKuttaV3 
{
    ArrayList<GBody> gbodies =  new ArrayList<GBody>();
    ArrayList<GBody> ignoredGbodies = new ArrayList<GBody>();
    public double update(double h1, int dt, ArrayList<GBody> gbodies2)
    {    
        this.gbodies.clear();
        this.gbodies.addAll(gbodies2);
        this.ignoredGbodies.clear();
        
        for (GBody gb : this.gbodies)
            if (gb.mergedObject != null)
                ignoredGbodies.add(gb);
        
        this.gbodies.removeAll(ignoredGbodies);
        double h2 = h1 * 0.5;
        for (int d=0; d<dt; d++) 
        {
            TrackPlotter.storeTracksDeleted(ignoredGbodies);
            resetA();
            int k;
            for (k=0; k<3; k++)
            {
                updateAccelerationV3(k);
                for (GBody body1 : gbodies) 
                {
                    updateK(body1.r[0], body1.v, body1.r, k, (k==3) ? h1:h2);    
                    updateK(body1.v[0], body1.a, body1.v, k, (k==3) ? h1:h2);
                }
            }
            updateAccelerationV3(k);
            for (GBody body1 : gbodies)
            {
                TrackPlotter.storeTrack(body1);
                body1.r[0].add (rk4mw(body1.v,h1));    
                body1.v[0].add (rk4mw(body1.a,h1));
            }
            ArrayList<GBody> collidet = CollisionDetection.check(gbodies);
            ignoredGbodies.addAll(collidet);
            gbodies.removeAll(collidet);
        }
        return h1*dt;
    }
    
    private void resetA()
    {
        for (GBody body :gbodies)
            for (int i=0; i<4; i++)
                body.a[i].reset();
    }
    
    private void updateAccelerationV3(int k)
    {
        for (int i = 0; i < (gbodies.size()); i++)
        {
            GBody body1 = gbodies.get(i);
            for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
            {
                GBody body2 = gbodies.get(i2);
                Physic.updateAccelerationRK3(body1, body2, k);
            }
        }
    }    
    
    private void updateK(VECTOR3 pos0, VECTOR3[] v0, VECTOR3[] pos, int k, double dt)
    {    
        pos[k+1].x = pos0.x + v0[k].x * dt;
        pos[k+1].y = pos0.y + v0[k].y * dt;
        pos[k+1].z = pos0.z + v0[k].z * dt;        
    }
    
    private VECTOR3 rk4mw(VECTOR3[] k, double dt)
    {
        double dt6 = dt / 6;
        return new VECTOR3
        (
            (k[0].x + k[1].x*2 + k[2].x*2 + k[3].x) * dt6,
            (k[0].y + k[1].y*2 + k[2].y*2 + k[3].y) * dt6,
            (k[0].z + k[1].z*2 + k[2].z*2 + k[3].z) * dt6
        );
    }
}
3 Verlet:
Code:
package gl_sceneProject.src.sceneProject;

import java.util.ArrayList;

public class CalculatorVerlet 
{
    ArrayList<GBody> gbodies = new ArrayList<GBody>();
    ArrayList<GBody> ignoredGbodies = new ArrayList<GBody>();
    boolean initialized = false;
    private void init(double h, int dt, ArrayList<GBody> gbodies)
    {
        updateAccelerationV3(0);
        for (GBody body : gbodies) 
        {
            body.r[1] = body.r[0].plus(body.v[0].mul(h)).plus(body.a[0].mul(0.5).mul(-h*h));
        }        
        initialized = true;
    }
    
    public void unInit()
    {
        //Estimate v;
        initialized = false;
    }
    
    double updateVerlet(double h, int dt, ArrayList<GBody> gbodies2)
    {        
        
        this.gbodies.clear();
        this.gbodies.addAll(gbodies2);
        this.ignoredGbodies.clear();
        
        for (GBody gb : this.gbodies)
            if (gb.mergedObject != null)
                ignoredGbodies.add(gb);
        
        this.gbodies.removeAll(ignoredGbodies);

        
        if (!initialized) init(h, dt, gbodies);
        int k = 2;
        
        for (int i = 0; i < dt; i++)    
        {
            TrackPlotter.storeTracksDeleted(ignoredGbodies);
            resetA();
            updateAccelerationV3(1);
            for (GBody body : gbodies) 
            {
                TrackPlotter.storeTrack(body);
                body.r[k] = (body.r[k-1].mul(2)).minus(body.r[k-2].plus(body.a[k-1].mul(-h*h)));
                body.r[0] = body.r[1];
                body.r[1] = body.r[2];
            }
            ArrayList<GBody> collidet = CollisionDetection.checkVerlet(gbodies,h);
            ignoredGbodies.addAll(collidet);
            gbodies.removeAll(collidet);
        }
        return h*dt;
    }    

    private void updateAccelerationV3(int k)
    {
        for (int i = 0; i < (gbodies.size()); i++)
        {
            GBody body1 = gbodies.get(i);
            for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
            {
                GBody body2 = gbodies.get(i2);
                Physic.updateAccelerationRK3(body1, body2, k);
            }
        }
    }
    
    private void resetA()
    {
        for (GBody body :gbodies)
            body.a[1].reset();
    }
}
4 Physics:
Code:
package gl_sceneProject.src.sceneProject;

public class Physic 
{
    static VECTOR3 getForce(VECTOR3 virtualPos1, double m1, GBody body2)
    {
        VECTOR3 er = body2.r[0].minus(virtualPos1);
        double distance = er.length();
        er = er.normalize();    
        er.multhis(SystemParameters.GRAV_CONST * m1 * body2.m  / (Math.pow(distance, 2)));
        return er;
    }
    
    static VECTOR3 getAcceleration(VECTOR3 f, double m)
    {
        return f.divi(m);
    }
    
    static VECTOR3 getForce(VECTOR3 er, double distance, double m1, double m2)
    {
        return er.mul(SystemParameters.GRAV_CONST * m1 * m2 / (Math.pow(distance, 2)));  //2a. f1
    }
    
    static VECTOR3 getForceVector(VECTOR3 pos1, VECTOR3 pos2)
    {    
        return new VECTOR3(pos1.x - pos2.x, pos1.y - pos2.y, pos1.z - pos2.z).normalize();
    }
    
    static void updateAccelerationRK3(GBody bodyFrom, GBody bodyTo, int k)
    {
        VECTOR3 V3_r = bodyTo.r[k].minus(bodyFrom.r[k]);
        double distance = V3_r.length();
        if (distance != 0)
        {
            VECTOR3 helper = V3_r.mul(SystemParameters.GRAV_CONST / (distance*distance*distance));
            bodyFrom.a[k].add( helper.mul(  bodyTo.m ));
            bodyTo.a[k].add( helper.mul( -bodyFrom.m ));
        }
    }    
}
5. Collition detection
Code:
package gl_sceneProject.src.sceneProject;

import java.util.ArrayList;

public class CollisionDetection 
{
	static public ArrayList<GBody> checkVerlet(ArrayList<GBody> gbodies, double h)
	{
		ArrayList<GBody> deleted = new ArrayList<GBody>();
		for (int i = 0; i < (gbodies.size()); i++)
		{				
			GBody body1 = gbodies.get(i);	        
			if (body1.mergedObject != null) continue;
			for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
			{				
				GBody body2 = gbodies.get(i2);
				if (body2.mergedObject != null) continue;
				GBody body3 = collisionDetectionVerlet(body1, body2, h);
				if (body3 != null) deleted.add(body3);
			}
		}
		return deleted;
	}
	
	static public ArrayList<GBody> check(ArrayList<GBody> gbodies)
	{
		ArrayList<GBody> deleted = new ArrayList<GBody>();
		for (int i = 0; i < (gbodies.size()); i++)
		{				
			GBody body1 = gbodies.get(i);	        
			if (body1.mergedObject != null) continue;
			for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
			{				
				GBody body2 = gbodies.get(i2);
				if (body2.mergedObject != null) continue;
				GBody body3 = collisionDetection(body1, body2);
				if (body3 != null) deleted.add(body3);
			}
		}
		return deleted;
	}
	
	static private GBody collisionDetection (GBody body1, GBody body2)
	{
        if(body1.r[0].distance(body2.r[0]) <= body1.rad + body2.rad)
        {
        	if (body1.m < body2.m)
        	{
	    		body2.v[0] = body1.v[0].mul(body1.m).plus(body2.v[0].mul(body2.m)).divi(body1.m + body2.m);
	    		body2.r[0] = body1.r[0].mul(body1.m).plus(body2.r[0].mul(body2.m)).divi(body1.m + body2.m);
	    		body2.setMass(body1.m + body2.m);
	    		body2.setDence((body1.getCalculatedDence()+body2.getCalculatedDence())/2);
	    		body1.setMass(0);
	    		body1.mergedObject = body2.getName();
	    		return body1;
        	}
        	else
        	{
	    		body1.v[0] = body2.v[0].mul(body2.m).plus(body1.v[0].mul(body1.m)).divi(body2.m + body1.m);
	    		body1.r[0] = body2.r[0].mul(body2.m).plus(body1.r[0].mul(body1.m)).divi(body2.m + body1.m);
	    		body1.setMass(body2.m + body1.m);
	    		body1.setDence((body2.getCalculatedDence()+body1.getCalculatedDence())/2);
	    		body2.setMass(0);
	    		body2.mergedObject = body1.getName();
	    		return body2;        		
        	}        	
        }		
        else return null;
	}
	static private GBody collisionDetectionVerlet (GBody body1, GBody body2, double h)
	{
        if(body1.r[0].distance(body2.r[0]) <= body1.rad + body2.rad)
        {
        	if (body1.m < body2.m)
        	{
	        	body1.v[0] = (body1.r[1].minus(body1.r[0])).mul(1/h);	//Estimate V
	        	body2.v[0] = (body2.r[1].minus(body2.r[0])).mul(1/h);	//Estimate V
	    		body2.v[0] = body1.v[0].mul(body1.m).plus(body2.v[0].mul(body2.m)).divi(body1.m + body2.m);
	    		body2.r[0] = body1.r[0].mul(body1.m).plus(body2.r[0].mul(body2.m)).divi(body1.m + body2.m);
	    		body2.setMass(body1.m + body2.m);
	    		body2.setDence((body1.getCalculatedDence()+body2.getCalculatedDence())/2);
	    		body1.setMass(0);
	    		body1.mergedObject = body2.getName();
	    		body2.r[1] = body2.r[0].plus(body2.v[0].mul(h)).plus(body2.a[0].mul(0.5).mul(-h*h));
	    		body2.r[2] = (body2.r[1].mul(2)).minus(body2.r[0].plus(body2.a[1].mul(-h*h)));
	    		//TrackPlotter.storeTrack(body1);
	    		//TrackPlotter.storeTrack(body2);
	    		return body1;
        	}
        	else
        	{
	        	body2.v[0] = (body2.r[1].minus(body2.r[0])).mul(1/h);	//Estimate V
	        	body1.v[0] = (body1.r[1].minus(body1.r[0])).mul(1/h);	//Estimate V
	    		body1.v[0] = body2.v[0].mul(body2.m).plus(body1.v[0].mul(body1.m)).divi(body2.m + body1.m);
	    		body1.r[0] = body2.r[0].mul(body2.m).plus(body1.r[0].mul(body1.m)).divi(body2.m + body1.m);
	    		body1.setMass(body2.m + body1.m);
	    		body1.setDence((body2.getCalculatedDence()+body1.getCalculatedDence())/2);
	    		body2.setMass(0);
	    		body2.mergedObject = body1.getName();
	    		body1.r[1] = body1.r[0].plus(body1.v[0].mul(h)).plus(body1.a[0].mul(0.5).mul(-h*h));
	    		body1.r[2] = (body1.r[1].mul(2)).minus(body1.r[0].plus(body1.a[1].mul(-h*h)));
	    		return body2;        	
	    	}
        }		
        else return null;
	}
}
I will public the whole project sources if it's finished.
 
Last edited:

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
Hi, 'MikeB'

Feel free to re-post on the Wiki.

Cheers.
 
Last edited:
Top