Numerical study of coupled oscillator system using the classical Euler-Lagrange equations

In classical mechanics, there are many useful techniques that can be used in solving physical systems. Newtonian mechanics is a technique based on vector concepts (i.e., forces). In this technique, one has to write down all forces acting on the system, and then, applying Newton’s second law. Another useful technique based on scalars quantities (i.e., kinetic and potential energies) has been developed. This technique, known as Lagrangian mechanics, helped physicists in solving a large number of complicated systems, where Newtonian techniques cannot be applied effectively in many areas in physics, such as classical mechanics [1, 2], nonlinear systems, and control theory [3, 4].

To apply Lagrangian mechanics, one has to write down the kinetic and potential energies of the system, and as a result, the Lagrangian L(q̇, q, t) is obtained, where q represents the generalized coordinate, q̇ is the generalized velocity, and t is the time. Applying the following equation ∂L∂q−ddt∂L∂q˙=0 {{\partial L} \over {\partial q}} - {d \over {dt}}{{\partial L} \over {\partial \dot q}} = 0 to the Lagrangian L(q̇, q, t), we obtained a set of differential equations known as Euler–Lagrange equations (ELE) or equations of motion. For more details, one can refer to [1]. The ELE can be solved analytically, and thus, the motion of the system is specified. In many cases, the ELE cannot be solved analytically. Therefore, we have to seek for suitable numerical techniques to apply, and in this case, some initial conditions are needed.

On the other hand, it is important to remember that all mechanical and electrical vibration problems reduce, in the case of small oscillations, to problems involving one or several coupled oscillators. Problems involving vibrations of strings, membranes, elastic solids, electrical, or acoustical can be reduced to problems of coupled oscillators; see for examples [5, 6] and references within them.

Many numerical methods have been introduced in literature to solve differential equations [7, 8, 9, 10, 11, 12]. These numerical methods are based on MATLAB, MATHEMATICA, MATCAD software and other computerized programs. For example, MATLAB contains many built-in codes that can be used in such cases (e.g., ode45). These numerical techniques enable us to solve many systems numerically, even when the analytical solution is not known. In addition, the ode45 code depends on a Runge-Kutta method, which can be used efficiently in solving different kinds of differential equations. One can refer to [13, 14] and references therein to construct a good idea about applying the Runge-Kutta method on solving differential equations.

On the other hand, many ideas have been carried out by applying MATLAB and built-in codes in solving some differential equations. For example in [15, 16, 17], the authors solved some differential equations describing physical systems by applying Lagrangian mechanics and obtaining the ELEs, which are in general differential equations that need to be solved numerically for certain initial conditions.

The outline of this paper is arranged as follows: In Section 2, the physical description of the system is presented, and the Lagrangain of the coupled oscillator is obtained. In Section 3, the numerical method used is presented. Section 4 provides the numerical solutions of the derived equations of motion for different initial conditions. The results were discussed in Section 5. We finally end the paper with conclusions.

Let us consider the system shown in Figure 1, which consists of two coupled carts of masses m₁ and m₂, attached to their respective walls by two identical springs (with force constant k) and to each other by another spring (with force constant k₂) [1]. It is of interest to mention that in some books, this system is known as: diatomic molecules.

Figure 1

Assuming all springs to be massless, we aim to obtain the ELE for the system, and then we require obtaining a numerical solution for these equations under different initial conditions.

First of all, the kinetic energy (T) of the system reads: (1) T=12m1x˙12+12m2x˙22. T = {1 \over 2}{m_1}\dot x_1^2 + {1 \over 2}{m_2}\dot x_2^2.

In order to write the potential energy (V) of the system, we have to identify the extensions of the three springs shown in Figure 1 as (x₁,x₁ − x₂, and −x₂). As a result, the potential energy reads: (2) V=12(k+k2)x12−k2x1x2+12(k2+k)x22. V = {1 \over 2}\left({k + {k_2}} \right)x_1^2 - {k_2}{x_1}{x_2} + {1 \over 2}\left({{k_2} + k} \right)x_2^2.

The classical Lagrangian can be immediately written as: (3) L=T−V=12m1x˙12+12m2x˙22−12(k+k2)x12+k2x1x2−12(k2+k)x22. \matrix{L \hfill & {= T - V = {1 \over 2}{m_1}\dot x_1^2 + {1 \over 2}{m_2}\dot x_2^2 - {1 \over 2}\left({k + {k_2}} \right)x_1^2} \hfill \cr {} \hfill & {+ {k_2}{x_1}{x_2} - {1 \over 2}\left({{k_2} + k} \right)x_2^2.} \hfill \cr}

Our final step in this section is in obtaining the equations of motion. As such, applying our Lagrangian (eq. 3) to the following equation ∂L∂q−ddt∂L∂q˙=0 {{\partial L} \over {\partial q}} - {d \over {dt}}{{\partial L} \over {\partial \dot q}} = 0 , with q = x₁, x₂, we got the following ELE:

For q = x₁, we have: (4) m1x¨1=−(k−k2)x1+k2x2. {m_1}{\ddot x_1} = - \left({k - {k_2}} \right){x_1} + {k_2}{x_2}.

Now, we aim to obtain the Hamilton’s equations of motion. The first point is introducing the following generalized momenta: (6) px1=∂L∂x˙1=mx˙1. {p_{{x_1}}} = {{\partial L} \over {\partial {{\dot x}_1}}} = m{\dot x_1}. (7) px2=∂L∂x˙2=mx˙2. {p_{{x_2}}} = {{\partial L} \over {\partial {{\dot x}_2}}} = m{\dot x_2}.

As a result, the Hamiltonian H = ∑ q̇_ip_xi − L, with I = 1, 2 reads: (8) H=px122m1+px222m1+12(k+k2)x12−k2x1x2+12(k2+k)x22. \matrix{H \hfill & {= {{p_{{x_1}}^2} \over {2{m_1}}} + {{p_{{x_2}}^2} \over {2{m_1}}} + {1 \over 2}\left({k + {k_2}} \right)x_1^2 - {k_2}{x_1}{x_2}} \hfill \cr {} \hfill & {+ {1 \over 2}\left({{k_2} + k} \right)x_2^2.} \hfill \cr}

Finally, the Hamilton’s equations can be obtained as: (9) x˙1=∂H∂Px1⇒m1x¨1=−(k+k2)x1+k2x2. {\dot x_1} = {{\partial H} \over {\partial {P_{{x_1}}}}} \Rightarrow {m_1}{\ddot x_1} = - \left({k + {k_2}} \right){x_1} + {k_2}{x_2}. (10) x˙2=∂H∂Px2⇒m2x¨2=k2x1−(k2+k)x2. {\dot x_2} = {{\partial H} \over {\partial {P_{{x_2}}}}} \Rightarrow {m_2}{\ddot x_2} = {k_2}{x_1} - \left({{k_2} + k} \right){x_2}.

In the section below, we aimed to obtain numerical solution for eqs. (9, 10) for some initial conditions.

A standard solver for ordinary differential equations (ODEs) in the simulation of MATLAB’s framework is the ode45 function. ODE45 executes the Runge-Kutta method with a variable time step for optimal computation [18, 19, 20]. A fourth-order Runge-Kutta algorithm is useful in solving non-linear differential equations, taking into consideration some given initial conditions [21]. Any problem in the following manner (11) dxdt=f(t,x). {{dx} \over {dt}} = f\left({t,x} \right).

With x(t_o) = x_o

Where t is the independent variable, x is a vector of dependent variables to be found, and f (t, x) is a function of t, and x can be solved using ode45. Specifying f (t, x) and the initial conditions enable us to solve the problem [22].

Bellow a short explanation is presented to explain how Runge-Kutta works [23].

In general, given that: (12) y'=f(x,y). y' = f\left({x,y} \right). (13) y(xn)=yn. y\left({{x_n}} \right) = {y_n}.

Then, one can compute in turn: (14) k1=hf(xn,yn). {k_1} = hf\left({{x_n},{y_n}} \right). (15) k2=hf(xn+h2,yn+k12). {k_2} = hf\left({{x_n} + {h \over 2},{y_n} + {{{k_1}} \over 2}} \right). (16) k3=hf(xn+h2,yn+k22). {k_3} = hf\left({{x_n} + {h \over 2},{y_n} + {{{k_2}} \over 2}} \right). (17) k4=hf(xn+h,yn+k3). {k_4} = hf\left({{x_n} + h,{y_n} + {k_3}} \right). (18) yn+1=yn+16(k1+2k2+2k3+k4). {y_{n + 1}} = {y_n} + {1 \over 6}\left({{k_1} + 2{k_2} + 2{k_3} + {k_4}} \right).

