Multi-Degree of Freedom Systems - Examples 1-5

Example 1 - Derive the differential equation of motion for the system shown in the following figure.


$$\begin{align} & {{m}_{1}}{{{\ddot{x}}}_{1}}=-k{{x}_{1}}-5k\left( {{x}_{1}}-{{x}_{3}} \right)-k\left( {{x}_{1}}-{{x}_{2}} \right)+{{F}_{1}}\left( t \right) \\ & {{m}_{2}}{{{\ddot{x}}}_{2}}=-k\left( {{x}_{2}}-{{x}_{1}} \right)-k\left( {{x}_{2}}-{{x}_{3}} \right)+{{F}_{2}}\left( t \right) \\ & \underline{{{m}_{3}}{{{\ddot{x}}}_{3}}=-5k\left( {{x}_{3}}-{{x}_{1}} \right)-k\left( {{x}_{3}}-{{x}_{2}} \right)+{{F}_{3}}\left( t \right)} \\ & {{m}_{1}}{{{\ddot{x}}}_{1}}+\left( 7{{x}_{1}}-{{x}_{2}}-5{{x}_{3}} \right)k={{F}_{1}}\left( t \right) \\ & {{m}_{2}}{{{\ddot{x}}}_{2}}+\left( -{{x}_{1}}+2{{x}_{2}}-{{x}_{3}} \right)k={{F}_{2}}\left( t \right) \\ & {{m}_{3}}{{{\ddot{x}}}_{3}}+\left( -5{{x}_{1}}-1{{x}_{2}}+7{{x}_{3}} \right)k={{F}_{3}}\left( t \right) \\ \end{align}$$

Or in matrix form
$$\left[ \begin{matrix} {{m}_{1}} & 0 & 0 \\ 0 & {{m}_{2}} & 0 \\ 0 & 0 & {{m}_{3}} \\ \end{matrix} \right]\left( \begin{matrix} {{{\ddot{x}}}_{1}} \\ {{{\ddot{x}}}_{2}} \\ {{{\ddot{x}}}_{3}} \\ \end{matrix} \right)+k\left[ \begin{matrix} 7 & -1 & -5 \\ -1 & 2 & -1 \\ -5 & -1 & 7 \\ \end{matrix} \right]\left( \begin{matrix} {{x}_{1}} \\ {{x}_{2}} \\ {{x}_{3}} \\ \end{matrix} \right)=\left( \begin{matrix} {{F}_{1}}\left( t \right) \\ {{F}_{2}}\left( t \right) \\ {{F}_{3}}\left( t \right) \\ \end{matrix} \right)$$


Runge-Kutta Method for SDOF Systems

In this method the approximate formula used for obtaining x_i+1 from x_i is made to coincide with Taylor’s series expansion of x and x_i+1 up to terms of order (∆t)^n. The Taylor’s series expansion of x(t) at t+∆t is given by:

$$x\left( t+\Delta t \right)=x\left( t \right)+\dot{x}\left( t \right)\Delta t+\ddot{x}\left( t \right)\frac{{{\left( \Delta t \right)}^{2}}}{2!}+\left( t \right)\frac{{{\left( \Delta t \right)}^{3}}}{3!}+\left( t \right)\frac{{{\left( \Delta t \right)}^{4}}}{4!}+....$$


The previous equation requires higher-order derivatives while Runge-Kutta requires only first derivative and not higher than that. In order to solve the equation

$$m\ddot{x}+c\dot{x}+kx=F$$


we need to reduce it to two first-order equations. But before we start to implement R-K method let’s rearrange previous second-order differential equation:

$$\begin{align} & m\ddot{x}=F-c\dot{x}-kx \\ & \ddot{x}=\frac{1}{m}\left[ F-c\dot{x}-kx \right]=f\left( x,\dot{x},t \right) \\ \end{align}$$


Now we can define

$$\begin{align} & {{x}_{1}}=x \\ & {{x}_{2}}=\dot{x} \\ \end{align}$$


and then we can write previously rearranged second order differential equation in the following form. 

$$\begin{align} & {{{\dot{x}}}_{1}}={{x}_{2}} \\ & {{{\dot{x}}}_{2}}=f\left( {{x}_{1}},{{x}_{2}},t \right) \\ \end{align}$$



The following recurrence formula is used to find values of X(t) at different grid points t_i according to the fourth order Runge-Kutta method. 
$${{\vec{X}}_{i+1}}={{\vec{X}}_{i}}+\frac{1}{6}\left[ {{{\vec{K}}}_{1}}+2{{{\vec{K}}}_{2}}+2{{{\vec{K}}}_{3}}+{{{\vec{K}}}_{4}} \right]$$


