Finding the stability from the eigenvalues of two second-order DEs

I'm analysing the stability of a double physical pendulum, and have determined the dimensionless system of equations to be \begin{align} \label{eq:dimensionless} \begin{aligned} & \ddot{\phi} \cos{(\theta - \phi)} + \dot{\phi}^2 \sin{(\theta - \phi)} + \alpha \ddot{\theta} + \beta \sin{\theta} = 0 \end{aligned} \\[5pt] \begin{aligned} & \ddot{\theta} \cos{(\theta - \phi)} - \dot{\theta}^2 \sin{(\theta - \phi)} + \gamma \ddot{\phi} + \sin{\phi} = 0 \end{aligned} \end{align} where $\alpha$, $\beta$ and $\gamma$ are defined to be \begin{align} \label{eq:dim_parameters} \alpha &= \frac{\big ( I_1 + (m_1 + m_2) L_{1}^2 \big )}{m_2 L_{1} L_{2}} \\[5pt] \beta &= \frac{L_1 (m_1 + m_2)}{m_2 L_{2}} \\[5pt] \gamma &= \frac{(I_2 + m_2 L_{2}^2)}{m_2 L_{1} L_{2}} \, . \end{align} Setting $\ddot{\phi} = \ddot{\theta} = \dot{\phi} = \dot{\theta} = 0$, we find that the stationary points of the system are at $(0, 0), \, (0, \pi), \, (\pi, 0), \, (\pi, \pi)$.

I want to find the stability of the system at these points. Here's how I attempted to find them, but am slightly unsure if my working is correct:

I linearised the system of equations \begin{align} \begin{aligned} & \ddot{\phi} + \alpha \ddot{\theta} + \beta \theta = 0 \\[5pt] & \ddot{\theta} + \phi + \gamma \ddot{\phi} = 0 \, . \end{aligned} \end{align} Rearranging and substituting the equations into each other \begin{align} \begin{aligned} \nonumber & \ddot{\theta} = \frac{\beta \gamma}{1 - \alpha \gamma} \theta - \frac{1}{1 - \alpha \gamma} \phi \end{aligned} \\[5pt] \begin{aligned} \nonumber & \ddot{\phi} = \frac{\alpha}{1 - \alpha \gamma} \phi - \frac{\beta}{1 - \alpha \gamma} \theta \, . \end{aligned} \end{align} Expressing this as a matrix equation \begin{align} \begin{aligned} \label{eq:matrix} \frac{d}{dt} \begin{bmatrix} \theta \\ \dot{\theta} \\ \phi \\ \dot{\phi} \end{bmatrix} &= \begin{bmatrix} 0 & 1 & 0 & 0 \\ \frac{\beta \gamma}{1 - \alpha \gamma} & 0 & -\frac{1}{1 - \alpha \gamma} & 0 \\ 0 & 0 & 0 & 1 \\ -\frac{\beta}{1 - \alpha \gamma} & 0 & \frac{\alpha}{1 - \alpha \gamma} & 0 \end{bmatrix} \begin{bmatrix} \theta \\ \dot{\theta} \\ \phi \\ \dot{\phi} \end{bmatrix} \, . \end{aligned} \end{align} Finding the characteristic equation for the eigenvalues of the matrix \begin{align} \begin{aligned} \nonumber & \text{det}(A - \lambda I) = \begin{vmatrix} -\lambda & 1 & 0 & 0 \\ \frac{\beta \gamma}{1 - \alpha \gamma} & -\lambda & -\frac{1}{1 - \alpha \gamma} & 0 \\ 0 & 0 & -\lambda & 1 \\ -\frac{\beta}{1 - \alpha \gamma} & 0 & \frac{\alpha}{1 - \alpha \gamma} & -\lambda \end{vmatrix} = 0 \\[5pt] & \implies \lambda^4 - \frac{\alpha + \beta \gamma}{1 - \alpha \gamma} \lambda^2 + \frac{\alpha \beta \gamma - \beta}{(1 - \alpha \gamma)^2} = 0 \\[5pt] & \implies \lambda^4 - \frac{\alpha + \beta \gamma}{1 - \alpha \gamma} \lambda^2 - \frac{\beta}{1 - \alpha \gamma} = 0 \end{aligned} \end{align} Solving for the eigenvalues using the quadratic formula gives \begin{align} \begin{aligned} \lambda^2 &= \frac{(\alpha + \beta \gamma) \pm \sqrt{(-(\alpha + \beta \gamma))^2 - 4 (1 - \alpha \gamma) (-\beta)}}{2 (1 - \alpha \gamma)} \\[5pt] &= \frac{(\alpha + \beta \gamma) \pm \sqrt{\alpha^2 - 2 \alpha \beta \gamma + \beta^2 \gamma^2 + 4 \beta}}{2 (1 - \alpha \gamma)} \\[5pt] &= \frac{(\alpha + \beta \gamma) \pm \sqrt{(\alpha - \beta \gamma)^2 + 4 \beta}}{2 (1 - \alpha \gamma)} \end{aligned} \end{align}

I'm not sure about this working. The stability at these points should not depend on the parameters, surely. I feel that $(0, 0)$ should always be stable.

Have I not found the Jacobian correctly? Or am I, perhaps, supposed to use the fact that the parameters are positive? Any help or clarification would be much appreciated.

• Why do you use this complicated approach? Write instead the potential energy of the system and study its Hessian matrix at the equilibrium points. That is equivalent to exploit Liapunov's theorem rather than linearization procedure... Apr 2, 2017 at 6:10

After continually trying to get this method to work, I still can't seem to get the correct result through my method. I may have made a mistake in my working somewhere, or am incorrectly analysing the signs in the different parameter regimes.

Regardless, I took Valter's advice and analysed the second derivative of the system's potential energy with respect to small displacements around the fixed points. This was much easier to analyse, to say the least.