Das Runge-Kutta-Verfahren und die Integration der
              ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
             Newtonschen Bewegungsgleichungen für N Massenpunkte
             ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀

      1) Das Runge-Kutta-Verfahren für explizite DGL-Systeme 1. Ordnung
         --------------------------------------------------------------

      Wir gehen von einem explizitem DGL-System 1. Ordnung aus mit N zu
      bestimmenden Funktionen xi(t), i=1,N, mit der einen unabhängigen
      Variablen t, die meist die Zeit bedeutet.  Gegeben sind dazu N
      Funktionen fi, i=1,N, die die ersten Zeitableitungen xi' in
      Abhängigkeit von den xi und der Zeit t beschreiben.  Ferner seien
      die Werte aller Funktionen xi zur Zeit t = to bekannt.  Die
      Aufgabe besteht nun darin, zum Zeitpunkt t1 = to + δt die
      Funktionswerte zu berechnen.  Vektoren stellen wir mit
      Großbuchstaben dar, sonst verwenden wir Kleinbuchstaben.  Wir
      haben also:

                 ┌     ┐                 ┌      ┐
                 │x1(t)│                 │x1'(t)│
                 │x2(t)│         dX      │x2'(t)│
          X(t) = │  ∙  │ ,  X' = ──(t) = │  ∙   │ ,
                 │  ∙  │         dt      │  ∙   │
                 │xN(t)│                 │xN'(t)│
                 └     ┘                 └      ┘

                        ┌       ┐
                        │f1(X;t)│
                        │f2(X;t)│
          X' = F(X;t) = │   ∙   │  ,  fi, i=1,..N, gegebene Funktionen,
                        │   ∙   │
                        │fN(X;t)│
                        └       ┘

      so daß sich das DGL-System auch kurz schreiben läßt in der Form:

          X' = F(X;t)  mit  Xo = X(to) ,

      mit Xo als der sog.  Anfangsbedingung.

      Zum Berechnen der Funktionswerte xi an der Stelle t1 := t0 + δt
      mit vorgegebener Schrittweite δt berechnet man nun nacheinander
      folgende Vektoren:

          KA = F( Xo          ; to       )
          KB = F( Xo + δt/2∙KA; to + δt/2)
          KC = F( Xo + δt/2∙KB; to + δt/2)
          KD = F( Xo +   δt∙KC; to + δt  )  ,

          K  = 1/6∙(KA + 2∙(KB + KC) + KD)  ,

          X(t1) = X(to + δt) = X(to) + K∙δt .

      Damit hat man die Funktionswerte an der Stelle t1 bestimmt.  Von
      t1 ausgehend kann man nun die Funktionswerte an den Stellen t2,
      ... usw.  berechnen.  Das ist das Runge-Kutta-Verfahren.


      2) Die Integration der Bewegungsgleichungen von N Massenpunkten
         ------------------------------------------------------------

      Seien Xi, i=1,N die Ortsvektoren der Massenpunkte, ihre Massen
      seien mi, i=1,N.  Der Beschleunigungsvektor Ai des i-ten
      Massenpunktes ist dann nach Newton (G:  Gravitationskonstante):

                   N       Xj(t) - Xi(t)
          Ai(t) =  Σ G∙mj∙────────────────   ,
                  j=1                    3
                  j╪i     │Xj(t) - Xi(t)│

      mit der Abkürzung

                       1
          Dij = ───────────────  ,
                │Xj(t) - Xi(t)│

      kompakter geschrieben:

                   N         3
          Ai(t) =  Σ G∙mj∙Dij ∙(Xj(t) - Xi(t)) .
                  j=1
                  j╪i

      Ai ist also nur von den Ortsvektoren Xj, j=1,N, abhängig.  Wir
      haben es mit einem DGL-System 2. Ordnung zu tun.  Für das
      Runge-Kutta-Verfahren, das für DGL-Systeme 1. Ordnung ist, kann
      man ein System 2. Ordnung in ein System 1. Ordnung (mit doppelt so
      vielen Gleichungen) nach folgendem Schema umschreiben.  Es wird
      gesetzt:

                          │               ┌     ┐          ┌     ┐
          X1' = V1        │               │X1(t)│          │V1(t)│
              ∙           │               │  ∙  │          │  ∙  │
              ∙           │               │  ∙  │          │  ∙  │
              ∙           │               │  ∙  │          │  ∙  │
          XN' = VN        │ bzw.:  Y(t) = │XN(t)│, Y'(t) = │VN(t)│,
          V1' = X1" ≡ A1  │               │V1(t)│          │A1(t)│
              ∙           │               │  ∙  │          │  ∙  │
              ∙           │               │  ∙  │          │  ∙  │
              ∙           │               │  ∙  │          │  ∙  │
          vN' = xN" ≡ AN  │               │VN(t)│          │AN(t)│
                          │               └     ┘          └     ┘

      so daß wir nun ein DGL-System 1. Ordnung bezüglich des Vektors Y
      erhalten, das nach Runge-Kutta gelöst werden kann:

                dY
           Y' = ── = F(Y)  ,
                dt

      dabei kommt also in diesem Fall die Zeit als unabhängige Variable
      nicht explizit auf der rechten Seite vor, was das Schema des
      Runge-Kutta-Verfahrens geringfügig vereinfacht.

      Man beachte, daß die Ai nur von den Ortsvektoren Xj, j=1,N
      abhängen.  Als Anfangsvektor Y0 gibt man die Ortsvektoren Xi und
      die Geschwindigkeiten Vi vor.

      Geben wir zur Zeit to die Xi, Vi und damit Y(to) vor, so können
      wir nach dem Runge-Kutta-Verfahren für den Zeitpunkt t1 = to + δt
      offenbar Y(t1) berechnen.  Für N Massenpunkte haben wir ein
      6N-dimensionales lineares DGL-System.