double-pendulum.svg

<svg height="100%" width="100%" viewBox="-2 -2 4 4" xmlns="http://www.w3.org/2000/svg">
  <title>Purely algebraic simulation of a double pendulum</title>
  <line x1="0" y1="0" stroke="grey" stroke-width="0.05" />
  <line stroke="grey" stroke-width="0.05" />
  <circle r="0.1" fill="black" />
  <circle r="0.1" fill="black" />
  <script>
function rk4(dt, x, f) {
  "use strict";
  let from = Array.from,
      a = from(f(from(x,  $    => $         )), $ => $*dt),
      b = from(f(from(x, ($,i) => $ + a[i]/2)), $ => $*dt),
      c = from(f(from(x, ($,i) => $ + b[i]/2)), $ => $*dt),
      d = from(f(from(x, ($,i) => $ + c[i]  )), $ => $*dt);
  return from(x, (_,i) => (a[i] + 2*b[i] + 2*c[i] + d[i])/6);
}

function setPendulumsPos(a,b,c,d) {
  const strings = document.getElementsByTagName("line"),
        balls = document.getElementsByTagName("circle");
  strings[0].setAttribute("x2", a);
  strings[0].setAttribute("y2", -b);
  strings[1].setAttribute("x1", a);
  strings[1].setAttribute("y1", -b);
  strings[1].setAttribute("x2", a+c);
  strings[1].setAttribute("y2", -(b+d));

  balls[0].setAttribute("cx", a);
  balls[0].setAttribute("cy", -b);
  balls[1].setAttribute("cx", a+c);
  balls[1].setAttribute("cy", -(b+d));
}

var abcd = [0,1,0,1],
    q    = [0.01,0,0,0];

setPendulumsPos(...abcd);


var previousTimestamp;
(function animate(timestamp) {
  if (previousTimestamp !== undefined) {
    let dq = rk4(
      (timestamp - previousTimestamp)/1000,
      q,
      ddot(...abcd)
    );
    q = Array.from(q, ($,i) => $ + dq[i]);
    [abcd, q] = psi(...abcd, ...q);
    setPendulumsPos(...abcd);
  }
  previousTimestamp = timestamp;
  window.requestAnimationFrame(animate);    
})();

function phi(a,b,l,ldot) {
  return [
    [-(2*b*l)/(1 + l**2) + (a*(1 - l**2))/(1 + l**2), (2*a*l)/(1 + l**2) + (b*(1 - l**2))/(1 + l**2)],
    [(4*b*l**2*ldot)/(1 + l**2)**2 - (2*a*l*(1 - l**2)*ldot)/(1 + l**2)**2 - (2*b*ldot)/(1 + l**2) - (2*a*l*ldot)/(1 + l**2), -(4*a*l**2*ldot)/(1 + l**2)**2 - (2*b*l*(1 - l**2)*ldot)/(1 + l**2)**2 + (2*a*ldot)/(1 + l**2) - (2*b*l*ldot)/(1 + l**2)]
  ]
}

function psi(a,b,c,d,lambda,mu,lambdadot,mudot) {
  let $phi = [
    phi(a,b,lambda,lambdadot),
    phi(c,d,mu,mudot)
  ],
      $a = $phi[0][0][0],
      $b = $phi[0][0][1],
      $c = $phi[1][0][0],
      $d = $phi[1][0][1],
      v1x = $phi[0][1][0],
      v1y = $phi[0][1][1],
      v2x = $phi[1][1][0],
      v2y = $phi[1][1][1];

  return [
    [$a, $b, $c, $d], 
    [
      0, 0,
      (-$b*v1x + $a*v1y)/(2*($a**2+$b**2)),
      (-$d*v2x + $c*v2y)/(2*($d**2+$c**2))
    ]
  ];
}

