Inspired by (and copied directly from) Dan Shiffman’s Coding Train video, I ported an algorithm for the motion of a double pendulum to run on the server so that it could be used to modulate parameters (GitHub link below).
This demo drives some FM parameters with the x and y positions of the masses.
That’s quite interesting, we should collect all those cool algorithms, translate into sc code, and put in a repo.
A very nice example !
Probably very wrong, but outputs something similar:
https://www.myphysicslab.com/springs/single-spring-en.html
(
k = 0.5; // Spring constant
m = 0.7; // Mass
b = 0.3; // Damping coefficient
~dxdt = { |xx| xx[1] }; // position
~dvdt = { |xx| (((k * -1) * xx[0]) / m) - (b * xx[1] / m) }; // velocity
~odes = { |y, t| [~dxdt.(y), ~dvdt.(y)] }; // both
/*~eulerMethod = { |func, y, t, tEnd, h|
var derivatives, results;
while { t <= tEnd } {
derivatives = func.value(y, t);
y = y + (derivatives * h);
t = t + h;
results = results.add([t, y[0], y[1]]);
};
results;
};*/
~rk4Method = { |func, y, t, tEnd, h|
var k1, k2, k3, k4;
var results = [];
while { t <= tEnd } {
k1 = h * func.(y, t);
k2 = h * func.(y + (k1 * 0.5), t + (h * 0.5));
k3 = h * func.(y + (k2 * 0.5), t + (h * 0.5));
k4 = h * func.(y + k3, t + h);
y = y + ((1/6) * (k1 + (2*k2) + (2*k3) + k4));
t = t + h;
results = results.add([t, y[0], y[1]]);
};
results;
};
// function, [initialPosition, initialVelocity], startTime, endTime, stepSize
~results = ~rk4Method.(~odes, [-1, 2], 0, 40, 0.025);
)
( //position
a = ~results.collect { |item| item[1] };
a.plot;
)
(//velocity
b = ~results.collect { |item| item[2] };
b.plot;
)
my favourite metaphor!
ooh nice. i’ll check this one out too!
I was not fishing for creds, I was just my usual excited puppy self 
very cool. reminded me of the difference of exploring bifurcation points of non-linear functions instead of using pure randomness according to this paper by tom mudd Nonlinear Dynamics In Musical Interactions - Open Research Online. maybe with the MLPRegressor. Probably also cool to use this to drive the x and y axis of the MLP to control a bunch of Synth paramters, or both.
BTW, Fb1_ODE from miSCellaneous_lib is an interface for numerical resolving of ODE systems (with a reliable symplectic procedure as default but it allows to choose other procedures as well). Pendulum examples are included in the help file as well, though not the double pendulum.
SMC paper: Fb1_ODE – an Interface for Synthesis and Sound Processing with Ordinary Differential Equations in SuperCollider
Video, summary and demo: https://www.youtube.com/watch?v=AcuKXxw4vrE
Before writing the script with odes, I should have known someone had already done this robustly. Thanks for sharing.
Nice!! I also have one https://github.com/aleksandarkoruga/public-extensions/tree/main/DoublePendulum this is a compiled uGen, source code is in another repo on the same github.
Cheers!
I also really like this one: myPhysicsLab Chaotic Pendulum
But I’m having some trouble implementing it (code below). It is acting quite chaotically (not in a good way…)
I think it has to do with the time step and feedback nature of the algorithm. I think there’s too much feedback happening causing it to literally spiral out of control.
Any one have any insights?
Code on the server:
/*
Formula from: https://www.myphysicslab.com/pendulum/chaotic-pendulum-en.html
*/
s.boot;
(
fork({
var win, uv, thetaDraw = 0, rDraw = 0, mDraw = 0;
var synth;
~cb = Bus.control(s,2);
s.sync;
~cb.setn([pi/2,0]);
Window.closeAll;
win = Window("Chaotic Pendulum",Rect(500,500,700,700));
uv = UserView().drawFunc_{
var origin = Point(uv.bounds.width/2,uv.bounds.height * 0.5);
var scale = 200;
var circleWidth, halfCircleWidth;
var pt1 = origin + Point(sin(thetaDraw) * rDraw * scale,cos(thetaDraw) * rDraw * scale);
Pen.line(origin,pt1);
Pen.stroke;
circleWidth = mDraw * scale * 0.2;
halfCircleWidth = circleWidth / 2;
Pen.circle(Rect(pt1.x-halfCircleWidth,pt1.y-halfCircleWidth,circleWidth,circleWidth));
Pen.fill;
};
win.layout = VLayout(
Button().states_([["Random Position"]]).action_{~cb.setn([0,0])},
uv
);
~chaoticPendulum = {
arg len = 1, gravity = 1, mass = 1, forceAmplitude = 1.5, damping = 0.5, forceFreq = 0.67;
var theta, thetaP, thetaPP;
var time = Sweep.kr;
var preplus, num, den, torque, quotient;
# theta, thetaP = In.kr(~cb,2);
preplus = -1 * (gravity/len) * sin(theta);
torque = forceAmplitude * cos(forceFreq * time);
num = (-1 * damping * thetaP) + torque;
den = mass * len * len;
quotient = num/den;
thetaPP = preplus + quotient;
// [thetaPP,thetaP,theta].poll;
thetaP = thetaP + thetaPP;
theta = theta + thetaP;
[quotient,thetaPP].poll;
Out.kr(~cb,[theta,thetaP]);
SendReply.kr(Impulse.kr(30),"/chaoticPendulum",[theta,len,mass]);
[time,preplus,torque,num,den,quotient,thetaPP,thetaP,theta]
};
synth = ~chaoticPendulum.play;
OSCdef(\chaoticPendulum,{
arg msg;
thetaDraw = msg[3];
rDraw = msg[4];
mDraw = msg[5];
defer{uv.refresh};
},"/chaoticPendulum");
win.onClose_({synth.free});
win.front;
},AppClock);
)
I also tried implementing in the language, just to trouble shoot, but it has the same results!:
/*
Formula from: https://www.myphysicslab.com/pendulum/chaotic-pendulum-en.html
*/
s.boot;
(
fork({
var win, uv, time = 0, timeStep = 0.025, theta = 0, thetaP = 0;
Window.closeAll;
win = Window("Chaotic Pendulum",Rect(500,500,700,700));
uv = UserView().drawFunc_{
var origin = Point(uv.bounds.width/2,uv.bounds.height * 0.5);
var scale = 200;
var circleWidth, halfCircleWidth, pt1;
var len = 1, gravity = 1, mass = 1, forceAmplitude = 1.47, damping = 0.5, forceFreq = 0.67;
var thetaPP;
var preplus, num, den, torque, quotient;
var dd, mlsq;
dd = -1 * (gravity / len) * sin(theta);
mlsq = mass * len * len;
dd = dd + (-1 * (damping/mlsq) * thetaP);
dd = dd + ((forceAmplitude/mlsq) * cos(forceFreq * time));
thetaP = thetaP + dd;
theta = theta + thetaP;
// draw
pt1 = origin + Point(sin(theta) * len * scale,cos(theta) * len * scale);
Pen.line(origin,pt1);
Pen.stroke;
circleWidth = mass * scale * 0.2;
halfCircleWidth = circleWidth / 2;
Pen.circle(Rect(pt1.x-halfCircleWidth,pt1.y-halfCircleWidth,circleWidth,circleWidth));
Pen.fill;
};
win.layout = VLayout(
Button().states_([["Random Position"]]).action_{"random position not yet implemented".postln;},
uv
);
win.front;
inf.do{
uv.refresh;
time = time + timeStep;
time.postln;
theta.postln;
timeStep.wait;
};
},AppClock);
)
I think the mistake is in how the ODE is solved:
thetaP = thetaP + dd;
theta = theta + thetaP;
This is the forward Euler method, but the time step is 1 second. Like all numerical ODE solvers, the time step has to be sufficiently small. I think this would fix it:
thetaP = thetaP + (dd * dt);
theta = theta + (thetaP * dt);
where dt is a small value like 0.001. MyPhysicsLab uses RK4, but if Euler works then don’t worry about it.
Thanks @nathan! that did it.
I had tried to scaled down the time step, but now I realized I was only scaling down the acceleration (as in the double pendulum example above), thinking it would propagate to the velocity and position, but alas. Again some FM for a demo:
/*
Formula from: https://www.myphysicslab.com/pendulum/chaotic-pendulum-en.html
*/
s.boot;
(
fork({
var win, uv, thetaDraw = 0, rDraw = 0, mDraw = 0;
var synth;
~cb = Bus.control(s,2);
s.sync;
~cb.setn([pi/2,0]);
Window.closeAll;
win = Window("Chaotic Pendulum",Rect(500,500,700,700));
uv = UserView().drawFunc_{
var origin = Point(uv.bounds.width/2,uv.bounds.height * 0.5);
var scale = 200;
var circleWidth, halfCircleWidth;
var pt1 = origin + Point(cos(thetaDraw) * rDraw * scale,sin(thetaDraw) * rDraw * scale);
Pen.line(origin,pt1);
Pen.stroke;
circleWidth = mDraw * scale * 0.2;
halfCircleWidth = circleWidth / 2;
Pen.circle(Rect(pt1.x-halfCircleWidth,pt1.y-halfCircleWidth,circleWidth,circleWidth));
Pen.fill;
};
win.layout = VLayout(
Button().states_([["Random Position"]]).action_{~cb.setn([rrand(0,2pi),0])},
uv
);
~chaoticPendulum = {
arg len = 1, gravity = 1, mass = 1, forceAmplitude = 1.15, damping = 0.5, forceFreq = 0.667;
var theta, thetaP, thetaPP;
var time = Sweep.kr;
var sig;
var preplus, num, den, torque, quotient, dt = 0.003;
# theta, thetaP = In.kr(~cb,2);
preplus = -1 * (gravity/len) * sin(theta);
torque = forceAmplitude * cos(forceFreq * time);
num = (-1 * damping * thetaP) + torque;
den = mass * len * len;
quotient = num/den;
thetaPP = preplus + quotient;
thetaP = thetaP + (thetaPP * dt);
theta = theta + (thetaP * dt);
Out.kr(~cb,[theta,thetaP]);
SendReply.kr(Impulse.kr(30),"/chaoticPendulum",[theta,len,mass]);
sig = SinOsc.ar(theta.fold(0,2pi).linlin(0,2pi,80,4000) * 2.pow(SinOsc.ar(thetaP.linlin(-1,1,20,200)) * 5 * thetaPP));
sig = sig * 0.1;
sig.dup;
};
synth = ~chaoticPendulum.play;
OSCdef(\chaoticPendulum,{
arg msg;
thetaDraw = msg[3];
rDraw = msg[4];
mDraw = msg[5];
defer{uv.refresh};
},"/chaoticPendulum");
win.onClose_({synth.free});
win.front;
},AppClock);
)