Now, everything becomes in first-order form, so the list of commands is needed in the code [24].

ODE45 is used to implement “first try” solution for most problems [25, 26]. For more details about the Runge-Kutta method, one can refer to [23, 27].

Repetition and conditional structures are used to repeat sequences in MATLAB, and from these structures are: if and if-else statements, do while, for each, while, and for loops.

In literature, there are many examples of repetition and conditional structures [23, 27].

In Section 4, we propose the numerical solution for eqs. (4 and 5) using loop structure. We consider some initial conditions, and the numerical solution for these particular conditions has been obtained.

Below, we consider some selected cases to obtain the numerical solution of eqs. (4, 5):

Case 1: Symmetric Case

In this case, we considered the two initial conditions x₁(0) = −x₂(0) = 1 and ẋ₁(0) = − ẋ₂(0) = 0.5; on the other hand, the other parameters were assumed to be k = 1, 2, with k₂ = 2k, while m₁ = 3 and m₂ = 7, see Figures 2, 3:

Figure 2

Figure 3

Case 2: Anti-symmetric Case

In this case, we considered the following two initial conditions x₁(0) = x₂(0) = 1 and ẋ₁(0) = ẋ₂(0) = 0.5, and the other parameters were assumed to be k = 1, and 2, with k₂ = 2k, while m₁ = 3 and m₂ = 7, see Figures 4, 5:

Figure 4

Figure 5

Case 3: Weak Coupled Oscillator (k₂ ≺≺ k)

Finally, we considered the weak coupled oscillator case, which is of special importance. In this case, the following parameters are considered k = 1, and 2, with k₂ = 0.01, while m₁ = 3 and m₂ = 7, with the following two initial conditions x₁(0) = 1, x₂(0) = 0 and ẋ₁(0) = ẋ₂(0) = 0, see Figures 6, 7:

Figure 6

Figure 7

In this paper, the motion of the coupled oscillator system has been studied numerically using the ode45 code, and the solutions x₁(t) and x₂(t) have been plotted for the several considered cases as in Figures 2–7.

Figures 2, 3 present the positions x₁(t) and x₂(t) of masses m₁ and m₂ against the time, respectively, for the symmetric case. The results in the same figures show that the periodicity of the oscillations has disappeared. Also, the influence on mass m₂ is clearer than that on m₁, because it has a larger value. In addition, because m₁ and m₂ were displaced initially in opposite directions (i.e., x₁(0) = −x₂(0)) and then released, this is the main reason for disappearing the periodicity of the oscillations.

Figures 4, 5 present the positions x₁(t) and x₂(t) of masses 1 and 2 against the time, respectively, for the antisymmetric case. Again, the periodicity of the oscillation has disappeared, but less than that in Figures 2, 3. The reason here refers to the fact that the two masses m₁ and m₂ were displaced initially in the same direction (i.e., x₁(0) = x₂(0)) and then released.

For Cases 1 and 2, it is expected that the periodicity of the oscillations will be destroyed. This is due to the fact that the spring-connected the two masses together have a larger stiffness than the two other springs.

In Figures 6, 7, the positions x₁(t) and x₂(t) of masses 1 and 2, respectively, is displaced against time for the weak coupled oscillator, where the stiffness of the spring connecting the two masses is very small compared with the other springs stiffness. In our case, we consider k₂ = 0.01, whereas k = 1 is shown in Figure 6, and k = 2 is shown in Figure 7. It is clear from these two figures that mass m₁ undergoes a simple harmonic oscillations, with a constant amplitude and a constant period, while mass m₂ nearly vibrates but with very small oscillations around the equilibrium position; this is due to the fact that m₂ has a larger value than m₁.

Finally, one can see from the figures that if the stiffness of the spring connecting the two masses is larger than the other springs’ stiffness, then this will not allow the system to undergo simple harmonic motion. As a result, the periodicity of the motion is destroyed.

In conclusion, the periodicity of the oscillator disappeared, especially in the symmetric case. While in the antisymmetric case, the periodicity is nearly still there and not destroyed. In the weak coupled case, one can see from Figures 6, 7 that mass m₂ does not move, because it is larger than m₁.

After releasing the system, the displacement in the symmetric and anti-symmetric cases exceeds the initial displacement. This can be explained by the effects of springs on each other. Nonetheless, for the weak coupled oscillator, mass m₁ undergoes a simple harmonic motion and not affected by mass m₂.