function ddot(a,b,c,d) {
  return ([l,m,lp,mp]) => [
    lp,
    mp,
     -(((1 + l**2)*(1 + m**2)*((2*a*d*(l - m + l**2*m - l*m**2) + 2*b*c*(m - l**2*m + l*(-1 + m**2)) + a*c*(1 + 4*l*m - m**2 + l**2*(-1 + m**2)) + b*d*(1 + 4*l*m - m**2 + l**2*(-1 + m**2)))*(2*d*(1 + l**2)**3*m*(1 + m**2) + c*(1 + l**2)**3*(-1 + m**4) - 4*b*c*lp**2*(1 + m**2)*(-1 - 6*l*m + 2*l**3*m + m**2 - 3*l**2*(-1 + m**2)) - 4*lp**2*(1 + m**2)*(b*d*(2*m - 6*l**2*m + 3*l*(-1 + m**2) - l**3*(-1 + m**2)) + a*(-(d*(-1 - 6*l*m + 2*l**3*m + m**2 - 3*l**2*(-1 + m**2))) + c*(2*m - 6*l**2*m + 3*l*(-1 + m**2) - l**3*(-1 + m**2)))) + 4*m*mp**2 + 12*l**2*m*mp**2 + 12*l**4*m*mp**2 + 4*l**6*m*mp**2) - 2*(1 + l**2)**2*(4*l*lp**2*(1 + m**2)**3 + 2*b*(1 + l**2)*(-((c - 3*c*m**2 + d*m*(-3 + m**2))*mp**2) + l**2*(c - 3*c*m**2 + d*m*(-3 + m**2))*mp**2 + l*(1 + 3*m**4 + m**6 - 2*d*mp**2 - 6*c*m*mp**2 + 2*c*m**3*mp**2 + m**2*(3 + 6*d*mp**2))) + a*(1 + l**2)*(-1 - 3*m**4 - m**6 + 2*d*mp**2 + 6*c*m*mp**2 - 2*c*m**3*mp**2 - 4*l*(c - 3*c*m**2 + d*m*(-3 + m**2))*mp**2 - 3*m**2*(1 + 2*d*mp**2) + l**2*(1 + 3*m**4 + m**6 - 2*d*mp**2 - 6*c*m*mp**2 + 2*c*m**3*mp**2 + m**2*(3 + 6*d*mp**2))))))/(4*(1 + l**2)**4*(1 + m**2)**4 - 2*(1 + l**2)**2*(1 + m**2)**2*(2*a*d*(l - m + l**2*m - l*m**2) + 2*b*c*(m - l**2*m + l*(-1 + m**2)) + a*c*(1 + 4*l*m - m**2 + l**2*(-1 + m**2)) + b*d*(1 + 4*l*m - m**2 + l**2*(-1 + m**2)))**2)),
    (8*b*d*lp**2*(1 + m**2)**3*(m - l**2*m + l*(-1 + m**2)) + 2*c**2*(1 + l**2)*(2*b**2*(m*(-1 + 3*m**2) + l**4*m*(-1 + 3*m**2) + l**3*(-1 + 10*m**2 - 5*m**4) + 2*l**2*m*(4 - 7*m**2 + m**4) + l*(1 - 10*m**2 + 5*m**4)) + a**2*(-2*l**2*m*(7 - 16*m**2 + m**4) + m*(3 - 4*m**2 + m**4) + l**4*m*(3 - 4*m**2 + m**4) - 2*l*(1 - 10*m**2 + 5*m**4) + 2*l**3*(1 - 10*m**2 + 5*m**4)) - a*b*(1 - 10*m**2 + 5*m**4 + 4*l*m*(5 - 10*m**2 + m**4) - 4*l**3*m*(5 - 10*m**2 + m**4) - 6*l**2*(1 - 10*m**2 + 5*m**4) + l**4*(1 - 10*m**2 + 5*m**4)))*mp**2 - 2*(1 + l**2)**3*m*(1 + m**2)**2*(d + d*m**2 + 2*mp**2) + 2*b**2*d*(1 + l**2)*(1 + 4*l*m - m**2 + l**2*(-1 + m**2))*(-(d*m*(-3 + m**2)*mp**2) + d*l**2*m*(-3 + m**2)*mp**2 + l*(1 + 3*m**4 + m**6 - 2*d*mp**2 + m**2*(3 + 6*d*mp**2))) + 2*a**2*d*(1 + l**2)*(l - m + l**2*m - l*m**2)*(-1 - 3*m**4 - m**6 + 2*d*mp**2 - 4*d*l*m*(-3 + m**2)*mp**2 - 3*m**2*(1 + 2*d*mp**2) + l**2*(1 + 3*m**4 + m**6 - 2*d*mp**2 + m**2*(3 + 6*d*mp**2))) + a*d*(4*lp**2*(1 + m**2)**3*(1 + 4*l*m - m**2 + l**2*(-1 + m**2)) + b*(1 + l**2)*(-1 + 2*m**6 + m**8 + 2*d*mp**2 + 10*d*m**4*mp**2 - 2*m**2*(1 + 10*d*mp**2) - 6*l**2*(-1 + 2*m**6 + m**8 + 2*d*mp**2 + 10*d*m**4*mp**2 - 2*m**2*(1 + 10*d*mp**2)) + l**4*(-1 + 2*m**6 + m**8 + 2*d*mp**2 + 10*d*m**4*mp**2 - 2*m**2*(1 + 10*d*mp**2)) - 8*l*(m + m**7 - 5*d*m*mp**2 + m**5*(3 - d*mp**2) + m**3*(3 + 10*d*mp**2)) + 8*l**3*(m + m**7 - 5*d*m*mp**2 + m**5*(3 - d*mp**2) + m**3*(3 + 10*d*mp**2)))) + c*(1 - 4*b*lp**2 + 2*m**2 - 8*b*lp**2*m**2 - 2*m**6 + 8*b*lp**2*m**6 - m**8 + 4*b*lp**2*m**8 - 2*b**2*d*mp**2 + 20*b**2*d*m**2*mp**2 - 10*b**2*d*m**4*mp**2 - 4*b**2*l**5*m*(1 + m**6 - 10*d*mp**2 + m**4*(3 - 2*d*mp**2) + m**2*(3 + 20*d*mp**2)) - l**6*(-1 + 2*m**6 + m**8 + 2*b**2*d*mp**2 + 10*b**2*d*m**4*mp**2 - 2*m**2*(1 + 10*b**2*d*mp**2)) + 4*b*l*m*(-4*lp**2*(1 + m**2)**3 + b*(1 + m**6 - 10*d*mp**2 + m**4*(3 - 2*d*mp**2) + m**2*(3 + 20*d*mp**2))) + l**4*(-3*(-1 + m**2)*(1 + m**2)**3 + 2*b**2*(-2 + 4*m**6 + 2*m**8 + 5*d*mp**2 + 25*d*m**4*mp**2 - 2*m**2*(2 + 25*d*mp**2))) + l**2*(-3*(-1 + m**2)*(1 + m**2)**3 - 4*b*lp**2*(-1 + m**2)*(1 + m**2)**3 + 2*b**2*(-2 + 4*m**6 + 2*m**8 + 5*d*mp**2 + 25*d*m**4*mp**2 - 2*m**2*(2 + 25*d*mp**2))) + a**2*(1 + l**2)*(-1 + 2*m**6 + m**8 + 2*d*mp**2 + 10*d*m**4*mp**2 - 2*m**2*(1 + 10*d*mp**2) + l**4*(-1 + 2*m**6 + m**8 + 2*d*mp**2 + 10*d*m**4*mp**2 - 2*m**2*(1 + 10*d*mp**2)) - 4*l*(m + m**7 - 10*d*m*mp**2 + m**5*(3 - 2*d*mp**2) + m**3*(3 + 20*d*mp**2)) + 4*l**3*(m + m**7 - 10*d*m*mp**2 + m**5*(3 - 2*d*mp**2) + m**3*(3 + 20*d*mp**2)) - 2*l**2*(-1 + 2*m**6 + m**8 + 6*d*mp**2 + 30*d*m**4*mp**2 - 2*m**2*(1 + 30*d*mp**2))) - 2*a*(4*lp**2*(1 + m**2)**3*(l - m + l**2*m - l*m**2) + b*(1 + l**2)*(m + m**7 - 10*d*m*mp**2 + m**5*(3 - 2*d*mp**2) + m**3*(3 + 20*d*mp**2) + 2*l*(-1 + 2*m**6 + m**8 + 4*d*mp**2 + 20*d*m**4*mp**2 - 2*m**2*(1 + 20*d*mp**2)) - 2*l**3*(-1 + 2*m**6 + m**8 + 4*d*mp**2 + 20*d*m**4*mp**2 - 2*m**2*(1 + 20*d*mp**2)) - 6*l**2*(m + m**7 - 10*d*m*mp**2 + m**5*(3 - 2*d*mp**2) + m**3*(3 + 20*d*mp**2)) + l**4*(m + m**7 - 10*d*m*mp**2 + m**5*(3 - 2*d*mp**2) + m**3*(3 + 20*d*mp**2))))))/((1 + l**2)*(1 + m**2)*(-2*(1 + l**2)**2*(1 + m**2)**2 - 4*a*b*c**2*(m*(-1 + m**2) - 6*l**2*m*(-1 + m**2) + l**4*m*(-1 + m**2) + l*(1 - 6*m**2 + m**4) - l**3*(1 - 6*m**2 + m**4)) + 4*a*b*d**2*(m*(-1 + m**2) - 6*l**2*m*(-1 + m**2) + l**4*m*(-1 + m**2) + l*(1 - 6*m**2 + m**4) - l**3*(1 - 6*m**2 + m**4)) + 2*a*b*c*d*(1 - 6*m**2 + m**4 - 16*l*m*(-1 + m**2) + 16*l**3*m*(-1 + m**2) - 6*l**2*(1 - 6*m**2 + m**4) + l**4*(1 - 6*m**2 + m**4)) + a**2*(2*d*(l - m + l**2*m - l*m**2) + c*(1 + 4*l*m - m**2 + l**2*(-1 + m**2)))**2 + b**2*(2*c*(m - l**2*m + l*(-1 + m**2)) + d*(1 + 4*l*m - m**2 + l**2*(-1 + m**2)))**2))
  ]
}
  </script>


</svg>