Help with symplectic integrators

Question

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 integrator 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!

Thelmathem

Basically, I just want to integrate an N-body problem. I suppose the parameters then are the bodies' positions, momenta, and masses. – Thelmathem

I went under the assumption that general n-body problems cannot be solved symbolically, that it's the reason why numerical integrators (such as RK4 and symplectic integrators) are used. If you mean setting up the problem with the appropriate differential equations, don't worry about it. It took me a while to get the RK4 integrator working as well, but it has more to do with my ability to transfer mathematical equations to code (ie. it's possible to be able to do it, but also not be able to write code to do it). – Thelmathem

I blush. I read you question all to quickly and simply saw "symbolic" where you wrote "symplectic". My apologies, but all my comments (now deleted as to the wrong point) were based on this misapprehension. – Tripalmitin

Oh, ok, no worries then. The problem I'm having is just getting a symplectic integrator to code. I'm getting nowhere with my attempts so I just wanted to see how other people do it. – Thelmathem

Answer 1

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.

Creosol

Thanks, that definitely helps! What I needed was having the equations in papers described in code, and the link you provided does that. – Thelmathem

@Thelmathem As you saw, source code examples for symplectic algorithms are scarce on the web. When you are done, consider posting your code somewhere to help others in need. – Creosol

I can definitely appreciate how scare examples are. Although a large number of papers are written on symplectic algorithms used to solve various problems, it seems that those same scientists don't put the code for their algorithms on the web. I will post a link later if I can get a working symplectic algorithm. – Thelmathem

Answer 2

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);
}

Consuetudinary

Answer 3

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).

Melissiamelita