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