Simulation of gravity on large bodies

[poo2thegeek]

For my own personal project I want to simulate our solar system on a computer program.
I understand that the force applied on a body due to gravity can be calculated by dividing the product of the two masses by the square of the distance between them, then multiplying this value by the gravitational constant.
I assume that the direction this force is applied from would be a vector equal to the normalised value of the first body's position subtracted from the the second mass, with the direction being negated for the second mass.

But how would I then calculate the new position of the body.

For my example the sun can be static as the movement it undergoes is too little to bother calculating.
Mercury must have a velocity and a direction down which it is traveling. The suns gravity would attract the planet, but it's inertia would keep in from falling into the sun. I assume I must calculate the the angle tangent to the sun, but how would I go about doing this in 3 dimension?

 

[Simon Bridge]

Welcome to PF;
That's a pretty good project and you should learn a lot by playing with it. Our exact solar system is a very big project, you should practice just on two bodies to start with.

But you really just want a central system and ignoring the interaction between planets?
See: http://home.comcast.net/~szemengtan/ClassicalMechanics/SingleParticle.pdf
... s1.7 "Central Force Problem".

for a planet mass m (and Sun mass M) at position $\vec r$ the important equation is: $$\frac{d^2}{dt^2}\vec r = -\frac{GM}{r^3}\vec r$$ ... the force is entirely radial. Notice how the mass m of the planet does not matter here? This is the famous result that objects fall at the same rate.

The planet's velocity is $\vec v = \frac{d}{dt}\vec r$

The motion of the planet will be in the plane formed by the v and r vectors.
If you pick the plane in advance, it becomes a 2D problem.

For very small time increments - the new position is close to:
$\vec r(t+\delta t) = \vec r(t)+ \vec v(t) t + \frac{1}{2}\vec a(t) t^2$

But be warned - the system is very sensitive to rounding errors ... so even with the correct numbers your simulation is likely to go wrong quite fast. It's usually better to computer the entire orbit first, then sample it at a small time-interval to build an animation.

 

[poo2thegeek]

Welcome to PF;
But be warned - the system is very sensitive to rounding errors ... so even with the correct numbers your simulation is likely to go wrong quite fast. It's usually better to computer the entire orbit first, then sample it at a small time-interval to build an animation.

What if I were to later simulate an object that is effected by the gravitational force of say, the sun, but not to a strong enough extent that it fulfils an orbit?

Thank you you for the link and the information, I will look into it further :)

 

[TurtleMeister]

I used the information here to write my first n-body program. I have since had a lot of fun with it making modifications and doing various experiments. Although it was meant as a Java introductory programming assignment, the info can be used for any language.

Also, the assignment is 2d only for the xy plane. You might want to start out that way, but it's pretty easy to add the z axis.

 

[Simon Bridge]

What if I were to later simulate an object that is effected by the gravitational force of say, the sun, but not to a strong enough extent that it fulfils an orbit?

That just means you give it a higher kinetic energy than if you wanted a captured object.
These conditions are all covered in the pdf in the link - which also has an "introduction to classical mechanics".

And now you have a useful programmers link.
Between them you should have a lot of fun :)

Like TurtleMeister says, do 2d first since the third dimension is easy to add - especially if the planet-planet forces are zero. Gradually tweak it up as you get used to the dynamics.

 

[K^2]

A few words of warning. Gravity is an extremely difficult problem to tackle numerically if you desire precision. It does not take a whole lot of physics knowledge, but it requires a lot of numerical computing knowledge. You have a gravitational force, you compute acceleration, use that to update velocity, and finally, update position. This is called Euler Method. But numerically, you hit a snag almost right away. Gravity at new position is slightly different from initial position. Doesn't seem like much, but after enough time steps, you'll start accumulating an error. Typically, your object starts gaining velocity. In short term, this causes orbital precession that shouldn't be there, and in the long term, it causes your entire system to fly apart.

In itself, this is nothing out of the ordinary. Almost any numerical integration problem will have errors due to finite quantization. And it's relatively easy to improve on this using more advanced integration techniques. If Euler doesn't work, your next step would typically be to try Verlet Integration. If your interactions are due to harmonic potentials, like simple springs, for example, you'll get rid of almost all systematic errors and gain a huge boost in precision. Not with gravity, though. You can then go to higher order methods, like RK4. And again, you will see an improvement, but for a gravity problem, not enough of it.

