Double Spring Pendulum

It's a small variation of a simple physics problem, the double pendulum. In this case, the wires are not rigid, but instead, they're springs, therefore, double spring  pendulum.
The problem consists in finding the motion equations of this system. Here's a representation of the system:

I chose to solve this problem with a Lagrangian approach. Of course it'd be possible with Newton, but it would be so much more "painful".
For this reason, I started by writing the Lagrangian of the system. There are 4 degrees of freedom: $\alpha$, $\beta$, $a$ and $b$. With those identified, we can write the kinetic energy followed by the potential energy.

$K = \frac{1}{2} m_{1}v_{1}^2 + \frac{1}{2} m_{2}v_{2}^2$

$\overrightarrow{v_{1}}=\overrightarrow{\dot{a}}+\overrightarrow{\dot{\alpha}a}$
$\overrightarrow{v_{2}}=\overrightarrow{v_{1}}+\overrightarrow{v'}= \overrightarrow{v_{1}} + \left ( \overrightarrow{\dot{b}}+\overrightarrow{\dot{\beta}b} \right ) $

In this case $\overrightarrow{v_{1}}$ and $\overrightarrow{v_{2}}$ stand for the velocity of each mass, and $\overrightarrow{v'}$ for the relative velocity between the second mass with respect to the first.
It will be useful to find the expression for the square of each. The first one is very easy, due to the fact that the vectors are perpendicular. The second one, unfortunately, isn't that straight forward, but I'll spare you the math in between.

$v_{1}^2=\dot{a}^2+\dot{\alpha}^2a^2$

$v_{2}^2=\dot{a}^2+\dot{\alpha}^2a^2+\dot{b}^2+\dot{\beta}^2b^2 + \\
\ \ \ \ \ \ \ + 2\ sin(\alpha - \beta)\left (\dot{a}\dot{\beta}b - \dot{b}\dot{\alpha}a \right ) + 2\ cos(\alpha - \beta)\left (\dot{a}\dot{b} + \dot{\alpha}a\dot{\beta}b \right )$

The potential energy follows as:

$U=m_{1}gy_{1}+m_{2}gy_{2}+\frac{1}{2}k_{1}\Delta l_{1} ^2+\frac{1}{2}k_{2}\Delta l_{2} ^2$

$y_{1}=-a\ cos\ \alpha$

$y_{2}=-(a\ cos\ \alpha + b\ cos\ \beta)$

$\Delta l_{1} = a - l_{1}$

$\Delta l_{2} = b - l_{2}$

Where $l_{1}$ and $l_{2}$ are the springs' free lengths.

Now that we have defined the Lagrangian, we can calculate the Euler–Lagrange equations of motion:

$$\frac{\partial L}{\partial \alpha}-\frac{\mathrm{d} }{\mathrm{d} t}\left ( \frac{\partial L}{\partial{\alpha}'} \right )=0$$

$$\frac{\partial L}{\partial \beta}-\frac{\mathrm{d} }{\mathrm{d} t}\left ( \frac{\partial L}{\partial{\beta}'} \right )=0$$

$$\frac{\partial L}{\partial a}-\frac{\mathrm{d} }{\mathrm{d} t}\left ( \frac{\partial L}{\partial{a}'} \right )=0$$

$$\frac{\partial L}{\partial b}-\frac{\mathrm{d} }{\mathrm{d} t}\left ( \frac{\partial L}{\partial{b}'} \right )=0$$

At this part, I decided to use the tools of Wolfram Mathematica to solve this equations and find the solution of this system. If you want to see those steps, click the spoiler bellow.

 

Otherwise, the result is:
$$\alpha''=-\frac{1}{a} \left (2\,a'\,\alpha' + g\,sin(\alpha)+\frac{k_{2}}{m_{1}}\,(b-l_{2})\,sin(\alpha-\beta) \right )$$
$$\beta''=-\frac{1}{b} \left(\,2\,b'\,\beta' - \frac{k_{1}}{m_{1}}\,(a-l_{1})\,sin(\alpha-\beta)\,\right )$$
$$a''=\alpha'^2a+g\,cos(\alpha)-\frac{k_{1}}{m_{1}}\,(a-l_{1})+\frac{k_{2}}{m_{1}}\,(b-l_{2})\,cos(\alpha-\beta)$$
$$b''=\beta'^2\,b-\frac{k_{2}}{m_{2}}\,(b-l_{2})\,-\frac{k_{2}}{m_{1}}\,(b-l_{2})+\frac{k_{1}}{m_{1}}\,(a-l_{1})\,cos(\alpha-\beta)$$
Having these equations, I chose to make a little simulation in Python to test the results. In the next gif, you can see how it worked (for some arbitrary initial values):

To check the Python code, click the spoiler button below:

import math
from visual import *

def calc_a_2 (alpha_0,alpha_1, beta_0, beta_1, a_0, a_1, b_0, b_1, l1, l2, k1, k2, m1, m2, g):
    return ((k1*l1+g*m1*math.cos(alpha_0)-k2*l2*math.cos(alpha_0-beta_0)+k2*b_0*math.cos(alpha_0-beta_0)+a_0*(-k1+m1*(alpha_1**2)))/m1)
    
def calc_b_2 (alpha_0,alpha_1, beta_0, beta_1, a_0, a_1, b_0, b_1, l1, l2, k1, k2, m1, m2, g):
    return ((k2*l2*m1+k2*l2*m2-k1*l1*m2*math.cos(alpha_0-beta_0)+k1*m2*a_0*math.cos(alpha_0-beta_0)-b_0*(k2*(m1+m2)-m1*m2*(beta_1**2)))/(m1*m2))

def calc_alpha_2 (alpha_0,alpha_1, beta_0, beta_1, a_0, a_1, b_0, b_1, l1, l2, k1, k2, m1, m2, g):
    return (-((g*m1*math.sin(alpha_0)-k2*l2*math.sin(alpha_0-beta_0)+k2*b_0*math.sin(alpha_0-beta_0)+2*m1*a_1*alpha_1)/(m1*a_0)))

def calc_beta_2 (alpha_0,alpha_1, beta_0, beta_1, a_0, a_1, b_0, b_1, l1, l2, k1, k2, m1, m2, g):
    return ((-k1*l1*math.sin(alpha_0-beta_0)+k1*a_0*math.sin(alpha_0-beta_0)-2*m1*b_1*beta_1)/(m1*b_0))

alpha_0 = math.pi/2
beta_0 = 0
alpha_1 = 0
beta_1 = 0
alpha_2 = 0
beta_2 = 0

a_0 = 1
b_0 = 1
a_1 = 0
b_1 = 0
a_2 = 0
b_2 = 0

l1 = 1
l2 = 1
m1 = 1
m2 = 2
k1 = 200
k2 = 300
g = 9.81

dt = 0.001
t_max = 60
t=0


scene.autoscale = 0 # stop it from zooming in and out
scene.title = 'Double pendulum'
scene.range = (3, 3, 1)

#Center
ball0 = sphere(pos=vector(0, 0, 0), radius=0.05, color=color.cyan)

#Balls
ball1 = sphere(pos=vector(a_0*math.sin(alpha_0), -a_0*math.cos(alpha_0),0), radius=0.12, color=color.blue, make_trail=True, retain=int(1.0/dt))
ball2 = sphere(pos=vector(a_0*math.sin(alpha_0) + b_0*math.sin(beta_0), -a_0*math.cos(alpha_0) - b_0*math.cos(beta_0), 0), radius=0.12, color=color.red, make_trail=True, retain=int(1.0/dt))

#Strings
arm1 = cylinder(pos=(0,0,0), axis=(a_0*math.sin(alpha_0), -a_0*math.cos(alpha_0),0), radius=.03, color=color.white)
arm2 = cylinder(pos=(a_0*math.sin(alpha_0), -a_0*math.cos(alpha_0),0), axis=(b_0*math.sin(beta_0), -b_0*math.cos(beta_0), 0), radius=.03, color=color.white)

#while t<t_max :
while 1:
    rate(1/dt)
    
    if(round(t, int(-math.log10(dt)))%1 == 0):
        print "%ds" %t
        
    alpha_2 = calc_alpha_2 (alpha_0,alpha_1, beta_0, beta_1, a_0, a_1, b_0, b_1, l1, l2, k1, k2, m1, m2, g)
    beta_2 = calc_beta_2 (alpha_0,alpha_1, beta_0, beta_1, a_0, a_1, b_0, b_1, l1, l2, k1, k2, m1, m2, g)
    a_2 = calc_a_2 (alpha_0,alpha_1, beta_0, beta_1, a_0, a_1, b_0, b_1, l1, l2, k1, k2, m1, m2, g)
    b_2 = calc_b_2 (alpha_0,alpha_1, beta_0, beta_1, a_0, a_1, b_0, b_1, l1, l2, k1, k2, m1, m2, g)

    alpha_1 = alpha_1 + alpha_2*dt
    beta_1 = beta_1 + beta_2*dt
    a_1 = a_1 + a_2*dt
    b_1 = b_1 + b_2*dt

    alpha_0 = alpha_0 + alpha_1*dt
    beta_0 = beta_0 + beta_1*dt
    a_0 = a_0 + a_1*dt
    b_0 = b_0 + b_1*dt

    alpha_0 = alpha_0%(2*math.pi)
    beta_0 = beta_0%(2*math.pi)
    
    ball1.pos.x = a_0*math.sin(alpha_0)
    ball1.pos.y = -a_0*math.cos(alpha_0)
    arm1.axis = (ball1.pos.x, ball1.pos.y, 0)
    
    ball2.pos.x = ball1.pos.x + b_0*math.sin(beta_0)
    ball2.pos.y = ball1.pos.y - b_0*math.cos(beta_0)
    arm2.pos = (ball1.pos.x, ball1.pos.y, 0)
    arm2.axis = (ball2.pos.x - ball1.pos.x, ball2.pos.y - ball1.pos.y, 0)

#    if (t==0):
#        raw_input()
    
    t=t+dt

Hope that this was helpful to you in any way! If any doubt remains, feel free to contact me or to comment bellow.