Where:

$$\begin{align} & {{{\vec{K}}}_{1}}=h\vec{F}({{{\vec{X}}}_{i}},{{t}_{i}}) \\ & {{{\vec{K}}}_{2}}=h\vec{F}({{{\vec{X}}}_{i}}+\frac{1}{2}{{{\vec{K}}}_{1}},{{t}_{i}}+\frac{1}{2}h) \\ & {{{\vec{K}}}_{3}}=h\vec{F}({{{\vec{X}}}_{i}}+\frac{1}{2}{{{\vec{K}}}_{2}},{{t}_{i}}+\frac{1}{2}h) \\ & {{{\vec{K}}}_{4}}=h\vec{F}({{{\vec{X}}}_{i}}+{{{\vec{K}}}_{3}},{{t}_{i+3}}) \\ \end{align}$$

Central difference method for SDOF systems

Let’s say that the vibrations of the system are described by the following governing equation:

$$m\frac{{{d}^{2}}x}{d{{t}^{2}}}+c\frac{dx}{dt}+kx=F\left( t \right)$$

Let the duration over which the solution of the previous equation is required be divided into n equal parts of interval h = ∆t each. To obtain a satisfactory solution, we must select a time step ∆t that is smaller than critical time step (∆t_cri)^2. Let the initial conditions be given by x(t=0) = x_0 and dx(t=0)/dt = dx_0/dt.

Replacing the derivatives by the central differences and writing differential equation at grid point I gives:


The previous formula is called the recurrence formula and this formula permits to calculate the displacement of the mass (x_i+1) if we know the previous history of displacements at t_i and t_i-1 as well as the present external force F_i. Applying the previous formula continuously we’ll get the complete time history of the behavior of the system. Certain care has to be exercised in applying previous equation when i = 0 since both x_0 and x_-1 are needed in finding x_1, and the initial conditions provide only the values of x_0 and dx_0/dt, we need to find the value of x_-1.

Substituting the known values of x_0 and dx_0/dt into differential equation which describes the motion of the system we get the following:

$$\begin{align} & m\frac{{{d}^{2}}x}{d{{t}^{2}}}+c\frac{dx}{dt}+kx=F\left( t \right) \\ & x(t=0)={{x}_{0}} \\ & \dot{x}(t=0)={{{\dot{x}}}_{0}} \\ & m{{{\ddot{x}}}_{0}}+c{{{\dot{x}}}_{0}}+k{{x}_{0}}=F(t=0) \\ & m{{{\ddot{x}}}_{0}}=F(t=0)-c{{{\dot{x}}}_{0}}-k{{x}_{0}}/:m \\ & {{{\ddot{x}}}_{0}}=\frac{1}{m}\left( F(t=0)-c{{{\dot{x}}}_{0}}-k{{x}_{0}} \right) \\ \end{align}$$

Now it’s time to apply the central difference approximation using the formulas:

$$\begin{align} & {{{\dot{x}}}_{i}}=\frac{1}{2h}\left( {{x}_{i+1}}-{{x}_{i-1}} \right) \\ & {{{\ddot{x}}}_{i}}=\frac{1}{{{h}^{2}}}\left( {{x}_{i+1}}-2{{x}_{i}}+{{x}_{i-1}} \right) \\ \end{align}$$

For i = 0 the formula can be rewritten as:

$$\begin{align} & {{{\dot{x}}}_{0}}=\frac{1}{2h}\left( {{x}_{1}}-{{x}_{-1}} \right) \\ & {{{\ddot{x}}}_{0}}=\frac{1}{{{h}^{2}}}\left( {{x}_{1}}-2{{x}_{0}}+{{x}_{-1}} \right) \\ \end{align}$$

Inserting them into differential equation we obtain value of x_-1:

$${{x}_{-1}}={{x}_{0}}-\Delta t{{\dot{x}}_{0}}+\frac{{{\left( \Delta t \right)}^{2}}}{2}{{\ddot{x}}_{0}}$$

Examples of Two Degrees of Freedom Systems 1-5

Example 1 - Find the natural frequencies of the system shown in the following figure with: 


$${{m}_{1}}=m,{{m}_{2}}=2m,{{k}_{1}}=k,{{k}_{2}}=2k.$$

Determine the response of the system when k = 1000 N/m, m = 20 kg  and when initial values of displacements of the masses m_1 and m_2 are 1 and -1, respectively.

 Solution: First step is to derive the differential equation which describe the motion of the system. 

