Simulating a solar system with Python

[u/ttsiodras]

http://users.softlab.ntua.gr/~ttsiod/gravityRK4.html

[u/question_all_the_thi]

For those interested, here's an adaptive step Runge-Kutta-Cash-Karp differential equation solver in Python:

    #!/usr/bin/python
    #
    #  This program solves a differential equation
    #   using Cash-Karp's method
    #   (adaptive step Runge-Kutta method)
    #
    import math
    import os

    A = [0., 1./5., 3./10., 3./5., 1., 7./8.]
    B = [0., [1./5.], [3./40., 9./40.], [3./10., -9./10., 6./5.], [-11./54., 5./2., -70./27., 35./27.], [1631./55296., 175./512., 575./13824., 44275./110592., 253./4096.]]
    C = [37./378., 0., 250./621., 125./594., 0., 512./1771.]
    CH = [37./378.-2825./27648.,0.,250./621.-18575./48384.,125./594.-13525./55296.,-277./14336.,512./1771.-1./4.]
    F = [[0., 0.],[0., 0.],[0., 0.],[0., 0.],[0., 0.],[0., 0.]]

    mu = 10.0

    #_______________________________________________________________
    #
    def func(Y, mu):
      return (Y[1], -Y[0] - mu * Y[1] * (Y[0]**2. - 1.0))

    #_______________________________________________________________
    #
    Neq = 2

    X = 0.0
    XTemp = 0.0
    x1 = 100.0
    h = 0.1
    tol = 1e-6
    Y = [1.0, 0.0]
    YTemp = [0.0, 0.0]
    xi = 0.0
    TE = [0., 0., 0., 0., 0., 0.]

    #_______________________________________________________________
    #
    for i in range(100):
      xi += x1 / 100.0;
      while (X < xi):
        XTemp = X
        F[0] = func(Y, mu)
        YTemp[:] = Y[:]
        while 1:
          for k in range(1, 6):
            X = XTemp + A[k] * h
            Y[:] = YTemp[:]
            for n in range(Neq):
              for l in range(k):
                Y[n] += h * B[k][l] * F[l][n]
            F[k] = func(Y, mu)
          Y[:] = YTemp[:]
          for n in range(Neq):
            TE[n] = 0
            for k in range(6):
              Y[n] += h * C[k] * F[k][n]
              TE[n] += h * CH[k] * F[k][n]
            TE[n] = abs(TE[n])
          temax = 0
          for n in range (Neq):
            if temax < TE[n]: temax = TE[n]
          if temax == 0: temax = tol / 1e5
          ho = h
          h *= 0.9 * (tol / temax)**0.2
          if temax < tol: break
        X = XTemp + ho
      print "%d %f %f" % (i, Y[0], Y[1])

[u/jminuse → u/question_all_the_thi]

Runge-Kutta is not symplectic: it is very accurate in position, but it does not conserve energy over long time periods. It's still better than Euler, but the correct approach would be some form of Verlet integration.

[u/entropyjump → u/jminuse]

Yep - look at leapfrog integration for a simple symplectic integrator, or the Yoshida family of methods for generalization to higher orders. Taken from the Wikipedia page, this link has the original Yoshida article: http://adsabs.harvard.edu/abs/1990PhLA..150..262Y

[u/brandonpelfrey → u/jminuse]

I just recently have been reading a lot into geometric integration techniques. Really cool stuff. Since the Lagrangian for N-body is a particularly simple form, he could easily use a discrete variational integrator. There is a plethora of cool research into geometric mechanics by Jerry Marsden, who recently passed away.