Help with symplectic integrators

I'm trying to develop a physics simulation and I want to implement a fourth-order symplectic integration method. The problem is that I must be getting the math wrong, since my simulation is not working at all when using the symplectic integrator (as compared to a fourth-order Runge-Kutta i开发者_运维技巧ntegrator that works reasonably well for the simulation). I've been googling this forever and all I can find are scientific articles on the subject. I tried to adapt the method used in the articles, but am having no luck. I want to know if anyone has source code for a simulation that makes use of symplectic integrators, preferably to simulate a gravitational field but any symplectic integrator would do. What language the source is in doesn't matter too much, but I would appreciate a language using C-style syntax. Thanks!

As you asked for source code: From HERE you can download MATLAB and FORTRAN code for symplectic methods for Hamiltonian systems and symmetric methods for reversible problems. And a lot of other methods for dealing with diff equations too.

And in THIS paper you can find the description of the algorithms.

Edit

If you use Mathematica this may help too.

Here is the source code for a fourth order composition method based on the Verlet scheme. A linear regression of $\log_{10}(\Delta t)$ vs. $\log_{10}(Error)$ will show a slope of 4, as expected from theory (see the graph below). More information can be found here, here or here.

#include <cmath>
#include <iostream>

using namespace std;

const double total_time = 5e3;

// Parameters for the potential.
const double sigma = 1.0;
const double sigma6 = pow(sigma, 6.0);
const double epsilon = 1.0;
const double four_epsilon = 4.0 * epsilon;

// Constants used in the composition method.
const double alpha = 1.0 / (2.0 - cbrt(2.0));
const double beta = 1.0 - 2.0 * alpha;


static double force(double q, double& potential);

static void verlet(double dt,
                   double& q, double& p,
                   double& force, double& potential);

static void composition_method(double dt,
                               double& q, double& p,
                               double& f, double& potential);


int main() {
  const double q0 = 1.5, p0 = 0.1;
  double potential;
  const double f0 = force(q0, potential);
  const double total_energy_exact = p0 * p0 / 2.0 + potential;

  for (double dt = 1e-2; dt <= 5e-2; dt *= 1.125) {
    const long steps = long(total_time / dt);

    double q = q0, p = p0, f = f0;
    double total_energy_average = total_energy_exact;

    for (long step = 1; step <= steps; ++step) {
      composition_method(dt, q, p, f, potential);
      const double total_energy = p * p / 2.0 + potential;
      total_energy_average += total_energy;
    }

    total_energy_average /= double(steps);

    const double err = fabs(total_energy_exact - total_energy_average);
    cout << log10(dt) << "\t"
         << log10(err) << endl;
  }

  return 0;
}

double force(double q, double& potential) {
  const double r2 = q * q;
  const double r6 = r2 * r2 * r2;
  const double factor6  = sigma6 / r6;
  const double factor12 = factor6 * factor6;

  potential = four_epsilon * (factor12 - factor6);
  return -four_epsilon * (6.0 * factor6 - 12.0 * factor12) / r2 * q;
}

void verlet(double dt,
            double& q, double& p,
            double& f, double& potential) {
  p += dt / 2.0 * f;
  q += dt * p;
  f = force(q, potential);
  p += dt / 2.0 * f;
}

void composition_method(double dt,
                        double& q, double& p,
                        double& f, double& potential) {
  verlet(alpha * dt, q, p, f, potential);
  verlet(beta * dt, q, p, f, potential);
  verlet(alpha * dt, q, p, f, potential);
}

I am in the field of accelerator physics (synchrotron light sources), and in modelling electrons moving through magnetic fields, we use symplectic integrators on a regular basis. Our basic workhorse is a 4th order symplectic integrator. As noted in the comments above, unfortunately these codes are not so well standardized or easilly available.

One open source Matlab based tracking code is called Accelerator Toolbox. There is a Sourceforge project called atcollab. See a messy wiki here https://sourceforge.net/apps/mediawiki/atcollab/index.php?title=Main_Page

For the integrators, you can look here: https://sourceforge.net/p/atcollab/code-0/235/tree/trunk/atintegrators/ The integrators are written in C with MEX link to Matlab. Because the electrons are relativistic the kinetic and potential terms look a little different than in the non-relativistic case, but one can still write the Hamiltonian as H=H1+H2 where H1 is a drift and H2 is a kick (say from quadrupole magnets or other magnetic fields).