$$\begin{align} & {{m}_{1}}{{{\ddot{x}}}_{1}}+({{k}_{1}}+{{k}_{2}}){{x}_{1}}-{{k}_{2}}{{x}_{2}}=0, \\ & {{m}_{2}}{{{\ddot{x}}}_{2}}+{{k}_{2}}{{x}_{2}}-{{k}_{2}}{{x}_{1}}=0 \\ \end{align}$$

Solution of previous equation can be expressed in the following form:


$$\begin{align} & {{x}_{i}}(t)={{X}_{i}}\cos (\omega t+\phi ),\text{ i =1}\text{,2} \\ & {{x}_{1}}(t)={{X}_{1}}\cos (\omega t+\phi ) \\ & {{x}_{2}}(t)={{X}_{2}}\cos (\omega t+\phi ) \\ \end{align}$$

Inserting the solutions in previous differential equation we get:


$$\begin{align} & -{{m}_{1}}{{\omega }^{2}}{{X}_{1}}\cos (\omega t+\phi )+({{k}_{1}}+{{k}_{2}}){{X}_{1}}\cos (\omega t+\phi )-{{k}_{2}}{{X}_{2}}\cos (\omega t+\phi )={0}/{:}\;\cos (\omega t+\phi ), \\ & -{{m}_{2}}{{\omega }^{2}}{{X}_{2}}\cos (\omega t+\phi )+{{k}_{2}}{{X}_{2}}\cos (\omega t+\phi )-{{k}_{2}}{{X}_{1}}\cos (\omega t+\phi )={0}/{:}\;\cos (\omega t+\phi ) \\ & \left\{ -{{m}_{1}}{{\omega }^{2}}+({{k}_{1}}+{{k}_{2}}) \right\}{{X}_{1}}-{{k}_{2}}{{X}_{2}}=0, \\ & -{{k}_{2}}{{X}_{1}}+\left\{ -{{m}_{2}}{{\omega }^{2}}+{{k}_{2}} \right\}{{X}_{2}}=0 \\ & \det \left[ \begin{matrix} -{{m}_{1}}{{\omega }^{2}}+({{k}_{1}}+{{k}_{2}}) & -{{k}_{2}} \\ -{{k}_{2}} & -{{m}_{2}}{{\omega }^{2}}+{{k}_{2}} \\ \end{matrix} \right]=0 \\ & {{m}_{1}}{{m}_{2}}{{\omega }^{4}}-{{k}_{2}}{{m}_{1}}{{\omega }^{2}}-({{k}_{1}}+{{k}_{2}}){{m}_{2}}{{\omega }^{2}}+{{k}_{2}}({{k}_{1}}+{{k}_{2}})-k_{2}^{2}=0 \\ & {{m}_{1}}{{m}_{2}}{{\omega }^{4}}-{{k}_{2}}{{m}_{1}}{{\omega }^{2}}-({{k}_{1}}+{{k}_{2}}){{m}_{2}}{{\omega }^{2}}+{{k}_{2}}({{k}_{1}}+{{k}_{2}})-k_{2}^{2}={0}/{:{{m}_{1}}{{m}_{2}}}\; \\ & {{\omega }^{4}}-\left( \frac{{{k}_{1}}+{{k}_{2}}}{{{m}_{1}}}+\frac{{{k}_{2}}}{{{m}_{2}}} \right){{\omega }^{2}}+\frac{{{k}_{1}}{{k}_{2}}}{{{m}_{1}}{{m}_{2}}}=0 \\ \end{align}$$

Now it’s time to solve previously derived quadratic equation:


$$\begin{align} & {{\omega }^{4}}-\left( \frac{{{k}_{1}}+{{k}_{2}}}{{{m}_{1}}}+\frac{{{k}_{2}}}{{{m}_{2}}} \right){{\omega }^{2}}+\frac{{{k}_{1}}{{k}_{2}}}{{{m}_{1}}{{m}_{2}}}=0 \\ & \omega _{1,2}^{2}=\frac{\left( \frac{{{k}_{1}}+{{k}_{2}}}{{{m}_{1}}}+\frac{{{k}_{2}}}{{{m}_{2}}} \right)\pm \sqrt{{{\left( \frac{{{k}_{1}}+{{k}_{2}}}{{{m}_{1}}}+\frac{{{k}_{2}}}{{{m}_{2}}} \right)}^{2}}-4\frac{{{k}_{1}}{{k}_{2}}}{{{m}_{1}}{{m}_{2}}}}}{2}= \\ & \omega _{1,2}^{2}=\left( \frac{{{k}_{1}}+{{k}_{2}}}{2{{m}_{1}}}+\frac{{{k}_{2}}}{2{{m}_{2}}} \right)\pm \frac{1}{2}\sqrt{{{\left( \frac{{{k}_{1}}+{{k}_{2}}}{{{m}_{1}}}+\frac{{{k}_{2}}}{{{m}_{2}}} \right)}^{2}}-4\frac{{{k}_{1}}{{k}_{2}}}{{{m}_{1}}{{m}_{2}}}} \\ & \omega _{1,2}^{2}=\left( \frac{{{k}_{1}}+{{k}_{2}}}{2{{m}_{1}}}+\frac{{{k}_{2}}}{2{{m}_{2}}} \right)\pm \sqrt{\frac{1}{4}{{\left( \frac{{{k}_{1}}+{{k}_{2}}}{{{m}_{1}}}+\frac{{{k}_{2}}}{{{m}_{2}}} \right)}^{2}}-\frac{{{k}_{1}}{{k}_{2}}}{{{m}_{1}}{{m}_{2}}}} \\ \end{align}$$