Typically, you start getting decent results once you switch to implicit methods. Gauss-Legendre Method, for example, gives you pretty good precision for the chosen time step. But even that won't be enough if you try to simulate the Solary System. There are more advanced methods that I'm not sufficiently familiar with, that one would use to get sufficiently good results to actually make predictions for trajectories of asteroids and comets.That said, don't feel too discouraged. You won't get anything with scientific precision, but as far as playing with simple systems and learning what does and does not work, it's a pretty good problem. It's also a good introduction into integration methods. If you can program the very basic, Euler Method integration of trajectories, then you are entirely capable of modifying that code to run Verlet Integration. It's very simple, and you'll see how much difference it makes.

The site that TurtleMeister linked uses leapfrog instead of Verlet. They are similar. In general, Verlet is a little better. But in this particular case, you are unlikely to see much difference. Go for whichever makes more sense, or try both.

 

[AlephZero]

A few words of warning. Gravity is an extremely difficult problem to tackle numerically if you desire precision. It does not take a whole lot of physics knowledge, but it requires a lot of numerical computing knowledge.

That depends how you look at it. If you want accurate results, your numerical method has to be consistent with conservation of energy, and conservation of angular monentum. If the energy is wrong, your orbits will turn into spirals, either inwards or outwards. If angular momentum isn't conserved, you might get a repeating orbit, but the period will be wrong.

In other words, interpreting the numerical method in terms of physics, not just math, can be a good way to see what's good or bad about it.

The buzzword for all this is "symplectic" - but don't worry if Googliing that doesn't find anything you can understand. Just learn by doing!

 

[K^2]

As I pointed out, even a Verlet, which is a symplectic integrator, results in energy buildup over time. It's not as bad as naive explicit Euler, but it's going to make precise computations of trajectories in Solar system impossible. This gets even worse if you are trying to compute trajectory of a small body, such as asteroid or a comet, influenced by much larger bodies. From perspective of asteroid/comet, you are dealing with a time-dependent potential, and energy is not conserved to begin with.

If you'd like to give it a try, you can get positions of all major players in the Sol System from JPL's Horizons and try integrating a trajectory of some comet for which you can also find data. You'll find significant disagreement unless you go to a high order implicit method.

 

[poo2thegeek]

I used TutrleMeister's method most, as it was by far the easiest to implement within a computer program. The difference in acceleration and movement is calculated in by a BigDecimal, with a rounding of 128 decimal places. This should theoretically mean that the majority of any numerical errors would come from the rounding of the Gravitational Constant, which I have set as 6.67384E-11
If anyone has a more accurate/exact value for the constant then could they please give me a link to it?

Thanks everyone for all this information! I haven't done a large amount of this mathematics before as I am still quite young in my education, but I'll try my best to learn. Thanks again all!

 

[sophiecentaur]

It is very straightforward to write a routine that will give you a mass, orbiting around a much larger mass - to a level of accuracy that makes it 'look ok'. You just do this step by step, from one point to the next, adding in the effect of the gravitational force to the initial velocity. I did it years ago on a Psion 3 and it looked very convincing and was stable for hours (precessing a bit, eventually running away or collapsing but even that could be caught with a suitable test or better precision).

 

[poo2thegeek]

I've created this simulation in 2D and it works pretty well. But when I bring it into 3d I have a problem based on the speed. I have calculated at the maths based on the accurate values for the sun, and for mercury (just two planets for simplicity for now), and then when I draw the objects I just multiply the values by a rather small number to make the objects easy to draw on one screen.
However when I do this in 3D each frame seems to jerk by a rather large amount.
I assume I just need to add a multiplier into the algorithm at some point. But I've no idea where this multiplier would go.

