Among other things, Christopher James Huff wrote:

> There is no way to compute an exact solution for 3 or more
> gravitationally interacting bodies. However, you can compute a very
> close approximation. I did a simulation (in C++) that was able to
> simulate the orbits of the planets from Mercury to Jupiter for several
> Jovian years (several dozen Earth years), with very little drift. (I
> simply haven't entered the data for the other planets)
> There are some rather complex algorithms needed to do this with both
> accuracy and efficiency, however. And there will always be cumulative
> errors, the simulation will diverge from reality over time, there is no
> way around that.

Exactly. In fact, such systems have a chaotic behaviour. By "chaotic" I 
don't mean they're undeterministic, but they show a strong dependence on 
initial conditions and divergence. If you don't mimic the *exact* initial 
conditions and don't solve *exactly* the movement equations, after some 
(maybe short) time, the solution you get is completely different from the 
real solution. It's the same thing that happens with weather forecast, 
although the equations look simpler.

>> Jelby, in your second post,  you offer an interesting idea.
>> >  p2 = p1 + v1*dt + 1/2*a1*(dt^2)
>> But is fundamentally any different from going to a smaller time slice? 
>> And
> 
> This basically assumes that acceleration will be constant over the time
> interval, and takes that into account for the new point. I take it you
> were originally using p2 = p1 + v1*dt?
> This does give a more accurate result, but is still not exact. You can
> get even better results by doing this one step, then backing up and
> using the average of the acceleration at the original position and the
> first approximation of the new position to compute the real new
> position, thus taking into account the fact that the acceleration
> changes over the time step. This isn't as good as you can get either,
> but it's easy to compute.

Yep, like all others this is only an approximation. It is designed to be 
cheap in computation power and storage. Position is computed using a simple 
second-order approximation. Velocity is computed using a first-order 
approximation, but taking the average of the acceleration. There are other 
algorithms, called "predictor-corrector", which basically do what you 
suggest.

Incidentally, the program I mentioned (Moldy) uses one of these algorithms. 
Even if it's written in C(++?), the algorithm is outlined in its manual, 
it's worth a look, I believe.

>> I was hoping for a solution where dt= one povray frame.
> 
> That is probably a bad idea. The accuracy of the simulation will then
> depend on the number of frames in the animation. You may get very
> different results when you do a longer render.

In the field of molecular dynamics, you usually don't store the data every 
time-step. To get a precise enough solution you need small dt, but then the 
configurations are too similar to be interesting. I'd aim for a frame every 
100*dt, more or less, and see what happens.

-- 
light_source{9+9*x,1}camera{orthographic look_at(1-y)/4angle 30location
9/4-z*4}light_source{-9*z,1}union{box{.9-z.1+x clipped_by{plane{2+y-4*x
0}}}box{z-y-.1.1+z}box{-.1.1+x}box{.1z-.1}pigment{rgb<.8.2,1>}}//Jellby

Post a reply to this message

In article <404178f6$1@news.povray.org>,
 "Greg M. Johnson" <gregj;-)565### [at] aolcom> wrote:

> Take a "sun".
> Take any mass and throw it in any direction not exactly at its center at any
> velocity.
> 
> If the mass doesn't sink into the sun,  will it:
> 
> i)  undertake an "orbit"  "right there" in the sense that it will undertake
> an elliptical orbit which intersects that point of release,

In a 2-mass system, just the sun and the orbiting mass, it will either 
do exactly that or escape the system entirely.


> OR
> ii) will a particle "destined" for a stable orbit sometimes "drift" out
> towards the radius of its stable orbit.

2 body systems are stable. Maybe you're mistakenly thinking of 
elliptical orbits as unstable?


> If I had  been much, much more patient, might some of those orbits which
> looked like they were decaying be just wandering off to something "close
> enough" to their own stable orbit?

If a 2-body system drifts significantly, it is your simulation which is 
unstable. In the real world, none of the planets in our solar system 
follow perfect elliptical orbits, and they do drift around and transfer 
energy to each other, but the level of instability is too small to 
matter. When the sun goes into its red giant stage and swallows up the 
inner planets, they will still be in very similar orbits to what they 
have today.

The p2 = p1 + v1*dt + 1/2*a1*(dt^2) equation will give a more accurate 
simulation than p2 = p1 + v1*dt, but for my simulation of the inner 5 
planets, I found it insufficient. The two-step method I described was 
sufficient, the simulation was stable for at least several decades 
without trying very hard.

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: <chr### [at] tagpovrayorg>
http://tag.povray.org/

Post a reply to this message

In article <cjameshuff-186E57.00003729022004@news.povray.org>,
 Christopher James Huff <cja### [at] earthlinknet> wrote:

> > So I'm falling back to three ideas:
> > 1) it is in fact impossible
> 
> Exact solutions are impossible. "Close enough" solutions are quite 
> possible.

BTW, here's a demo of several different integrators and orbit types:
http://www.princeton.edu/~rvdb/JAVA/astro/galaxy/Galaxy.html

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: <chr### [at] tagpovrayorg>
http://tag.povray.org/

Post a reply to this message

Christopher James Huff <cja### [at] earthlinknet> wrote:

> BTW, here's a demo of several different integrators and orbit types:
> http://www.princeton.edu/~rvdb/JAVA/astro/galaxy/Galaxy.html

Seems to be quite buggy though. Switching between integrators gives 
inconsistent results: sometimes earth escapes the sun with RungeKutta4
at warp=500 (I guess it is really using Euler then), sometimes it is
stable as expected.
Not very good conditions for comparing the integrators.

