Skip to content

Instantly share code, notes, and snippets.

@travisdoesmath
Last active December 28, 2017 20:38
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/edcd5b3b02a9d924ce3ea0fb48b83f8e to your computer and use it in GitHub Desktop.
Save travisdoesmath/edcd5b3b02a9d924ce3ea0fb48b83f8e to your computer and use it in GitHub Desktop.
Double pendulum simulator

A double pendulum simulator using a Runge-Kutta integrator.

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.5*Math.PI + Math.random()/10,
theta2 = 1.0*Math.PI + Math.random()/10,
p1 = 0,
p2 = 0;
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;
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])];
}
f(Z)
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: #F00;
}
.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 + ")");
var canvas = d3.select("canvas");
var context = canvas.node().getContext('2d');
var scale = d3.scaleLinear().domain([0,1]).range([0,100])
var oldZ;
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 tracePath = .append("path").datum(Zs).attr("d", function(d) { return path(d); }).attr("class","trace");
var pendulum = g.append("g").attr("class","pendulum");
var shaft1 = pendulum.append("line").attr("class", "shaft").attr("id","shaft1").attr("x1",0).attr("y1", 0).attr("x2",0).attr("y2",0)
var shaft2 = pendulum.append("line").attr("class", "shaft").attr("id","shaft2").attr("x1",0).attr("y1", 0).attr("x2",0).attr("y2",0)
var bob1 = pendulum.append("circle").attr("class", "bob").attr("id","bob1").attr("cx",0).attr("cy",0).attr("r",10)
var bob2 = pendulum.append("circle").attr("class", "bob").attr("id","bob2").attr("cx",0).attr("cy",0).attr("r",10)
var update = function(Z, oldZ) {
var coords = getCoords(Z),
oldCoords = getCoords(oldZ);
//tracePath.datum(Zs).attr("d", function(d) { return path(d); }).attr("class","trace");
context.beginPath();
context.globalAlpha = 0.2
context.strokeStyle = "#00F";
context.moveTo(scale(oldCoords.x2) + width/2, scale(oldCoords.y2) + height/2);
context.lineTo(scale(coords.x2) + width/2, scale(coords.y2) + height/2);
context.stroke();
context.globalAlpha = 1;
d3.select("#shaft1")
.attr("x1", +0)
.attr("y1", +0)
.attr("x2", +scale(coords.x1))
.attr("y2", +scale(coords.y1))
d3.select("#shaft2")
.attr("x1", +scale(coords.x1))
.attr("y1", +scale(coords.y1))
.attr("x2", +scale(coords.x2))
.attr("y2", +scale(coords.y2))
d3.select("#bob1").attr("cx", scale(coords.x1)).attr("cy", scale(coords.y1))
d3.select("#bob2").attr("cx", scale(coords.x2)).attr("cy", scale(coords.y2))
}
var run = setInterval(function() { update(Z, oldZ); oldZ = Z; 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