BigDecimal delta_x = sun.x.subtract(x);
			BigDecimal delta_y = sun.y.subtract(y);
			
			BigDecimal distance = BigDecimalMath.sqrt(delta_x.pow(2).add(delta_y.pow(2)));
			
			BigDecimal force = (sun.mass.multiply(mass).multiply(Physics.BIG_G)).divide(distance.pow(2),MathContext.DECIMAL128);
			
			BigDecimal force_x = force.multiply(delta_x.divide(distance,MathContext.DECIMAL128));
			BigDecimal force_y = force.multiply(delta_y.divide(distance,MathContext.DECIMAL128));

			BigDecimal accelerationx = force_x.divide(mass,MathContext.DECIMAL128);
			BigDecimal accelerationy = force_y.divide(mass,MathContext.DECIMAL128);

			this.speedx = speedx.add(accelerationx);
			this.speedy = speedy.add(accelerationy);

			this.x = this.x.add(speedx);
			this.y = this.y.add(speedy);

I assumed I would multiply the 'speedx' and 'speedy' variables, but this seems to change the way the entity behaves without really changing it's jerkiness. I tried the same with then acceleration varaibles... I've no idea where else I could put the value in?

 

[D H]

I used TutrleMeister's method most, as it was by far the easiest to implement within a computer program. The difference in acceleration and movement is calculated in by a BigDecimal, with a rounding of 128 decimal places. This should theoretically mean that the majority of any numerical errors would come from the rounding of the Gravitational Constant, which I have set as 6.67384E-11
If anyone has a more accurate/exact value for the constant then could they please give me a link to it?

Thanks everyone for all this information! I haven't done a large amount of this mathematics before as I am still quite young in my education, but I'll try my best to learn. Thanks again all!

The method described in TurtleMeister's post claims to be using leapfrog, but it isn't. Look at the algorithm. It's just the symplectic Euler method, aka the Euler-Cromer method, aka a whole bunch of other names. The *only* thing that is worse is basic Euler, where you simply reverse the order in which position and velocity are updated compared to symplectic Euler.

There is absolutely no reason to use BigDecimals for Euler integration. It's extreme overkill. You might as well use floats (not even doubles). It's that lousy of a technique. It's important to understand how Euler's method works because it's at the base of many numerical integration techniques. Once you understand it, it's equally important to drop it for something better. Keep in mind that anything is better.

There also isn't much value in making your accelerations BigDecimals. If you do use BigDecimals, you need to use them for position and velocity. It's the additions x+vΔt and v+aΔt that kill the precision rather than the calculations of vΔt and aΔt. BTW, the code you posted in post #11 doesn't have a Δt. That factor is extremely important.

With regard to your question about G, don't use it. Don't use mass and don't use force either. Instead compute acceleration directly using $\ddot{\vec x} = -\,\displaystyle\frac{\mu_{\odot}}{||\vec x||^3} \vec x$, where mu ($\mu_{\odot}$) is the solar gravitational parameter, 132,712,440,018 km³/s² (or 1.32712440018×10²⁰ m³/s²). Conceptually, you can think of the standard gravitational parameter for some body as the product of G and the body's mass. That's not quite right, however. The mass of the Sun is calculated by dividing the solar gravitational parameter by G. The gravitational parameter is highly observable. For example, the uncertainty in $\mu_{\odot}$ is 9 km³/s². That's tiny compared to the uncertainty in G and the solar mass. You can see a list of somewhat up-to-date gravitational parameters for the Sun, all the planets, and some of the lesser bodies in the solar system on this wikipedia page: http://en.wikipedia.org/wiki/Standard_gravitational_parameter.

 

[tfr000]

A full-blown simulation of the Solar System is difficult. I have attempted it, with some success, but bits of it wobble away from where they should be - particularly the bits closer to the Sun.
If you want to pursue it, study numerical integration first. Search around for: Runge-Kutta, Cowell's method, Encke's method. Simulate something simple for starters - one object in an elliptical orbit around another, in 2D. Learn about orbits, and how to analyze them - you will need to know whether your two objects are acting like they are supposed to. When you are ready to go 3D multi-object, you will need to program a high-order integration method for accuracy. Search for: Gauss-Jackson (in particular, the Berry-Healy paper), Dormand-Prince, Runge-Kutta-Fehlberg. There are plenty of other methods, but these are some of the basic Multi-step and Single-step methods, and you can find them online.
For initial conditions, you will need to get data from JPL (their DE series) or IMCCE (their INPOP series). You will need to know some astronomy to understand these.
If you want to see a working simulation, look at SOLEX: http://chemistry.unina.it/~alvitagl/solex/