The values of X­1 and X2 remain to be determined. These values depend on the natural frequencies ω1 and ω2. We shall denote these values of X1 and X2 corresponding to ω1 as X1(1) and X2(1) and those corresponding to ω2 as X1(2) and X2(2). All we need to do now is to define the ratios

$$\begin{align} & {{r}_{1}}=\frac{X_{2}^{(1)}}{X_{1}^{(1)}}=\frac{-{{m}_{1}}\omega _{1}^{2}+{{k}_{1}}+{{k}_{2}}}{{{k}_{2}}}=\frac{{{k}_{2}}}{-{{m}_{2}}\omega _{1}^{2}+{{k}_{2}}} \\ & {{r}_{2}}=\frac{X_{2}^{(2)}}{X_{1}^{(2)}}=\frac{-{{m}_{1}}\omega _{2}^{2}+{{k}_{1}}+{{k}_{2}}}{{{k}_{2}}}=\frac{{{k}_{2}}}{-{{m}_{2}}\omega _{2}^{2}+{{k}_{2}}} \\ \end{align}$$

General solution of differential equation can be written in the following form:


$$\begin{align} & {{x}_{1}}(t)=X_{1}^{(1)}\cos ({{\omega }_{1}}t+{{\phi }_{1}})+X_{1}^{(2)}\cos ({{\omega }_{2}}t+{{\phi }_{2}}) \\ & {{x}_{2}}(t)={{r}_{1}}X_{2}^{(1)}\cos ({{\omega }_{1}}t+{{\phi }_{1}})+{{r}_{2}}X_{1}^{(2)}\cos ({{\omega }_{2}}t+{{\phi }_{2}}) \\ \end{align}$$


The next step is to find variables inside the general solution according to these formulas:

$$\begin{align} & X_{1}^{(1)}=\frac{1}{{{r}_{2}}-{{r}_{1}}}{{\left[ {{\left\{ {{r}_{2}}{{x}_{1}}(0)-{{x}_{2}}(0) \right\}}^{2}}+\frac{{{\left\{ -{{r}_{2}}{{{\dot{x}}}_{1}}(0)-{{{\dot{x}}}_{2}}(0) \right\}}^{2}}}{\omega _{1}^{2}} \right]}^{1/2}} \\ & X_{1}^{(2)}=\frac{1}{{{r}_{2}}-{{r}_{1}}}{{\left[ {{\left\{ -{{r}_{1}}{{x}_{1}}(0)-{{x}_{2}}(0) \right\}}^{2}}+\frac{{{\left\{ {{r}_{1}}{{{\dot{x}}}_{1}}(0)-{{{\dot{x}}}_{2}}(0) \right\}}^{2}}}{\omega _{2}^{2}} \right]}^{1/2}} \\ & {{\phi }_{1}}={{\tan }^{-1}}\left\{ \frac{-{{r}_{2}}{{{\dot{x}}}_{1}}(0)+{{{\dot{x}}}_{2}}(0)}{{{\omega }_{1}}\left[ {{r}_{2}}{{x}_{1}}\left( 0 \right)-{{x}_{2}}(0) \right]} \right\} \\ & {{\phi }_{2}}={{\tan }^{-1}}\left\{ \frac{{{r}_{1}}{{{\dot{x}}}_{1}}(0)+{{{\dot{x}}}_{2}}(0)}{{{\omega }_{2}}\left[ {{r}_{1}}{{x}_{1}}\left( 0 \right)-{{x}_{2}}(0) \right]} \right\} \\ \end{align}$$

For m1 = m, m2 = m,  k1 = k and k2 = 2k we get: 


$$\begin{align} & \omega _{1,2}^{2}=\left( \frac{{{k}_{1}}+{{k}_{2}}}{2{{m}_{1}}}+\frac{{{k}_{2}}}{2{{m}_{2}}} \right)\pm \sqrt{\frac{1}{4}{{\left( \frac{{{k}_{1}}+{{k}_{2}}}{{{m}_{1}}}+\frac{{{k}_{2}}}{{{m}_{2}}} \right)}^{2}}-\frac{{{k}_{1}}{{k}_{2}}}{{{m}_{1}}{{m}_{2}}}} \\ & \omega _{1,2}^{2}=\left( \frac{k+2k}{2m}+\frac{2k}{4m} \right)\pm \sqrt{\frac{1}{4}{{\left( \frac{k+2k}{m}+\frac{2k}{2m} \right)}^{2}}-\frac{2{{k}^{2}}}{2{{m}^{2}}}} \\ & \omega _{1,2}^{2}=\left( \frac{6k+2k}{4m} \right)\pm \sqrt{\frac{1}{4}{{\left( \frac{6k+2k}{2m} \right)}^{2}}-\frac{2{{k}^{2}}}{2{{m}^{2}}}} \\ & \omega _{1,2}^{2}=\left( \frac{2k}{m} \right)\pm \sqrt{\frac{1}{4}{{\left( \frac{4k}{m} \right)}^{2}}-\frac{2{{k}^{2}}}{2{{m}^{2}}}}=\left( \frac{2k}{m} \right)\pm \sqrt{\frac{4{{k}^{2}}}{{{m}^{2}}}-\frac{2{{k}^{2}}}{2{{m}^{2}}}}=\left( \frac{2k}{m} \right)\pm \sqrt{\frac{8{{k}^{2}}-2{{k}^{2}}}{2{{m}^{2}}}} \\ & \omega _{1,2}^{2}=\left( \frac{2k}{m} \right)\pm \sqrt{\frac{8{{k}^{2}}-2{{k}^{2}}}{2{{m}^{2}}}}=\left( \frac{2k}{m} \right)\pm \sqrt{\frac{3{{k}^{2}}}{{{m}^{2}}}} \\ & \omega _{1}^{2}=(2-\sqrt{3})\frac{k}{m};\omega _{2}^{2}=(2+\sqrt{3})\frac{k}{m} \\ \end{align}$$

For k = 1000 N/m m = 20kg:


$$\begin{align} & \omega _{1}^{2}=(2-\sqrt{3})\frac{k}{m}=(2-\sqrt{3})\frac{1000}{20}=(2-\sqrt{3})50=13.397 \\ & {{\omega }_{1}}=3.6603\text{ rad/sec} \\ & \omega _{2}^{2}=(2+\sqrt{3})\frac{k}{m}=(2+\sqrt{3})50=186.602 \\ & {{\omega }_{2}}=13.6603\text{ rad/sec} \\ \end{align}$$

 $$\begin{align} & {{r}_{1}}=\frac{X_{2}^{(1)}}{X_{1}^{(1)}}=\frac{{{k}_{2}}}{-{{m}_{2}}\omega _{1}^{2}+{{k}_{2}}}=\frac{2k}{-2m\omega _{1}^{2}+2k}=\frac{k}{-m\omega _{1}^{2}+k}=1.36604 \\ & {{r}_{2}}=\frac{X_{2}^{(2)}}{X_{1}^{(2)}}=\frac{{{k}_{2}}}{-{{m}_{2}}\omega _{2}^{2}+{{k}_{2}}}=\frac{k}{-m\omega _{2}^{2}+k}=-0.36602 \\ \end{align}$$

With: 

