# Double Pendulum as LFO

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.

8 Likes

That’s quite interesting, we should collect all those cool algorithms, translate into sc code, and put in a repo.

1 Like

A very nice example !

2 Likes

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 }; // position
~dvdt = { |xx| (((k * -1) * xx) / m) - (b * xx / 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;
};*/

~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;
};
// function, [initialPosition, initialVelocity], startTime, endTime, stepSize
~results = ~rk4Method.(~odes, [-1, 2], 0, 40, 0.025);
)

( //position
a = ~results.collect { |item| item };
a.plot;
)

(//velocity
b = ~results.collect { |item| item };
b.plot;
)

``````

my favourite metaphor!

3 Likes
2 Likes

i should have @-ed you @tremblap for inspiration as well!

1 Like

ooh nice. i’ll check this one out too!

1 Like

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.

2 Likes

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.

5 Likes

Before writing the script with odes, I should have known someone had already done this robustly. Thanks for sharing.

1 Like
1 Like

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!

2 Likes

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]);

[time,preplus,torque,num,den,quotient,thetaPP,thetaP,theta]
};

synth = ~chaoticPendulum.play;

OSCdef(\chaoticPendulum,{
arg msg;
rDraw = msg;
mDraw = msg;
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.

1 Like

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]);

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;