# RATTLE numerical integrator example

I want to understand how the RATTLE algorithm works. Can somebody give me an example (in pseudocode or using any programming language like python or matlab) of how would I implement a numerical integrator for the simple pendulum problem?

• It might help readers of this question to know that RATTLE is a variation of Verlet's algorithm used in molecular dynamics simulations. See physics.ujep.cz/~mlisal/md/rattle.pdf There are example Verlet implementations all over the web... – Brian Borchers Mar 29 '14 at 4:32

The following snippet of code is an implementation of RATTLE on a system with the constraint $g(x, y) = K x^2 + y^2 - 1 = 0$.

// Constraint.
double g(const double2& r) {
return K * r.x * r.x + r.y * r.y - 1.0;
}

// Gradient of constraint.
double2 G(const double2& r) {
return double2(2.0 * K * r.x, 2.0 * r.y);
}

void rattle(double2& q, double2& p, double& lambda, double h) {
// Declare auxiliary constants.
const double2 Gqprev = G(q);

// Deal with constraint on the configuration manifold.
q += h * p;
double lambda_r = 0.0;
// Solve using Newton's method.
for (size_t k = 1; k <= max_iters; k++) {
if (k == max_iters)
abort();

const double2 r = q - lambda_r * Gqprev;
const double phi = g(r);
const double dphi_dl = -dot(G(r), Gqprev);
const double update = phi / dphi_dl;

if (fabs(phi) < tol && fabs(update) < tol)
break;

lambda_r -= update;
}

q -= lambda_r * Gqprev;
p -= lambda_r / h * Gqprev;

// Deal with constraint on the tangent space.
const double2 Gq = G(q);
double lambda_v = dot(Gq, p) / dot(Gq, Gq);
p -= lambda_v * Gq;

lambda = (lambda_r + lambda_v) / 2.0;

assert(fabs(g(q)) < tol && fabs(dot(G(q), p)) < tol);
}


I use the notation from chapter 7 in the book Simulating Hamiltonian Dynamics.

For your convenience, I have uploaded a standalone program that you can compile, run, and play with at https://gist.github.com/jmbr/9857614#file-rattle-cpp The program produces output that can be used to plot the graphs shown below (an example trajectory and its corresponding plot of total energy as a function of time).