$$\begin{align} & {{x}_{1}}\left( 0 \right)=1,{{{\dot{x}}}_{1}}\left( 0 \right)=0,{{x}_{2}}\left( 0 \right)=-1,{{{\dot{x}}}_{2}}\left( 0 \right)=0, \\ & X_{1}^{(1)}=\frac{1}{{{r}_{2}}-{{r}_{1}}}{{\left[ {{\left\{ {{r}_{2}}{{x}_{1}}(0)-{{x}_{2}}(0) \right\}}^{2}}+\frac{{{\left\{ -{{r}_{2}}{{{\dot{x}}}_{1}}(0)-{{{\dot{x}}}_{2}}(0) \right\}}^{2}}}{\omega _{1}^{2}} \right]}^{1/2}}=-0.36602 \\ & X_{1}^{(2)}=\frac{1}{{{r}_{2}}-{{r}_{1}}}{{\left[ {{\left\{ -{{r}_{1}}{{x}_{1}}(0)-{{x}_{2}}(0) \right\}}^{2}}+\frac{{{\left\{ {{r}_{1}}{{{\dot{x}}}_{1}}(0)-{{{\dot{x}}}_{2}}(0) \right\}}^{2}}}{\omega _{2}^{2}} \right]}^{1/2}}=-1.36603 \\ & {{\phi }_{1}}={{\tan }^{-1}}\left\{ \frac{-{{r}_{2}}{{{\dot{x}}}_{1}}(0)+{{{\dot{x}}}_{2}}(0)}{{{\omega }_{1}}\left[ {{r}_{2}}{{x}_{1}}\left( 0 \right)-{{x}_{2}}(0) \right]} \right\}=0 \\ & {{\phi }_{2}}={{\tan }^{-1}}\left\{ \frac{{{r}_{1}}{{{\dot{x}}}_{1}}(0)+{{{\dot{x}}}_{2}}(0)}{{{\omega }_{2}}\left[ {{r}_{1}}{{x}_{1}}\left( 0 \right)-{{x}_{2}}(0) \right]} \right\}=0 \\ \end{align}$$
$$\begin{align} & {{x}_{1}}(t)=-0.36602cos3.6603t-1.36603\cos 13.6603t \\ & {{x}_{2}}(t)=-0.5cos3.6603t+0.5\cos 13.6603t \\ \end{align}$$

Finite Difference Method - Exercises 1-3

Exercise 1.1: The forward difference formulas make use of the values of the function to the right of the base grid point. Thus the first derivative at point i(t = t_i) is defined as:

$$\frac{dx}{dt}=\frac{x(t+\Delta t)-x(t)}{\Delta t}=\frac{{{x}_{i+1}}-{{x}_{i}}}{\Delta t}$$
Derive the forward difference formulas for:

$$\frac{{{d}^{2}}x}{d{{t}^{2}}},\frac{{{d}^{3}}x}{d{{t}^{3}}}\text{ and }\frac{{{d}^{4}}x}{d{{t}^{4}}}\text{ at }{{t}_{i}}.$$

Solution:

$${{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i}}=\frac{{{\left. \frac{dx}{dt} \right|}_{i+1}}-{{\left. \frac{dx}{dt} \right|}_{i}}}{\Delta t}=\frac{\left( \frac{{{x}_{i+2}}-{{x}_{i+1}}}{\Delta t} \right)-\left( \frac{{{x}_{i+1}}-{{x}_{i}}}{\Delta t} \right)}{\Delta t}=\frac{{{x}_{i+2}}-2{{x}_{i+1}}+{{x}_{i}}}{{{\left( \Delta t \right)}^{2}}}$$ $${{\left. \frac{{{d}^{3}}x}{d{{t}^{3}}} \right|}_{i}}=\frac{{{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i+1}}-{{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i}}}{\Delta t}=\frac{\left( \frac{{{x}_{i+3}}-2{{x}_{i+2}}+{{x}_{i+1}}}{{{\left( \Delta t \right)}^{2}}} \right)-\left( \frac{{{x}_{i+2}}-2{{x}_{i+1}}+{{x}_{i}}}{{{\left( \Delta t \right)}^{2}}} \right)}{\Delta t}=\frac{{{x}_{i+3}}-3{{x}_{i+2}}+3{{x}_{i+1}}-{{x}_{i}}}{{{\left( \Delta t \right)}^{3}}}$$ $$\begin{align} & {{\left. \frac{{{d}^{4}}x}{d{{t}^{4}}} \right|}_{i}}=\frac{{{\left. \frac{{{d}^{3}}x}{d{{t}^{3}}} \right|}_{i+1}}-{{\left. \frac{{{d}^{3}}x}{d{{t}^{3}}} \right|}_{i}}}{\Delta t}=\frac{\left( \frac{{{x}_{i+4}}-3{{x}_{i+3}}+3{{x}_{i+2}}-{{x}_{i+1}}}{{{\left( \Delta t \right)}^{3}}} \right)-\left( \frac{{{x}_{i+3}}-3{{x}_{i+2}}+3{{x}_{i+1}}-{{x}_{i}}}{{{\left( \Delta t \right)}^{3}}} \right)}{\Delta t}= \\ & {{\left. \frac{{{d}^{4}}x}{d{{t}^{4}}} \right|}_{i}}=\frac{{{x}_{i+4}}-4{{x}_{i+3}}+6{{x}_{i+2}}-4{{x}_{i+1}}+{{x}_{i}}}{{{\left( \Delta t \right)}^{4}}} \\ \end{align}$$

