Navigation Menu

Skip to content

Instantly share code, notes, and snippets.

@travisdoesmath
Last active November 24, 2020 16:26
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/1206c28b72589d7269b670df05730705 to your computer and use it in GitHub Desktop.
Save travisdoesmath/1206c28b72589d7269b670df05730705 to your computer and use it in GitHub Desktop.
(N>1)-tuple pendulums are chaotic

While doing the math for the triple pendulum, I kept losing track of terms, so I solved for the general case of N pendulums and set N = 3 (to keep the math simple, all masses and lengths are assumed to be 1). Here are a single, double, quadruple, octuple, and a sexdecuple (16-tuple) pendulum. Click to restart with a new starting angle.

let test = new TuplePendulum({N: 4, thetas:[1, 2, 3], omegas:[0.1,0.2,0.3]})
var nPendulums = 5
var pendulums = d3.range(nPendulums).map(x => new TuplePendulum({N: 2**x, thetas: [0.5 * Math.PI]}))
var fadeBackground = false;
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, nPendulums]);
svg.on('click', e => {
var mousePos = d3.mouse(svg.node());
reset(mousePos);
});
var canvas = d3.select("canvas");
var context = canvas.node().getContext('2d');
var scale = d3.scaleLinear().domain([0,1]).range([0,200])
var update = function() {
var oldCoords = pendulums.map(p => p.getCoords());
pendulums.forEach(p => p.time_step(0.005));
var coords = pendulums.map(p => p.getCoords());
draw(oldCoords, coords);
}
var trailOpacity = 1;
var maxThetaDelta = 0;
var opacityScale = d3.scaleLinear().domain([0, 2*Math.PI]).range([1, 0])
var draw = function(oldCoords, coords) {
if (maxThetaDelta < 2*Math.PI) {
if (fadeBackground) {
maxThetaDelta = Math.max(maxThetaDelta, Math.abs(d3.max(pendulums, d => d.theta1) - d3.min(pendulums, d => d.theta1)))
//trailOpacity -= maxThetaDelta / 1500;
trailOpacity = opacityScale(maxThetaDelta)
//trailOpacity = opacityScale(Math.abs(pendulums[nPendulums - 1].theta1 - pendulums[0].theta1))
}
canvas.style('opacity', trailOpacity);
}
var pendulum = g.selectAll(".pendulum").data(coords, (d, i) => i)
pendulum.enter()
.append("g").attr("class","pendulum")
.attr('stroke', (d, i) => color(i))
.attr('fill', (d, i) => color(i))
var shafts = pendulum.selectAll('.shaft').data((d, i) => d)
shafts.enter().append('line').attr("class", "shaft")
.merge(shafts)
.attr('test', d => console.log(d))
.attr("x1", d => scale(d.x1))
.attr("y1", d => scale(d.y1))
.attr("x2", d => scale(d.x2))
.attr("y2", d => scale(d.y2))
var bobs = pendulum.selectAll('.bob').data(d => d)
bobs.enter().append('circle').attr('class', 'bob').attr('r',3)
.merge(bobs)
.attr("cx", d => scale(d.x2))
.attr("cy", d => scale(d.y2))
for (let i = 0; i < coords.length - 1; i++) {
// context.beginPath();
// context.strokeStyle = color(i);
// context.lineWidth = 2;
// context.moveTo(scale(oldCoords[i][oldCoords[i].length - 1].x2) + width/2, scale(oldCoords[i][oldCoords[i].length - 1].y2) + height/2);
// context.lineTo(scale(coords[i][coords[i].length - 1].x2) + width/2, scale(coords[i][coords[i].length - 1].y2) + height/2);
// context.stroke();
}
}
var reset = function(mousePos) {
console.log(mousePos)
var theta = 0.5*Math.PI + Math.atan2(height/2 - mousePos[1], mousePos[0] - width/2)
// trailOpacity = 1;
// maxThetaDelta = 0;
pendulums = d3.range(nPendulums).map(x => new TuplePendulum({N: 2**x, thetas: [theta]}));
context.clearRect(0, 0, width, height);
}
var run = setInterval(() => { update() }, 2);
<html>
<head>
<style>
.shaft {
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/5.16.0/d3.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/mathjs/7.6.0/math.min.js"></script>
<script src="./tuple_pendulum.js"></script>
<!-- <script src="./double_pendulum.js"></script> -->
</head>
<body>
<canvas width="960" height="500"></canvas>
<svg width="960" height="500"></svg>
</body>
<script src="app.js"></script>
class TuplePendulum {
constructor(opts) {
this.N = 3;
this.G = -9.8;
this.thetas=[0.75*Math.PI];
this.omegas=[0];
if (opts) ['N','G','thetas','omegas'].map(k => opts[k] ? this[k] = opts[k] : null)
if (this.thetas.length != this.N) {
this.thetas = new Array(this.N).fill(0).map((x, i) => this.thetas[i % this.thetas.length])
}
if (this.omegas.length != this.N) {
this.omegas = new Array(this.N).fill(0).map((x, i) => this.omegas[i % this.omegas.length])
}
}
A(thetas = this.thetas) {
var M = [];
for (let i = 1; i <= this.N; i++) {
let row = []
for (let j = 1; j <= this.N; j++) {
row.push((this.N - Math.max(i, j) + 1)/(this.N - i + 1) * Math.cos(thetas[i-1] - thetas[j-1]))
}
M.push(row)
}
return M;
}
Ainv(thetas) {
return math.inv(this.A(thetas))
}
b(thetas = this.thetas, omegas = this.omegas) {
let v = []
for (let i = 1; i <= this.N; i++) {
let b_i = 0
for (let j = 1; j <= this.N; j++) {
let coef = (this.N - Math.max(i, j) + 1)/(this.N - i + 1)
let theta_i = thetas[i-1]
let theta_j = thetas[j-1]
let omega_j = omegas[j - 1]
let delta = coef * Math.sin(theta_j - theta_i) * omega_j ** 2
b_i += delta
}
b_i += this.G * Math.sin(thetas[i-1])
v.push(b_i)
}
return v;
}
lagrange_rhs([thetas, omegas]) {
var AinvB = math.multiply(this.Ainv(thetas), this.b(thetas, omegas))
return [omegas, AinvB]
}
time_step(dt) {
var k1 = this.lagrange_rhs([this.thetas, this.omegas])
var k2 = this.lagrange_rhs([math.add(this.thetas, k1[0].map(x => 0.5 * dt * x)), math.add(this.omegas, k1[1].map(x => 0.5 * dt * x))])
var k3 = this.lagrange_rhs([math.add(this.thetas, k2[0].map(x => 0.5 * dt * x)), math.add(this.omegas, k2[1].map(x => 0.5 * dt * x))])
var k4 = this.lagrange_rhs([math.add(this.thetas, k3[0].map(x => 1.0 * dt * x)), math.add(this.omegas, k3[1].map(x => 1.0 * dt * x))])
let theta_deltas = math.add(k1[0], k2[0].map(x => 2 * x), k3[0].map(x => 2 * x), k4[0]).map(x => x * dt/6)
let omega_deltas = math.add(k1[1], k2[1].map(x => 2 * x), k3[1].map(x => 2 * x), k4[1]).map(x => x * dt/6)
this.thetas = math.add(this.thetas, theta_deltas)
this.omegas = math.add(this.omegas, omega_deltas)
}
getCoords() {
let x1 = 0, y1 = 0, coords = [];
for (let i = 1; i <= this.thetas.length; i++) {
let x2 = 0
let y2 = 0
for (let j = 0; j < i; j++) {
x2 += Math.sin(this.thetas[j]) / this.N
y2 += Math.cos(this.thetas[j]) / this.N
}
coords.push({x1: x1, y1: y1, x2: x2, y2: y2})
x1 = x2
y1 = y2
}
return coords;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment