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

Question

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.

Answer 1

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