Exercise 1.2: The backward difference formulas make use of the values of the function to the right of the base grid point. Thus the first derivative at point i(t = t_i) is defined as:

$$\frac{dx}{dt}=\frac{x(t)-x(t-\Delta t)}{\Delta t}=\frac{{{x}_{i}}-{{x}_{i-1}}}{\Delta t}$$

Derive the forward difference formulas for:

$$\frac{{{d}^{2}}x}{d{{t}^{2}}},\frac{{{d}^{3}}x}{d{{t}^{3}}}\text{ and }\frac{{{d}^{4}}x}{d{{t}^{4}}}\text{ at }{{t}_{i}}.$$+
Solution:

$${{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i}}=\frac{{{\left. \frac{dx}{dt} \right|}_{i}}-{{\left. \frac{dx}{dt} \right|}_{i-1}}}{\Delta t}=\frac{\left( \frac{{{x}_{i}}-{{x}_{i-1}}}{\Delta t} \right)-\left( \frac{{{x}_{i-1}}-{{x}_{i-2}}}{\Delta t} \right)}{\Delta t}=\frac{{{x}_{i}}-2{{x}_{i-1}}+{{x}_{i-2}}}{{{\left( \Delta t \right)}^{2}}}$$ $${{\left. \frac{{{d}^{3}}x}{d{{t}^{3}}} \right|}_{i}}=\frac{{{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i}}-{{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i-1}}}{\Delta t}=\frac{\left( \frac{{{x}_{i}}-2{{x}_{i-1}}+{{x}_{i-2}}}{{{\left( \Delta t \right)}^{2}}} \right)-\left( \frac{{{x}_{i-1}}-2{{x}_{i-2}}+{{x}_{i-3}}}{{{\left( \Delta t \right)}^{2}}} \right)}{\Delta t}=\frac{{{x}_{i}}-3{{x}_{i-1}}+3{{x}_{i-2}}-{{x}_{i-3}}}{{{\left( \Delta t \right)}^{3}}}$$ $$\begin{align} & {{\left. \frac{{{d}^{4}}x}{d{{t}^{4}}} \right|}_{i}}=\frac{{{\left. \frac{{{d}^{3}}x}{d{{t}^{3}}} \right|}_{i}}-{{\left. \frac{{{d}^{3}}x}{d{{t}^{3}}} \right|}_{i-1}}}{\Delta t}=\frac{\left( \frac{{{x}_{i}}-3{{x}_{i-1}}+3{{x}_{i-2}}-{{x}_{i-3}}}{{{\left( \Delta t \right)}^{3}}} \right)-\left( \frac{{{x}_{i-1}}-3{{x}_{i-2}}+3{{x}_{i-3}}-{{x}_{i-4}}}{{{\left( \Delta t \right)}^{3}}} \right)}{\Delta t}= \\ & {{\left. \frac{{{d}^{4}}x}{d{{t}^{4}}} \right|}_{i}}=\frac{{{x}_{i}}-4{{x}_{i-1}}+6{{x}_{i-2}}-4{{x}_{i-3}}+{{x}_{i-4}}}{{{\left( \Delta t \right)}^{4}}} \\ \end{align}$$

Exercise 1.3: Derive the formula for the fourth derivative according to the central difference method.
Solution:

$${{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i}}=\frac{{{x}_{i+1}}-2{{x}_{i}}+{{x}_{i-1}}}{{{\left( \Delta t \right)}^{2}}}$$ $$\begin{align} & {{\left. \frac{{{d}^{4}}x}{d{{t}^{4}}} \right|}_{i}}=\frac{{{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i+1}}-2{{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i}}+{{\left. \frac{{{d}^{2}}x}{d{{t}^{2}}} \right|}_{i-1}}}{{{\left( \Delta t \right)}^{2}}}=\frac{\left( \frac{{{x}_{i+2}}-2{{x}_{i+1}}+{{x}_{i}}}{{{\left( \Delta t \right)}^{2}}} \right)-2\left( \frac{{{x}_{i+1}}-2{{x}_{i}}+{{x}_{i-1}}}{{{\left( \Delta t \right)}^{2}}} \right)+\left( \frac{{{x}_{i}}-2{{x}_{i-1}}+{{x}_{i-2}}}{{{\left( \Delta t \right)}^{2}}} \right)}{{{\left( \Delta t \right)}^{2}}} \\ & {{\left. \frac{{{d}^{4}}x}{d{{t}^{4}}} \right|}_{i}}=\frac{{{x}_{i+2}}-4{{x}_{i+1}}+6{{x}_{i}}-4{{x}_{i-1}}+{{x}_{i-2}}}{{{\left( \Delta t \right)}^{4}}} \\ \end{align}$$

