Skip to content

Instantly share code, notes, and snippets.

@travisdoesmath
Last active January 6, 2018 13:40
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save travisdoesmath/c0d08396a0bf8d5de8519ded00bb7866 to your computer and use it in GitHub Desktop.
Save travisdoesmath/c0d08396a0bf8d5de8519ded00bb7866 to your computer and use it in GitHub Desktop.
Double Pendulums are Chaotic

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>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment