# System of Four First Order Equations for Double Pendulum using RK-4 - Code Improvements

I have spent quite a bit of time implementing the double pendulum equations at the bottom of this web site using Runge-Kutta-4.

I am also quite aware of the built-in Runge-Kutta methods, but I need control over everything and the built in methods were not letting me do that (I could be wrong).

The code is as follows.

  ClearAll[x, u, y, v, f1, f2, f3, f4, pts1, pts2, pts3, pts4, pts5,
pts6, pts7, pts8]
g = 9.81; mass1 = 1.; mass2 = 1.; len1 = 1; len2 = 1;
f1[t_, x_, u_, y_, v_] := y;
f2[t_, x_, u_, y_, v_] := v;
f3[t_, x_, u_, y_,
v_] := (-g (2 mass1 + mass1) Sin[x] - mass2 g Sin[x - 2 u] -
2 Sin[x - u] mass2 (v^2 len2 +
y^2 len1 Cos[x - u]))/(len1 (2 mass1 + mass2 -
mass2 Cos[2 x - 2 u]));
f4[t_, x_, u_, y_,
v_] := (2 Sin[
x - u] (y^2 len1 (mass1 + mass2) + g (mass1 + mass2) Cos[x] +
v^2 len2 mass2 Cos[x - u]))/(len2 (2 mass1 + mass2 -
mass2 Cos[2 x - 2 u]));
x = 0.0 Pi; u = -0.25 Pi; y = 0; v = 0; t = 0; h = 0.01; n = 5000;

pts1 = {{x, y}}; pts2 = {{u, v}}; pts3 = {{x, u}}; pts4 = {{y,
v}}; pts5 = {{t, x}}; pts6 = {{t, u}}; pts7 = {{t, y}}; pts8 = {{t,
v}};
Do[
j1 = f1[t, x, u, y, v];
k1 = f2[t, x, u, y, v];
l1 = f3[t, x, u, y, v];
m1 = f4[t, x, u, y, v];
j2 = f1[t + h/2, x + h/2*j1, u + h/2*k1, y + h/2*l1, v + h/2*m1];
k2 = f2[t + h/2, x + h/2*j1, u + h/2*k1, y + h/2*l1, v + h/2*m1];
l2 = f3[t + h/2, x + h/2*j1, u + h/2*k1, y + h/2*l1, v + h/2*m1];
m2 = f4[t + h/2, x + h/2*j1, u + h/2*k1, y + h/2*l1, v + h/2*m1];
j3 = f1[t + h/2, x + h/2*j2, u + h/2*k2, y + h/2*l2, v + h/2*m2];
k3 = f2[t + h/2, x + h/2*j2, u + h/2*k2, y + h/2*l2, v + h/2*m2];
l3 = f3[t + h/2, x + h/2*j2, u + h/2*k2, y + h/2*l2, v + h/2*m2];
m3 = f4[t + h/2, x + h/2*j2, u + h/2*k2, y + h/2*l2, v + h/2*m2];
j4 = f1[t + h, x + h*j3, u + h*k3, y + h*l3, v + h*m3];
k4 = f2[t + h, x + h*j3, u + h*k3, y + h*l3, v + h*m3];
l4 = f3[t + h, x + h*j3, u + h*k3, y + h*l3, v + h*m3];
m4 = f4[t + h, x + h*j3, u + h*k3, y + h*l3, v + h*m3];
x = x + h*(j1 + 2*j2 + 2*j3 + j4)/6;
u = u + h*(k1 + 2*k2 + 2*k3 + k4)/6;
y = y + h*(l1 + 2*l2 + 2*l3 + l4)/6;
v = v + h*(m1 + 2*m2 + 2*m3 + m4)/6;
t = t + h;
{AppendTo[pts1, {x, y}]; AppendTo[pts2, {u, v}],
AppendTo[pts3, {x, u}], AppendTo[pts4, {y, v}],
AppendTo[pts5, {t, x}],
AppendTo[pts6, {t, u}] AppendTo[pts7, {t, y}],
AppendTo[pts8, {t, v}]},
{i, 1, n}]

ListPlot[pts1, Joined->True]
ListPlot[pts2, Joined->True]
ListPlot[pts3, Joined->True]
ListPlot[pts4, Joined->True]
ListPlot[pts5, Joined->True]
ListPlot[pts6, Joined->True]
ListPlot[pts7, Joined->True]
ListPlot[pts8, Joined->True]


I know it is not the prettiest coding and that is my question. It appears to be working as desired and seems to be pretty fast, but looks ugly!

Note: it was not written in the easiest fashion in order to allow changing the system. Ideally, I wanted to look into extending it to all the variants with How to deal with the condition that a function own many options?. There is also this beautiful Animation of double pendulum.

• “I am also quite aware of the built-in Runge-Kutta methods, but I need control over everything and the built in methods were not letting me do that (I could be wrong)” I'm afraid you're indeed wrong, as far as I can tell, the "ExplicitRungeKutta" method of NDSolve is quite open and transparent and it allows you to control over every thing. Have you read this and this answer? (They are just quoting from the document actually.) – xzczd Jan 1 '16 at 5:38
• Have you see reference.wolfram.com/language/example/DoublePendulum.html? Also related: (37265), (40122), (48059) – Michael E2 Jan 1 '16 at 5:42
• @MichaelE2: Yes, I did see that on the Wolfram site and even referenced one of the items you linked to in my post. I also saw the other two items you posted. Thanks – Moo Jan 1 '16 at 5:46
• @xzczd: That first link you sent is worth exploring and somehow I missed it when searching - I found several others that I had spent time on to make sure I wasn't rehashing old ground.. Thanks for that link. – Moo Jan 1 '16 at 5:48

Replace

j1 = f1[t, x, u, y, v];
k1 = f2[t, x, u, y, v];
l1 = f3[t, x, u, y, v];
m1 = f4[t, x, u, y, v];
j2 = f1[t + h/2, x + h/2*j1, u + h/2*k1, y + h/2*l1, v + h/2*m1];
k2 = f2[t + h/2, x + h/2*j1, u + h/2*k1, y + h/2*l1, v + h/2*m1];
l2 = f3[t + h/2, x + h/2*j1, u + h/2*k1, y + h/2*l1, v + h/2*m1];
m2 = f4[t + h/2, x + h/2*j1, u + h/2*k1, y + h/2*l1, v + h/2*m1];
j3 = f1[t + h/2, x + h/2*j2, u + h/2*k2, y + h/2*l2, v + h/2*m2];
k3 = f2[t + h/2, x + h/2*j2, u + h/2*k2, y + h/2*l2, v + h/2*m2];
l3 = f3[t + h/2, x + h/2*j2, u + h/2*k2, y + h/2*l2, v + h/2*m2];
m3 = f4[t + h/2, x + h/2*j2, u + h/2*k2, y + h/2*l2, v + h/2*m2];
j4 = f1[t + h, x + h*j3, u + h*k3, y + h*l3, v + h*m3];
k4 = f2[t + h, x + h*j3, u + h*k3, y + h*l3, v + h*m3];
l4 = f3[t + h, x + h*j3, u + h*k3, y + h*l3, v + h*m3];
m4 = f4[t + h, x + h*j3, u + h*k3, y + h*l3, v + h*m3];


with

{j1, k1, l1, m1} = Map[#[t, x, u, y, v]&, {f1,f2,f3,f4}];
{j2, k2, l2, m2} = Map[#[t+h/2, x+h/2*j1, u+h/2*k1, y+h/2*l1, v+h/2*m1]&, {f1,f2,f3,f4}];
{j3, k3, l3, m3} = Map[#[t+h/2, x+h/2*j2, u+h/2*k2, y+h/2*l2, v+h/2*m2]&, {f1,f2,f3,f4}];
{j4, k4, l4, m4} = Map[#[t+h, x+h*j3, u+h*k3, y+h*l3, v+h*m3]&, {f1,f2,f3,f4}];


Or you can use Through to accomplish the same thing.

{j1, k1, l1, m1} = Through[{f1,f2,f3,f4}[t, x, u, y, v]];
{j2, k2, l2, m2} = Through[{f1,f2,f3,f4}[t+h/2, x+h/2*j1, u+h/2*k1, y+h/2*l1, v+h/2*m1]];
{j3, k3, l3, m3} = Through[{f1,f2,f3,f4}[t+h/2, x+h/2*j2, u+h/2*k2, y+h/2*l2, v+h/2*m2]];
{j4, k4, l4, m4} = Through[{f1,f2,f3,f4}[t+h, x+h*j3, u+h*k3, y+h*l3, v+h*m3]];