Finite Difference Method

The idea of the finite difference method or shortly FDM is to use approximations to derivatives. By applying FDM governing differential equation of motion is associated boundary conditions. There are three different types of formulas in  FDM, and these are:

  • Forward difference formulas,
  • Backward difference formulas and
  • Central difference formulas
These formulas are used to derive finite difference equations and the most accurate one are central difference formulas. In FDM we replace the solution domain (over which the solution of the given differential equation is required) with a finite difference number of points, referred to as mesh or grid points, and seek to determine the values of the desired solution at these points. The grid points are usually considered to be equally spaced along each of the independent coordinates.


First thing is to apply Taylor’s series expansion, x_i+1 and x_i-1 can be expressed about grid point i as: 

$${{x}_{i+1}}={{x}_{i}}+h{{\dot{x}}_{i}}+\frac{{{h}^{2}}}{2}{{\ddot{x}}_{i}}+\frac{{{h}^{3}}}{6}{{}_{i}}+...$$
$${{x}_{i-1}}={{x}_{i}}-h{{\dot{x}}_{i}}+\frac{{{h}^{2}}}{2}{{\ddot{x}}_{i}}-\frac{{{h}^{3}}}{6}{{}_{i}}+...$$
Where:  x_i is x at time t = t_i and h = t_i+1 – t_i which is equal to ∆t. By taking first two terms in previous equations and subtracting them we get: 

$$\begin{align} & {{x}_{i+1}}={{x}_{i}}+h{{{\dot{x}}}_{i}}\Rightarrow {{x}_{i}}={{x}_{i+1}}-h{{{\dot{x}}}_{i}} \\ & {{x}_{i-1}}={{x}_{i}}-h{{{\dot{x}}}_{i}}\Rightarrow {{x}_{i-1}}={{x}_{i+1}}-h{{{\dot{x}}}_{i}}-h{{{\dot{x}}}_{i}} \\ & {{x}_{i-1}}={{x}_{i+1}}-2h{{{\dot{x}}}_{i}} \\ & 2h{{{\dot{x}}}_{i}}={{x}_{i+1}}-{{x}_{i-1}} \\ & {{{\dot{x}}}_{i}}=\frac{1}{2h}\left( {{x}_{i+1}}-{{x}_{i-1}} \right) \\ \end{align}$$
By taking terms up to the second derivative and adding them together we get:

$$\left. \begin{align} & {{x}_{i+1}}={{x}_{i}}+h{{{\dot{x}}}_{i}}+\frac{{{h}^{2}}}{2}{{{\ddot{x}}}_{i}} \\ & {{x}_{i-1}}={{x}_{i}}-h{{{\dot{x}}}_{i}}+\frac{{{h}^{2}}}{2}{{{\ddot{x}}}_{i}} \\ \end{align} \right\}+$$ $$\begin{align} & {{x}_{i+1}}+{{x}_{i-1}}={{x}_{i}}+{{x}_{i}}+h{{{\dot{x}}}_{i}}-h{{{\dot{x}}}_{i}}+\frac{{{h}^{2}}}{2}{{{\ddot{x}}}_{i}}+\frac{{{h}^{2}}}{2}{{{\ddot{x}}}_{i}} \\ & {{x}_{i+1}}+{{x}_{i-1}}=2{{x}_{i}}+{{h}^{2}}{{{\ddot{x}}}_{i}} \\ & {{h}^{2}}{{{\ddot{x}}}_{i}}={{x}_{i+1}}-2{{x}_{i}}+{{x}_{i-1}} \\ & {{{\ddot{x}}}_{i}}=\frac{1}{{{h}^{2}}}\left( {{x}_{i+1}}-2{{x}_{i}}+{{x}_{i-1}} \right) \\ \end{align}$$

Numerical integration Methods in Vibration Analysis

In mechanics a two body system problem is very easy to solve while a problem of three or more bodies in the system is very difficult. That’s why we have to use approximate methods. From mathematical point of view differential equation of motion that describe a vibrating system cannot be integrated in the closed form that’s why we have to use a numerical approach.
There are several numerical methods available for solving vibration problems such as:
  • Finite difference Method,
  • Central difference Method
  • Runge-Kutta Method
  • Houbolt Method
  • Wilson Method
  • Newmark Method

Numerical integration methods have following fundamental characteristics:
  1. They are not intended to satisfy the governing differential equation(s) at all-time t but only at discrete time intervals Δt apart
  2. Suitable type of variation of displacement x, velocity dx/dt and acceleration d^2x/dt^2 is assumed within each time interval Δt.