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!
Can improvements be made so it is still human readable?
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.
"ExplicitRungeKutta"
method ofNDSolve
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.) $\endgroup$ – xzczd Jan 1 '16 at 5:38