Difficulty using Excel to solve the double pendulum problem using RK4 to solve four simultaneous first order ODEs

Question

I have set up a spreadsheet in Excel to solve the double pendulum problem using RK4 on the four first order DEs. However it isn't working, as the energy in the system is increasing drastically, instead of staying constant. I have spent days (literally) checking the Lagrangian and the spreadsheet entries, and I can't find any faults. I wonder if I am mis-applying the RK4 procedure. I can't find any examples of this in Excel - only Matlab, which is quite different, and I don't have it! I'm desperate and so thought I would post here.

I have four first order DEs. The first two are simple W'=Y, X'=Z. Y' & Z' are both very long with terms in W, X, Y & Z. So, I have

W'=f(t,Y)

X'=f(t,Z)

Y'=f(t,X,Y,Z)

Z'=f(t,X,Y,Z)

I calculate

Wn+1 = Wn + h/6*(J1 + 2*J2 + 2*J3 + J4)

Xn+1 = Xn + h/6*(K1 + 2*K2 + 2*K3 + K4)

Yn+1 = Yn + h/6*(L1 + 2*L2 + 2*L3 + L4)

Zn+1 = Zn + h/6*(M1 + 2*M2 + 2*M3 + M4)

The values of these constants are worked out as follows: (Is this where I'm going wrong?)

J1 = f(tn, Yn)

J2 = f(tn + h/2, Yn + h/2*L1) 

J3 = f(tn + h/2, Yn + h/2*L2)

J4 = f(tn + h, Yn + h*L3)

The values for K1, K2, K3 & K4 are calculated in the same way using Zn & M1, M2 M3 & M4.

As Y is a function f(t, W, X, Y, Z), I calculate as follows:

L1 = f(tn, Wn, Xn, Yn, Zn)

L2 = f(tn + h/2, Wn + h/2*J1, Xn + h/2*K1, Yn + h/2*L1, Zn + h/2*M1)

L3 = f(tn + h/2, Wn + h/2*J2, Xn + h/2*K2, Yn + h/2*L2, Zn + h/2*M2)

L4 = f(tn + h, Wn + h*J3, Xn + h*K3, Yn + h*L3, Zn + h*M3).

The values for M1, M2, M3 & M4 used to calculate Z were computed in the exact same manner as for L1 etc. above.

If anyone can identify a flaw in my methodology and let me know I would be very grateful.

Answer 1

The goal is to use a fourth order Runge-Kutta Method to solve two converted $2nd$ order differential equations into four $1st$ order differential equations given by the double pendulum.

Depending on initial conditions, this system can show chaotic behavior, so play around with initial conditions (or seek them out on the web).

The system is given by:

$$ \left\{ \begin{array}{c} \theta_1' = \omega_1 \\ \theta_2' = \omega_2 \\ \omega_1' = \dfrac{-g(2 m_1 + m_2) \sin \theta_1 - m_2 g \sin(\theta_1 - 2 \theta_2) - 2 \sin(\theta_1 - \theta_2)m_2(\omega_2^2 L_2 + \omega_1^2 L_1 \cos(\theta_1 - \theta_2))}{L_1(2 m_1 + m_2 - m_2 \cos(2 \theta_1 - 2 \theta_2))} \\ \omega_2' = \dfrac{2 \sin(\theta_1 -\theta_2)(\omega_1^2 L_1(m_1 + m_2) + g(m_1 + m_2) \cos(\theta_1) + \omega_2^2 L_2 m_2 \cos(\theta_1 - \theta_2))}{L_2(2 m_1 + m_2 - m_2 \cos(2 \theta_1 - 2 \theta_2))} \end{array} \right. $$

Here is an implementation and it must be done in the order shown because of the dependencies. I changed some of the variable names, so you will have to adapt it, but it is exactly your problem. I used the equations from the web site you linked.

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-2u]
 -2Sin[x-u]mass2(v^2 len2 + y^2 len1 Cos[x-u]))/(len1(2 mass1
 + mass2 - mass2 Cos[2x-2u]))
f4[t, x, u, y, v] = (2Sin[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[2x-2u]))
x = 0.0 Pi
u = -0.25 Pi
y = 0
v = 0
t = 0
h = 0.01
n = 5000
Do from i = 1 to n
  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
End Do

Of course, there are much cleaner approaches in code, but since you are using a spreadsheet, this gives everything in a form like you posted. Here are various plots for each of the parameters with each other for the initial conditions (recall that this system can exhibit chaotic behavior based on the initial conditions).