Showing the sensitive dependence to initial conditions, forty double pendulums start in very nearly the same spot, but quickly diverge.
Inspired by this triple pendulum GIF
Showing the sensitive dependence to initial conditions, forty double pendulums start in very nearly the same spot, but quickly diverge.
Inspired by this triple pendulum GIF
Array.prototype.scalarMultiply = function(x) { return this.map(function(d) { return x * d; }); } | |
var l1 = 1, l2 = 1, | |
m1 = 1, m2 = 1, | |
G = 9.8; | |
var theta1 = 0.49*Math.PI, | |
theta2 = 1.0*Math.PI, | |
p1 = 0, | |
p2 = 0; | |
var nPendula = 40, | |
initialRange = 1e-4; | |
theta1dot = function(theta1, theta2, p1, p2) { | |
return (p1*l2 - p2*l1*Math.cos(theta1 - theta2))/(l1**2*l2*(m1 + m2*Math.sin(theta1 - theta2)**2)); | |
} | |
theta2dot = function(theta1, theta2, p1, p2) { | |
return (p2*(m1+m2)*l1 - p1*m2*l2*Math.cos(theta1 - theta2))/(m2*l1*l2**2*(m1 + m2*Math.sin(theta1 - theta2)**2)); | |
} | |
p1dot = function(theta1, theta2, p1, p2) { | |
var A1 = (p1*p2*Math.sin(theta1 - theta2))/(l1*l2*(m1 + m2*Math.sin(theta1 - theta2)**2)), | |
A2 = (p1**2*m2*l2**2 - 2*p1*p2*m2*l1*l2*Math.cos(theta1 - theta2) + p2**2*(m1 + m2)*l1**2)*Math.sin(2*(theta1 - theta2))/(2*l1**2*l2**2*(m1 + m2*Math.sin(theta1 - theta2)**2)**2); | |
return -(m1 + m2)*G*l1*Math.sin(theta1) - A1 + A2; | |
} | |
p2dot = function(theta1, theta2, p1, p2) { | |
var A1 = (p1*p2*Math.sin(theta1 - theta2))/(l1*l2*(m1 + m2*Math.sin(theta1 - theta2)**2)), | |
A2 = (p1**2*m2*l2**2 - 2*p1*p2*m2*l1*l2*Math.cos(theta1 - theta2) + p2**2*(m1 + m2)*l1**2)*Math.sin(2*(theta1 - theta2))/(2*l1**2*l2**2*(m1 + m2*Math.sin(theta1 - theta2)**2)**2); | |
return -m2*G*l2*Math.sin(theta2) + A1 - A2; | |
} | |
var Z = [theta1, theta2, p1, p2], | |
oldZ = Z; | |
deltas = d3.range(-initialRange/2, initialRange/2,initialRange/nPendula); | |
Zs = deltas.map(x=>[Z[0],Z[1]+x,Z[2],Z[3]]) | |
f = function(Z) { | |
return [theta1dot(Z[0], Z[1], Z[2], Z[3]), theta2dot(Z[0], Z[1], Z[2], Z[3]), p1dot(Z[0], Z[1], Z[2], Z[3]), p2dot(Z[0], Z[1], Z[2], Z[3])]; | |
} | |
RK4 = function(Z_n, f, tau) { | |
var Y1 = f(Z_n).scalarMultiply(tau); | |
var Y2 = f([Z_n[0] + 0.5*Y1[0], Z_n[1] + 0.5*Y1[1], Z_n[2] + 0.5*Y1[2], Z_n[3] + 0.5*Y1[3]]).scalarMultiply(tau); | |
var Y3 = f([Z_n[0] + 0.5*Y2[0], Z_n[1] + 0.5*Y2[1], Z_n[2] + 0.5*Y2[2], Z_n[3] + 0.5*Y2[3]]).scalarMultiply(tau); | |
var Y4 = f([Z_n[0] + Y3[0], Z_n[1] + Y3[1], Z_n[2] + Y3[2], Z_n[3] + Y3[3]]).scalarMultiply(tau); | |
return [ | |
Z_n[0] + Y1[0]/6 + Y2[0]/3 + Y3[0]/3 + Y4[0]/6, | |
Z_n[1] + Y1[1]/6 + Y2[1]/3 + Y3[1]/3 + Y4[1]/6, | |
Z_n[2] + Y1[2]/6 + Y2[2]/3 + Y3[2]/3 + Y4[2]/6, | |
Z_n[3] + Y1[3]/6 + Y2[3]/3 + Y3[3]/3 + Y4[3]/6, | |
] | |
} | |
getCoords = function(Z) { | |
return { | |
'x1':l1*Math.sin(Z[0]), | |
'y1':l1*Math.cos(Z[0]), | |
'x2':l1*Math.sin(Z[0]) + l2*Math.sin(Z[1]), | |
'y2':l1*Math.cos(Z[0]) + l2*Math.cos(Z[1]) | |
}; | |
} |
<html> | |
<head> | |
<style> | |
.bob { | |
fill: #000; | |
} | |
.shaft { | |
stroke: #000; | |
stroke-width: 2px; | |
} | |
.trace { | |
fill: none; | |
stroke: #F00; | |
stroke-width: 2px; | |
} | |
svg { | |
position:absolute; | |
top:0px; | |
left:0px; | |
} | |
canvas { | |
position:absolute; | |
top:0px; | |
left:0px; | |
} | |
</style> | |
<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/4.7.4/d3.js"></script> | |
<script src="./double_pendulum.js"></script> | |
</head> | |
<body> | |
<canvas width="960" height="500"></canvas> | |
<svg width="960" height="500"></svg> | |
</body> | |
<script> | |
var svg = d3.select("svg") | |
width = +svg.attr("width"), | |
height = +svg.attr("height"), | |
g = svg.append("g").attr("transform", "translate(" + width*.5 + "," + height*.5 + ")"); | |
color = d3.scaleSequential(d3.interpolateRainbow).domain([0, nPendula]); | |
var canvas = d3.select("canvas"); | |
var context = canvas.node().getContext('2d'); | |
var scale = d3.scaleLinear().domain([0,1]).range([0,100]) | |
var oldZs = Zs; | |
var path = d3.line() | |
.x(function(d) { return scale(l1*Math.sin(d[0])+l2*Math.sin(d[1])); }) | |
.y(function(d) { return scale(l1*Math.cos(d[0])+l2*Math.cos(d[1])); }) | |
var update = function(Zs, oldZs) { | |
var coords = Zs.map(Z=>getCoords(Z)), | |
oldCoords = oldZs.map(Z=>getCoords(Z)); | |
for (var i = coords.length - 1; i >= 0; i--) { | |
context.beginPath(); | |
context.strokeStyle = color(i); | |
context.lineWidth = 2; | |
context.moveTo(scale(oldCoords[i].x2) + width/2, scale(oldCoords[i].y2) + height/2); | |
context.lineTo(scale(coords[i].x2) + width/2, scale(coords[i].y2) + height/2); | |
context.stroke(); | |
} | |
var pendulum = g.selectAll(".pendulum").data(coords, function(d, i) { return i; }) | |
var pendulumEnter = pendulum.enter() | |
.append("g").attr("class","pendulum") | |
pendulumEnter.append("line").attr("class", "firstShaft shaft") | |
pendulumEnter.append("line").attr("class", "secondShaft shaft") | |
pendulumEnter.append("circle").attr("class", "firstBob bob").attr("r",3) | |
pendulumEnter.append("circle").attr("class", "secondBob bob").attr("r",3) | |
var shaft1 = pendulum.select(".firstShaft") | |
.attr("x1", 0) | |
.attr("y1", 0) | |
.attr("x2", function(d) { return +scale(d.x1); }) | |
.attr("y2", function(d) { return +scale(d.y1); }) | |
var shaft2 = pendulum.select(".secondShaft") | |
.attr("x1", function(d) { return +scale(d.x1); }) | |
.attr("y1", function(d) { return +scale(d.y1); }) | |
.attr("x2", function(d) { return +scale(d.x2); }) | |
.attr("y2", function(d) { return +scale(d.y2); }) | |
var bob1 = pendulum.select(".firstBob") | |
.attr("cx",function(d) { return +scale(d.x1); }) | |
.attr("cy",function(d) { return +scale(d.y1); }) | |
var bob2 = pendulum.select(".secondBob") | |
.attr("cx",function(d) { return +scale(d.x2); }) | |
.attr("cy",function(d) { return +scale(d.y2); }) | |
} | |
var run = setInterval(function() { update(Zs, oldZs); oldZs = Zs; Zs = Zs.map(Z=>RK4(Z, f, 0.005)); }, 2); | |
</script> |