1
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ 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... $\endgroup$ – Brian Borchers Mar 29 '14 at 4:32
5
$\begingroup$

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

Example of constrained trajectory Total energy as function of time

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.