Interesting thing about the double pendulum - it is very difficult to get the numerical integration correct, even using very high order integrators like RK4. Th

Interesting thing about the double pendulum - it is very difficult to get the numerical integration correct, even using very high order integrators like RK4. The simple Euler integrator used here has no hope of getting it right (not a criticism, just an observation!) The problem is that the double pendulum equations (in fact, any equations of motion derived from a Hamiltonian) have symplectic structure, i.e. they have conserved quantities. One of the conserved quantities is the total energy of the system. When you discretize the equations to simulate them on a computer, you lose this conservation property (due to discretization error - nothing to do with floating point). The structure of the equations means that the total energy becomes an increasing quantity in time, so you tend to see the simulated system "speed up" or become more energetic as the simulation progresses. The error builds up in the same direction over time - exactly what you don't want! The simplest example is the harmonic oscillator, with second-order equation of motion x'' = -x which gives the first order equations in terms of position x and momentum p x' = p p' = -x which conserve the total energy 0.5 * (x^2 + p^2). Discretizing these using a first-order forwards Euler scheme x(t+1) = x(t) + p(t) * dt p(t+1) = p(t) - x(t) * dt you can see that the total energy changes on each time step to x(t+1)^2 + p(t+1)^2 = x(t)^2 + 2 x(t) p(t) dt + x(t)^2 dt^2 + p(t)^2 - 2 x(t) p(t) dt + p(t)^2 dt^2 = (1 + dt^2) (x(t)^2 + p(t)^2) so the total energy increases by a factor of (1 + dt^2) each step. Over time, the total energy increases exponentially. The solution is to use a geometric or sympletic integrator which explicitly takes into account the symplectic structure, producing a set of discrete update equations which still conserve a total energy quantity. [0] http://en.wikipedia.org/wiki/Symplectic_integrator http://en.wikipedia.org/wiki/Symplectic_integrator