The dynamics of coupled bodies and oscillators is significant in mechanics, engineering, electronics as well as biological systems. They are mostly represented as nonlinear dynamical systems  . One of the most important stages in the analysis of any mechanical model is to establish and find the solution of the dynamical equations which are referred to as equations of motion  . The equations of motion are often derived by the Euler-Lagrange equations. The fundamental idea of the Lagrangean approach to mechanics is to reformulate the equations of motion in terms of the dynamical variables that describe the degree of freedom, and thereby incorporate constraint forces into the definition of the degrees of freedom rather than explicitly including them as forces in Newton’s second law. Of importance is the notion of stability of a given dynamical system, where we would be concerned with the stability of some critical point of the system. Indeed stability plays a central role in system engineering, especially in the field of control system and automation, with regards to both dynamics and control.
Chutiphon  suggested Lyapunov stability as a general and useful approach to analyze the stability of nonlinear systems. It has two approaches: indirect and direct methods. For the second method of Lyapunov (indirect method), the idea of system linearization around a given point is used and local stability within small stability regions is possibly achieved. Seyrania and Wang  studied the stability of periodic solutions of the harmonically excited Duffing’s equation with the direct application of the Lyapunov theorem. The damping coefficient and excitation amplitude are assumed to be small. The approximate methods were used to find the periodic solutions. They derived the stability conditions and found stable and unstable region on the frequency response curve.
Maliki and Nwoba  studied a mathematical model of a coupled system of harmonic oscillators using the generalized coordinates and Euler-Lagrange equation. Laplace transform was also used to get the analytical solution of the system. The stability analysis of the system was investigated by the direct method and it was observed that the coupled system is asymptotically stable for the strictly negative roots and strongly unstable for the positive roots.
Maliki and Okereke  investigated the stability analysis of certain third order linear and nonlinear ordinary differential equations. They employed the method of phase portrait analysis and showed, using simulation that the Hartman-Groβman theorem is verified, for a second order linearized system, which approximates the nonlinear system, preserving the topological features.
1.1. Statement of the Problem
We consider the problem of analyzing the dynamics of a triple pendulum as shown in Figure 1. Despite the complexity of the system we will obtain, with the help of the Euler Lagrange (E-L) equations, the equation of motion of each of the individual masses (assumed to be equal). The system is coupled and joined by three rigid rods of equal lengths, hinged to each of the masses.
1.2. Model Formulation
From the figure below we choose
as our generalized coordinates.
To compute the Lagrangean of the system, we first compute the total kinetic energy (K.E) and potential energy (P.E) of the system. However, we require the following expressions.
Differentiating the above coordinates with respect to time, we get;
The square of the velocities for each mass is given by;
After some algebraic manipulations, we get
The total kinetic energy of the pendulum is then:
The total potential energy of the pendulum is the sum of the potential energy of each mass;
Recall that the Lagrangean is given by
We now employ the E-L equations to obtain the equations of motion, i.e.;
For the first generalized coordinate j = 1, we have;
For the second generalized coordinate j = 2 in (1.20), we have;
For the third generalized coordinate j = 3 in (20), we have;
Therefore, the equations of motion are;
Assuming equal masses and equal lengths of rods, i.e.,
The equations become respectively;
We assume the generalized coordinates
representing the angular displacements are small so that;
The equations of motion for the coupled rigid body then become;
Following the above procedure, we have for Equation (28)
Similarly for Equation (29) we have;
The model equations for the given problem are summarized as;
1.3. Stability Analysis
For the purpose of stability analysis we must vectorize the above coupled system of differential equations.
Thus the vectorized system of equations for the coupled pendulums is;
We shall investigate the stability of the pendulum at the critical point, where
. This implies that
and by substitution our critical point is
Equations (34)-(39) can, for convenience) be written simply as;
represent the RHS of the system.
The Jacobian of the system is written
Computing the entries and evaluating at the critical point, we get
The eigenvalues of the matrix J are given by
This simplifies to give;
. For simplicity we take
, and naturally
We now employ the POLYROOT algorithm in MathCad  to solve the polynomial Equation (42).
Define a vector v of the coefficients beginning with the constant term, i.e.:
Clearly not all the eigenvalues have negative real parts, we therefore conclude that the critical point of the system is unstable.
2. Simulation of the System of ODEs
In MathCAD we define the vector of derivatives, viz.;
Additional arguments for the ODE solver are
Initial value of independent variable
End value of independent variable
Vector of initial function values
Number of solution values on [t0, t1]
Independent variable values
First solution function values
Second solution function values
Third solution function values
Fourth solution function values
Fifth solution function values
Sixth solution function values
Table 1 is the MathCad solution matrix for the system of differential equations describing the dynamics of the coupled pendulums. This solution matrix is obtained using the Runge-Kutta algorithm. Figures 2(a)-(c) depict the graphical profiles of the solution curves. We recall that the variables
are actually the generalized coordinates, representing the angular displacements from the vertical position, while their derivates with respect to time
are the angular velocities.
The graph of
(Figure 2(a)) starts from the origin, then performs irregular oscillations with increasing amplitude over time. Figure 2(b) and Figure 2(c) depict a similar variation for
Figure 2(g) is the profile of the combined graph of
against time. This is important because although the coupled pendulum performs irregular oscillations, at specific times the oscillations for each mass intersect, meaning they all pass through the same point. Furthermore, we observe that the respective amplitudes for each generalized coordinate is increasing. Figures 2(d)-(f) depict respectively the variation of the angular velocities
Table 1. Solution matrix for the ODE.
(a) (b) (c)
Figure 3. (a) Graph of u4 against u1; (b) Graph of u5 against u2; (c) Graph of u6 against u3.
. These velocities are clearly oscillatory and irregular with increasing amplitude.
Figures 3(a)-(c) are respectively the phase portrait for
. For each of the phase portraits we obtain a spiral source.
Coupled oscillators in general, are useful in the study of vibrations and as such computing the specific modes of vibration for oscillating systems is very important from a practical point of view, particularly in engineering. The study of coupled systems is useful in mechanics, electronics as well as biological systems. It is also pertinent for the biomechanical analysis of animals, humans or robotic systems. An early application of mathematical modeling in biological systems was first pioneered by Van der Pol in 1928. The Van der Pol oscillator  was used to explain some normal and pathological rhythms of the heart.
Our research work has clearly demonstrated the power of mathematical modeling, where the behaviour of a complex mechanical system can be captured in the form of a system of ordinary differential equations. Furthermore, with versatile mathematical software such as MathCAD, the system can be analyzed in detailed graphical format to give a deep understanding of the inherent dynamics as well as its stability properties.