-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdouble-pendulum.html
117 lines (98 loc) · 8.54 KB
/
double-pendulum.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
<!DOCTYPE html>
<html lang=en>
<head>
<meta charset="utf-8">
<title>An algebraic lagrangian for the double pendulum</title>
</head>
<body>
<p><img src="https://i.imgflip.com/5yqsyy.jpg"/></p>
<canvas id="double-pendulum" width=400 height=400>
Your browser can't display canvas.
</canvas>
<script>
'use strict;'
function phi(a,b,c,d,lambda,mu,lambdadot,mudot) {
return [
[-(2*b*lambda)/(1 + lambda**2) + (a*(1 - lambda**2))/(1 + lambda**2), (2*a*lambda)/(1 + lambda**2) + (b*(1 - lambda**2))/(1 + lambda**2)],
[-((2*b*lambda)/(1 + lambda**2)) + (a*(1 - lambda**2))/(1 + lambda**2) - (2*d*mu)/(1 + mu**2) + (c*(1 - mu**2))/(1 + mu**2), (2*a*lambda)/(1 + lambda**2) + (b*(1 - lambda**2))/(1 + lambda**2) + (2*c*mu)/(1 + mu**2) + (d*(1 - mu**2))/(1 + mu**2)],
[(4*b*lambda**2*lambdadot)/(1 + lambda**2)**2 - (2*a*lambda*(1 - lambda**2)*lambdadot)/(1 + lambda**2)**2 - (2*b*lambdadot)/(1 + lambda**2) - (2*a*lambda*lambdadot)/(1 + lambda**2), -((4*a*lambda**2*lambdadot)/(1 + lambda**2)**2) - (2*b*lambda*(1 - lambda**2)*lambdadot)/(1 + lambda**2)**2 + (2*a*lambdadot)/(1 + lambda**2) - (2*b*lambda*lambdadot)/(1 + lambda**2)],
[(4*b*lambda**2*lambdadot)/(1 + lambda**2)**2 - (2*a*lambda*(1 - lambda**2)*lambdadot)/(1 + lambda**2)**2 - (2*b*lambdadot)/(1 + lambda**2) - (2*a*lambda*lambdadot)/(1 + lambda**2) + (4*d*mu**2*mudot)/(1 + mu**2)**2 - (2*c*mu*(1 - mu**2)*mudot)/(1 + mu**2)**2 - (2*d*mudot)/(1 + mu**2) - (2*c*mu*mudot)/(1 + mu**2),-((4*a*lambda**2*lambdadot)/(1 + lambda**2)**2) - (2*b*lambda*(1 - lambda**2)*lambdadot)/(1 + lambda**2)**2 + (2*a*lambdadot)/(1 + lambda**2) - (2*b*lambda*lambdadot)/(1 + lambda**2) - (4*c*mu**2*mudot)/(1 + mu**2)**2 - (2*d*mu*(1 - mu**2)*mudot)/(1 + mu**2)**2 + (2*c*mudot)/(1 + mu**2) - (2*d*mu*mudot)/(1 + mu**2)]
]
}
function psi(a,b,c,d,lambda,mu,lambdadot,mudot) {
let $phi = phi(a,b,c,d,lambda,mu,lambdadot,mudot),
$a = $phi[0][0],
$b = $phi[0][1],
$c = $phi[1][0] - $a,
$d = $phi[1][1] - $b,
v1x = $phi[2][0],
v1y = $phi[2][1],
v2x = $phi[3][0],
v2y = $phi[3][1];
return [
[$a, $b, $c, $d],
[
0, 0,
(-$b*v1x + $a*v1y)/(2*($a**2+$b**2)),
(-$d*(v2x-v1x) + $c*(v2y-v1y))/(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))
]
}
function rk4(dt, x, D) {
let indices = [...Array(x.length).keys()],
map = f => indices.map(f),
dtx = $ => map( i => $[i]*dt ),
a = dtx(D(...x)),
b = dtx(D(...map(i => x[i] + a[i]/2))),
c = dtx(D(...map(i => x[i] + b[i]/2))),
d = dtx(D(...map(i => x[i] + c[i])));
return map(i => (a[i] + 2*b[i] + 2*c[i] + d[i])/6);
}
var canvas = document.getElementById("double-pendulum");
var context = canvas.getContext("2d");
function draw(...balls) {
let x = context,
s = (canvas.width + canvas.height)/2/4,
r = .1;
x.save();
x.clearRect(0, 0, canvas.width, canvas.height);
x.translate(canvas.width/2, canvas.height/2);
x.fillRect(0,0,r,r);
x.scale(s,-s);
x.beginPath();
x.lineWidth = 2/s;
x.moveTo(0,0);
for (let ball of balls)
x.lineTo(...ball);
x.stroke();
x.fillRect(-r,-r,2*r,2*r);
for (let ball of balls) {
x.beginPath();
x.arc(...ball,r,0,2*355/113);
x.fill();
}
x.restore();
}
var abcd = [0,1,1,0];
var x = [0,0,0,0];
(function animate(then) {
let now = Date.now()/1000,
dt = now - then,
dx = rk4(dt, x, ddot(...abcd)),
$psi = psi(...abcd,...[...x.keys()].map(i => x[i]+dx[i]));
abcd = $psi[0];
x = $psi[1];
draw(...phi(...abcd,...x).slice(0,2));
window.requestAnimationFrame(() => animate(now));
})(Date.now()/1000);
</script>
</body>
</html>