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.
Last active
July 23, 2024 01:32
-
-
Save travisdoesmath/1206c28b72589d7269b670df05730705 to your computer and use it in GitHub Desktop.
(N>1)-tuple pendulums are chaotic
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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); | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
<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> |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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