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