Lutz-Peter

Post a reply to this message

Okay great info thanks all.

I had heard that the gas giants were a Godsend in that they protected the
earth from wandering asteroids and wanted to model it. Hence my desire for
an actual working multi-body system.   Thanks.

Post a reply to this message

In article <4042a654$1@news.povray.org>,
 "Greg M. Johnson" <gregj;-)565### [at] aolcom> wrote:

> Okay great info thanks all.
> 
> I had heard that the gas giants were a Godsend in that they protected the
> earth from wandering asteroids and wanted to model it. Hence my desire for
> an actual working multi-body system.   Thanks.

That's one theory. Another theory is that Earth is so wet because of 
comets smacked into it by the gas giants.

They do pull in a great deal of debris that comes too close to them, but 
mostly stuff that would never go near Earth anyway. They also strongly 
deflect debris orbits within a much greater volume than the volume they 
"sweep clean", sometimes kicking it into the inner system.

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: <chr### [at] tagpovrayorg>
http://tag.povray.org/

Post a reply to this message

<"Greg M. Johnson" <gregj;-)565### [at] aolcom>> wrote:

> I had heard that the gas giants were a Godsend in that they protected the
> earth from wandering asteroids and wanted to model it. Hence my desire for
> an actual working multi-body system.

See my post "Euler vs Runge-Kutta4" in povray.binaries.images.
It seems this actually will solve the problem, if you're still interested.

Lutz-Peter

Post a reply to this message

In scientific simulations involving multiple particles, the energy is
achieved by specifying not the force between particles, but the energy
potential field around them. The force is then derived from the derivative
of this field. This can stop simulations going off to infinity, because one
can calculate the total energy of the system and then scale the velocities
of the particles appropriately so that the total energy is exactly what it
was at the start of the simulation.


proportional to -1/r. This is the formula to put into the energy calculation
when summing the potential energies.

Obviously there are still errors here - the "scaling" is just a cheat (has
no physical basis), but at least the error is spread through all the
particles, which have slightly different speeds, rather than in only one
particle which shoots off at a totally unrealistic speed.

-Chris

Post a reply to this message

In article <404269a7$1@news.povray.org>,
 Lutz-Peter Hooge <lpv### [at] gmxde> wrote:

> Christopher James Huff <cja### [at] earthlinknet> wrote:
> 
> > BTW, here's a demo of several different integrators and orbit types:
> > http://www.princeton.edu/~rvdb/JAVA/astro/galaxy/Galaxy.html
> 
> Seems to be quite buggy though. Switching between integrators gives 
> inconsistent results: sometimes earth escapes the sun with RungeKutta4
> at warp=500 (I guess it is really using Euler then), sometimes it is
> stable as expected.
> Not very good conditions for comparing the integrators.

True...I only tried a few examples before posting. It's slow 
too...slower than I think it should be, even for a Java program. I've 
done similar simulations with thousands of particles in Sapphire, which 
ought to be quite a bit slower than Java, and they ran quite smoothly.

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: <chr### [at] tagpovrayorg>
http://tag.povray.org/

Post a reply to this message

Greg M. Johnson 565### [at] aolcom> wrote:


> When I tried two years ago, I wasn't able to set up a correct model of the
> solar system, (perhaps it might help to say to do an "AI" model?),  without
> the planets either sinking into the sun or leaving the system altogether.  I
> think I tried pretty hard, and i could not find a gravitational constant (or
> a sufficiently small time slice) where particles didn't do one or the other.

The reason might be that you are numerically integrating a system of
differential equation using what is known as Euler's algorithm. The
central problem with that algorithm is that you assume the rate of
change to be constant during the time interval Delta t. Because you then
use x(t1) to calculate x(t2) (and x(t2) to calculate x(t3) and so on),
the initially small error becomes progressively larger.

One possible solution is to use Runge-Kutta's algorithm for the
integration. Instead of using dx(t1), this algorithm trys to approximate
an average rate of change. The simplest case is the 2nd order RK:

1) calculate dx(t1) and use it to calculate an approximation for x(t2)
2) calculate dx(t2)
3) take the arithmetic mean of dx(t1) and dx(t2) and use that to
calculate a final estimate of x(t2)

If you do a paper trace for, say, a parabula you will immediately see
the improvement that RK brings over Euler.

There are higher order RK, which minimise the error further. In
addition, you can use the difference between the results of an n-th and
an (n-1)-th order RK to estimate the error and controll the step width
Delta t to maintain an error smaller than a predefined value epsilon,
and at the same time reduce the computation time as much as possible.
But implementing this in PovRay is probably overkill. 

If you are interested and can read German, a good paper on this topic,
including an implementation of a 5(4)-th order RK in C(*) may be found
in 

@article{Hei-92,
        AUTHOR= {G. Heinzel},

von Differentialgleichungen},
        YEAR= {1992},
        JOURNAL= {c't},
        PAGES= {172-185},
        VOLUME= {8},
        NUMBER= {8},
        MONTH= {aug},
        ABSTRACT= {Schrittweitensteuerung durch Vergleich der Ergebnisse
von n-ter und n-1 ter Ordnung},
        LANGUAGE= {dt}}

(*) C: third letter of the roman alphabet. 1) Chemistry: Symbol for the
element number 12, carbon. 2) Physics: abbreviation for the unit of ->
charge (Coulomb) and symbol for -> capacity. 3) Technology: 
macroassembler for the (outdated) -> PDP11. Its deficiencies were the
cause of the so-called -> Software-crisis during the last millenium.
[;-)]

Post a reply to this message