To make it easier to obtain the equation of motions, the top mass \(m\) is modeled as attached to spring of stiffness \(k\) which is in turn attached to an infinitely stiff vertical massless beam. This way the vibration of the mass \(m~\)at the top can be more easily modeled.
Based on the above diagram, we now obtain the free body diagram as follows. In this, we assume that \(x_{2}>x_{1}\) and both as positive. Hence spring \(k\) attached to \(m~\) is in tension.
The top mass \(m\) vibrates in horizontal direction only. Hence this assumes the spring will remain horizontal and we must assume that \(x_{2}-x_{1}\) remain small for this model to be realistic.
From this free body diagram we see now that the reaction force \(F_{x}\) is equal to \(k\left ( x_{2}-x_{1}\right )\). (By resolving forces in the \(x\) direction for the massless beam).
Therefore \[ F_x = k\left ( x_{2}-x_{1}\right ) \] And the equation of motion for \(x_{2}\) is\begin{align} m\ddot{x}_{2} & =-k\left ( x_{2}-x_{1}\right ) \nonumber \\ m\ddot{x}_{2}+kx_{2}-kx_{1} & =0\tag{1} \end{align}
The equation of motion for the cart is\begin{align} 2m\ddot{x}_{1} & =-4kx_{1}+F_{x}\nonumber \\ 2m\ddot{x}_{1} & =-4kx_{1}+k\left ( x_{2}-x_{1}\right ) \nonumber \\ 2m\ddot{x}_{1}+5kx_{1}-kx_{2} & =0\tag{2} \end{align}
Writing (1) and (2) in matrix form\[ \fbox{$\begin{bmatrix} 2m & 0\\ 0 & m \end{bmatrix}\begin{Bmatrix} \ddot{x}_{1}\\ \ddot{x}_{2}\end{Bmatrix} +\begin{bmatrix} 5k & -k\\ -k & k \end{bmatrix}\begin{Bmatrix} x_{1}\\ x_{2}\end{Bmatrix} =\begin{Bmatrix} 0\\ 0 \end{Bmatrix} $}\] Or\[\begin{bmatrix} 4 & 0\\ 0 & 2 \end{bmatrix}\begin{Bmatrix} \ddot{x}_{1}\\ \ddot{x}_{2}\end{Bmatrix} +\begin{bmatrix} 1000 & -200\\ -200 & 200 \end{bmatrix}\begin{Bmatrix} x_{1}\\ x_{2}\end{Bmatrix} =\begin{Bmatrix} 0\\ 0 \end{Bmatrix} \] The first step is to find the eigenvalues (which are the square of the natural frequency) for the system.
Let \begin{align*} A & =M^{-1}K\\ & =\begin{bmatrix} 4 & 0\\ 0 & 2 \end{bmatrix} ^{-1}\begin{bmatrix} 1000 & -200\\ -200 & 200 \end{bmatrix} \end{align*}
But \begin{align*} \begin{bmatrix} 4 & 0\\ 0 & 2 \end{bmatrix} ^{-1} & =\frac{1}{\det \left ( M\right ) }\begin{bmatrix} 2 & 0\\ 0 & 4 \end{bmatrix} \\ & =\frac{1}{8}\begin{bmatrix} 2 & 0\\ 0 & 4 \end{bmatrix} \\ & =\begin{bmatrix} \frac{1}{4} & 0\\ 0 & \frac{1}{2}\end{bmatrix} \end{align*}
Hence\begin{align*} A & =\begin{bmatrix} \frac{1}{4} & 0\\ 0 & \frac{1}{2}\end{bmatrix}\begin{bmatrix} 1000 & -200\\ -200 & 200 \end{bmatrix} \\ & =\begin{bmatrix} 250 & -50\\ -100 & 100 \end{bmatrix} \end{align*}
Now we will find the eigenvalues of \(A\) (these will be the \(\omega _{n}^{2}\) values). To find the eigenvalues of \(A\), we solve \begin{align*} \det \left ( \left [ A\right ] -\lambda \left [ I\right ] \right ) & =0\\ \det \left ( \begin{bmatrix} 250 & -50\\ -100 & 100 \end{bmatrix} -\begin{bmatrix} \lambda & 0\\ 0 & \lambda \end{bmatrix} \right ) & =\\\begin{vmatrix} 250-\lambda & -50\\ -100 & 100-\lambda \end{vmatrix} & =\\ \left ( 250-\lambda \right ) \left ( 100-\lambda \right ) -5000 & =0\\ \lambda ^{2}-350\lambda +20\,000 & =0 \end{align*}
Hence\begin{align*} \lambda & =\frac{-b}{2a}\pm \frac{\sqrt{b^{2}-4ac}}{2a}\\ & =\frac{350}{2}\pm \frac{\sqrt{350^{2}-4\left ( 20\,000\right ) }}{2}\\ & =175\pm 103.08\\ & =\left \{ 71.92,278.08\right \} \end{align*}
Therefore, the eigenvalues are \begin{equation} \lambda =\omega _{n}^{2}=\left \{ 71.92,278.08\right \} \tag{3} \end{equation} The natural frequencies of the system are the sqrt of the eigenvalues. Therefore\begin{align*} \omega _{n} & =\left \{ \sqrt{71.92},\sqrt{278.08}\right \} \\ & =\left \{ 8.4806,16.676\right \} \end{align*}
Hence\begin{align*} \omega _{n\left ( 1\right ) } & =8.4806\text{ rad/sec}\\ \omega _{n\left ( 2\right ) } & =16.676\text{ rad/sec} \end{align*}
The next step is to find the eigenvectors. These are also called the shape vectors, or the \(u\) vectors. Each eigenvalue will generate one eigenvector. We need to solve\[ \left [ A\right ] \left \{ u\right \} =\lambda \left \{ u\right \} \] For each eigenvalue, we find the corresponding eigenvector.
For \(\lambda =71.92\), we obtain the equation\[\begin{bmatrix} 250 & -50\\ -100 & 100 \end{bmatrix}\begin{Bmatrix} u_{11}\\ u_{21}\end{Bmatrix} =71.92\begin{Bmatrix} u_{11}\\ u_{21}\end{Bmatrix} \] From first equation\[ 250u_{11}-50u_{21}=71.92u_{11}\] We always let \(u_{11}=1\). Therefore\begin{align*} 250-50u_{21} & =71.92\\ u_{21} & =\frac{250-71.92}{50}\\ & =3.5616 \end{align*}
Therefore, the first eigenvector is\[ \vec{u}_{1}=\begin{Bmatrix} 1\\ 3.5616 \end{Bmatrix} \] For \(\lambda =278.\allowbreak 08\), we obtain the equation\[\begin{bmatrix} 250 & -50\\ -100 & 100 \end{bmatrix}\begin{Bmatrix} u_{12}\\ u_{22}\end{Bmatrix} =278.08\begin{Bmatrix} u_{12}\\ u_{22}\end{Bmatrix} \] From first equation\[ 250u_{12}-50u_{22}=278.08u_{12}\] We always let \(u_{12}=1\). Hence\begin{align*} 250-50u_{22} & =278.08\\ u_{22} & =\frac{250-278.08}{50}\\ & =-0.561\,6 \end{align*}
Therefore, the second eigenvector is\[ \vec{u}_{2}=\begin{Bmatrix} 1\\ -0.561\,6 \end{Bmatrix} \] Therefore the modal matrix \(\left [ u\right ] \) is\[ \fbox{$u=\begin{bmatrix} 1 & 1\\ 3.5616 & -0.5616 \end{bmatrix} $}\] And \(\Omega \) matrix is\begin{align*} \Omega & =\begin{bmatrix} \omega _{n\left ( 1\right ) }^{2} & 0\\ 0 & \omega _{n\left ( 2\right ) }^{2}\end{bmatrix} \\ & =\begin{bmatrix} 71.92 & 0\\ 0 & 278.08 \end{bmatrix} \end{align*}
And the system of equations written in principle coordinates \(q\) is\begin{align*} \left \{ \ddot{q}\right \} +\left [ \Omega \right ] \left \{ q\right \} & =\left \{ 0\right \} \\\begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}\begin{Bmatrix} \ddot{q}_{1}\left ( t\right ) \\ \ddot{q}_{2}\left ( t\right ) \end{Bmatrix} +\begin{bmatrix} 71.92 & 0\\ 0 & 278.08 \end{bmatrix}\begin{Bmatrix} \ddot{q}_{1}\left ( t\right ) \\ \ddot{q}_{2}\left ( t\right ) \end{Bmatrix} & =\begin{Bmatrix} 0\\ 0 \end{Bmatrix} \end{align*}
which is now decoupled. The solution in normal coordinates is\begin{align*} \begin{Bmatrix} x_{1}\left ( t\right ) \\ x_{2}\left ( t\right ) \end{Bmatrix} & =A_{1}\begin{Bmatrix} u_{11}\\ u_{21}\end{Bmatrix} \cos \left ( \omega _{n\left ( 1\right ) }t-\phi _{1}\right ) +A_{2}\begin{Bmatrix} u_{12}\\ u_{22}\end{Bmatrix} \cos \left ( \omega _{n\left ( 2\right ) }t-\phi _{2}\right ) \\ & =A_{1}\begin{Bmatrix} 1\\ 3.5616 \end{Bmatrix} \cos \left ( 8.481t-\phi _{1}\right ) +A_{2}\begin{Bmatrix} 1\\ -0.561\,6 \end{Bmatrix} \cos \left ( 16.676t-\phi _{2}\right ) \end{align*}
This is derivation of the same equations of motions using energy method. (In this example, this method is much simpler to use to find equation of motions). The kinetic energy of the system is\[ T=\frac{1}{2}m\dot{x}_{2}^{2}+\frac{1}{2}\left ( 2m\right ) \dot{x}_{1}^{2}\] And the potential energy comes only from the springs, since we assumed the top mass \(m\) remain horizontal as it vibrates back and forth\[ U=\frac{1}{2}4kx_{1}^{2}+\frac{1}{2}k\left ( x_{2}-x_{1}\right ) ^{2}\] Therefore the Lagrangian is\begin{align*} \Gamma & =T-U\\ & =\frac{1}{2}m\dot{x}_{2}^{2}+m\dot{x}_{1}^{2}-\frac{1}{2}\left ( 4k\right ) x_{1}^{2}-\frac{1}{2}k\left ( x_{2}-x_{1}\right ) ^{2} \end{align*}
EQM for \(x_{1}\)\begin{align} \frac{d}{dt}\left ( \frac{\partial \Gamma }{\dot{x}_{1}}\right ) -\frac{\partial \Gamma }{x_{1}} & =0\nonumber \\ \frac{d}{dt}\left ( 2m\dot{x}_{1}\right ) -\left ( -4kx_{1}+k\left ( x_{2}-x_{1}\right ) \right ) & =0\nonumber \\ 2m\ddot{x}_{1}-\left ( -4kx_{1}+kx_{2}-kx_{1}\right ) & =0\nonumber \\ 2m\ddot{x}_{1}-\left ( -5kx_{1}+kx_{2}\right ) & =0\nonumber \\ 2m\ddot{x}_{1}+5kx_{1}-kx_{2} & =0\tag{1} \end{align}
EQM for \(x_{2}\)\begin{align} \frac{d}{dt}\left ( \frac{\partial \Gamma }{\dot{x}_{2}}\right ) -\frac{\partial \Gamma }{x_{2}} & =0\nonumber \\ \frac{d}{dt}\left ( m\dot{x}_{2}\right ) -\left ( -k\left ( x_{2}-x_{1}\right ) \right ) & =0\nonumber \\ m\ddot{x}_{2}-\left ( -kx_{2}+kx_{1}\right ) & =0\nonumber \\ m\ddot{x}_{2}+kx_{2}-kx_{1} & =0\tag{2} \end{align}
In Matrix form (1,2) becomes\[\begin{bmatrix} 2m & 0\\ 0 & m \end{bmatrix}\begin{Bmatrix} \ddot{x}_{1}\\ \ddot{x}_{2}\end{Bmatrix} +\begin{bmatrix} 5k & -k\\ -k & k \end{bmatrix}\begin{Bmatrix} x_{1}\\ x_{2}\end{Bmatrix} =\begin{Bmatrix} 0\\ 0 \end{Bmatrix} \] Which is the same exact result obtained earlier.
The Matlab code is the following
The output is
Definitions For stiffness matrix \(\left [ K\right ] \), element \(k_{ij}\) means: Apply unit displacement at location \(j\) and measure the force at location \(i\). While for flexibility matrix \(\left [ a\right ] \), its element \(a_{ij}\) means: Apply unit force at location \(j\) and measure the displacement at location \(i\).
To solve this problem, this part of handout is used
Since \(\left [ a\right ] \) is symmetric, only lower triangle part needs to be found (or upper triangle).
\[\begin{bmatrix} a_{11} & & \\ a_{21} & a_{22} & \\ a_{31} & a_{32} & a_{33}\end{bmatrix} \] To find \(a_{11}\), a unit force is put at location \(m_{1}\) and displacement at \(m_{1}\) is measured. To find \(a_{21}\), a unit force is put at location \(m_{1}\) and displacement at \(m_{2}\) is measured and so on. The formulas in the above hand out are used for this. To speed this process and make less error, a small function is written to do the computation. Here is the function and the result generated for \(a_{11},a_{21},a_{32,}a_{22,}a_{32,}a_{33}\)
Therefore, using this result, the lower triangle is\[\begin{bmatrix} \frac{9}{64} & & \\ \frac{1}{6} & \frac{1}{3} & \\ \frac{13}{192} & \frac{1}{6} & \frac{9}{64}\end{bmatrix} \frac{L^{3}}{EI}\] Hence by symmetry\[ \left [ a\right ] =\begin{bmatrix} \frac{9}{64} & \frac{1}{6} & \frac{13}{192}\\ \frac{1}{6} & \frac{1}{3} & \frac{1}{6}\\ \frac{13}{192} & \frac{1}{6} & \frac{9}{64}\end{bmatrix} \frac{L^{3}}{EI}\]