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.

Linear Equations and Matrix Inverse

Sympy has a Matrix class and associated function that allow the symbolic solution of systems of linear equations. We shall consider a simple example to show implementation of Matrix class.
>>> from sympy import symbols, Matrix
>>> x,y,z = symbols('x,y,z')
>>> A = Matrix(([3, 7], [4, -2]))
>>> print A
Matrix([[3, 7], [4, -2]])
>>> print A.inv()
Matrix([[1/17, 7/34], [2/17, -3/34]])
>>> b = Matrix((12*z, 5*z))
>>> print b
Matrix([[12*z], [5*z]])
>>> x = A.inv()*b
>>> print x
Matrix([[59*z/34], [33*z/34]])
>>> print x.subs({z:3.3}).evalf(4)
Matrix([[5.726], [3.203]])
>>> type(x)

Alternative method of solving the same problem is to construct the system as a matrix in augmented form.
>>> from sympy import Matrix, symbols, solve_linear_system
>>> from sympy import Matrix, symbols, solve_linear_system
>>> x,y,z = symbols('x,y,z')
>>> system = Matrix(([3, 7, 12*z], [4, -2, 5*z]))
>>> print system
Matrix([[3, 7, 12*z], [4, -2, 5*z]])
>>> sol = solve_linear_system(system, x, y)
>>> print sol
{x: 59*z/34, y: 33*z/34}
>>> type(sol)
>>> for k in sol.keys():
...     print k, '=', sol[k].subs({z:3.3}).evalf(4)
...
x = 5.726
y = 3.203
A third option is the solve() method, whose arguments include the individual symbolic equations, rather than any matrices.
>>> from sympy import symbols, solve, Eq
>>> x,y,z = symbols('x,y,z')
>>> solve((Eq(3*x+7*y, 12*z), Eq(4*x-2*y, 5*z)), x, y)
{x: 59*z/34, y: 33*z/34}
>>> solve((3*x+7*y-12*z, 4*x-2*y-5*z), x, y)

{x: 59*z/34, y: 33*z/34}

Series Expansion and Plotting

It is possible to expand many SymPy expressions as Taylor series. The series method makes this straightforward. At minimum, we must specify the expression and the variable in which to expand it. Optionally, we can also specify the point around which to expand, the maximum term number, and the direction of the expansion.
>>> from sympy import *
>>> x = Symbol('x')
>>> sin(x).series(x, 0)
x - x**3/6 + x**5/120 + O(x**6)
>>> series(sin(x), x, 0)
x - x**3/6 + x**5/120 + O(x**6)
>>> cos(x).series(x, 0.5, 10)
1.11729533119247 - 0.438791280945186*(x - 0.5)**2 + 0.0799042564340338*(x - 0.5)**3 + 0.0365659400787655*(x - 0.5)**4 - 0.00399521282170169*(x - 0.5)**5 - 0.00121886466929218*(x - 0.5)**6 + 9.51241148024212e-5*(x - 0.5)**7 + 2.17654405230747e- 5*(x - 0.5)**8 - 1.32116826114474e-6*(x - 0.5)**9 - 0.479425538604203*x + O((x - 0.5)**10, (x, 0.5))
In some cases, especially for numerical evaluation and plotting the results, it is necessary to remove the trailing O(n) therm. To do that just type .remove0(n)
>>> cos(x).series(x, 0.5, 10).removeO()
-0.479425538604203*x - 1.32116826114474e-6*(x - 0.5)**9 + 2.17654405230747e-5*(x - 0.5)**8 + 9.51241148024212e-5*(x - 0.5)**7 - 0.00121886466929218*(x - 0.5)**6 - 0.00399521282170169*(x - 0.5)**5 + 0.0365659400787655*(x - 0.5)**4 + 0.079904256 4340338*(x - 0.5)**3 - 0.438791280945186*(x - 0.5)**2 + 1.11729533119247
SymPy provide two inbuild function, Plot()  from the sympy.plotting module, and plot from sympy. mpmath.visualization. Generally for most purposes the Matplotlib should be the plotting tool of choice. (CHECK OUT THE POSTS MATPLOTLIB).
Here is an example of plotting Sympy Computation.
>>> from sympy import sin, series, Symbol
>>> import pylab
>>> x = Symbol('x')
>>> s10 = sin(x).series(x, 0, 10).removeO()
>>> s50 = sin(x).series(x, 0, 50).removeO()
>>> s = sin(x)
>>> xx = []
>>> y10 =[]
>>> y50 = []
>>> y = []
>>> for i in range(1000):
...     xx.append(i/100.0)
...     y10.append(float(s10.subs({x:i/100.0})))
...     y50.append(float(s50.subs({x:i/100.0})))
...     y.append(float(s.subs({x:i/100.0})))
...
>>> pylab.figure()

>>> pylab.plot(xx, y10, label='O(10)')
[]
>>> pylab.plot(xx, y50, label='O(50)')
[]
>>> pylab.plot(xx, y, label='sin(x)')
[]
>>> pylab.axis([0, 10, -4, 4])
[0, 10, -4, 4]
>>> pylab.xlabel('x')

>>> pylab.ylabel('f(x)')

>>> pylab.legned()
>>> pylab.savefig('sympy.pdf')
>>> pylab.savefig('sympy.png')
>>> pylab.show()


Ordinary Differential Equations

A good thing about Sympy that it has support for solving several kinds of ordinary differential equations via its dsolve command. We need to set up the ODE and pass it as the first argument, eq. The second argument is the function f(x) to solve for. An optimal third argument, hint, influences the method that dsolve uses: some methods are better-suited to certain classes of ODE’s, or will express the solution more simply, than others.
In order to set up ODE solver, we need a way to refer to the unknown function for which we a re solving as well as its derivatives. The Function and Derivative classes facilitate this:
>>> from sympy import Symbol, dsolve, Function, Derivative, Eq
>>> y = Function('y')
>>> x = Symbol('x')
>>> y_ = Derivative(y(x), x)
>>> print dsolve(y_ + 5*y(x), y(x))
y(x) == C1*exp(-5*x)
It is worth mentioning that dsolve has introduced a constant of integration C1. It will introduce as many constants as are required, and they will all be named Cn, whene n is an integer. Also the first argument to dsolve is taken to be equal to zero unless we use Eq() function to specify otherwise:
>>> print dsolve(y_ + 5*y(x), y(x))
y(x) == C1*exp(-5*x)
>>> print dsolve(Eq(y_ + 5*y(x), 0), y(x))
y(x) == C1*exp(-5*x)
>>> print dsolve(Eq(y_ + 5*y(x), 12), y(x))
y(x) == C1*exp(-5*x)/5 + 12/5
The results from dsolve are an instance of the Equality class. This has consequences when we wish to numerically evaluate the function and use the result elsewhere, because even after using subs() and evalf(), we still have an Equalitym not any sort of scalar. The way to evaluate the function to a number is via the rhs attribute of the Equality.
Note that here, we use z to store Equality returned by dsolve, even though it is an expression for a function called y(x), to emphasis the distinction between the Equality itself and the data it contains.
>>> z=dsolve(y_ + 5*y(x), y(x))
>>> print z
y(x) == C1*exp(-5*x)
>>> type(z)

>>> print z.rhs
C1*exp(-5*x)
>>> C1 = Symbol('C1')
>>> y3 = z.subs({C1:2, x:3})
>>> print y3
y(3) == 2*exp(-15)
>>> y3.evalf(10)
y(3) == 6.11804641e-7
>>> z.rhs.subs({C1:2, x:4}).evalf(10)
4.122307245e-9
>>> type(z.rhs.subs({C1:2, x:4}).evalf(10))

Sometimes the dsolve my return too general solution. One example is when there is possibility that some coefficients may be complex. If we know that, for example, they are always real and positive, we can provide dsolve this information to avoid the solution becoming unnecessarily complicated.
>> from sympy import *
>>> a, x = symbols('a, x')
>>> f = Function('f')
>>> print dsolve(Derivative(f(x), x, 2) - a**4*f(x), f(x))
f(x) == C1*exp(-x*sqrt(a**4)) + C2*exp(x*sqrt(a**4))
>>> print dsolve(Derivative(f(x), x, 4) - a**4*f(x), f(x))
f(x) == C1*exp(-x*(a**4)**(1/4)) + C2*exp(x*(a**4)**(1/4)) + C3*exp(-I*x*(a**4)**(1/4)) + C4*exp(I*x*(a**4)**(1/4))

Differentiation and Integration

Sympy is capable of carrying out differentiation and integration of many functions.
>>> from sympy import Symbol, exp, sin, sqrt, diff
>>> x = Symbol('x')
>>> y = Symbol('y')
>>> diff(sin(x),x)
cos(x)
>>> diff(sin(x),y)
0
>>> diff(10 + 2*x + 4*y + 10*x**2 + x**9, x)
9*x**8 + 20*x + 2
>>> diff(10 + 2*x + 4*y + 10*x**2 + x**9, y)
4
>>> diff(10 + 2*x + 4*y + 10*x**2 + x**9, x).subs(x,1)
31
>>> diff(10 + 2*x + 4*y + 10*x**2 + x**9, x).subs(x,1.5)
262.660156250000
>>> diff(exp(x),x)
exp(x)
>>> diff(exp(-x**2/2),x)
-x*exp(-x**2/2)
The Sympy diff() function takes a minimum of two arguments: the function to be differentiated and the variable with respect ot which the differentiation is performed. Higher derivatives may be calculated by specifying additional variables, or by adding and optional integer argument.
>>> diff(3*x**4,x)
12*x**3
>>> diff(3*x**4, x, x, x)
72*x
>>> diff(3*x**4, x, 3)
72*x
>>> diff(3*x**4*y**7, x , 2, y, 2)
1512*x**2*y**5
>>> diff(diff(3*x**4*y**7, x , x), y, y)
1512*x**2*y**5
Integration uses a similar syntax. For the indefinite case, specify the function and a variable with respect to which the integration is performed:
>>> from sympy import integrate
>>> integrate(x**2,x)
x**3/3
>>> integrate(x**2,y)
x**2*y
>>> integrate(sin(x),y)
y*sin(x)
>>> integrate(sin(x),x)
-cos(x)
>>> integrate(-x*exp(-x**2/2),x)
exp(-x**2/2)
We can calculate definite integrals by providing integrate() method with a tuple containing the variable of interest, the lower and the upper bounds. If several variables are specified, multiple integration is performed. When Sympy returns a result in the Rational class, it is possible to evaluate it to a floating-point representation at any desired precision.
>>> integrate(x*2, (x, 0, 1))
1
>>> integrate(x**2, x)
x**3/3
>>> integrate(x**2, x, x)
x**4/12
>>> integrate(x**2, x, x, y)
x**4*y/12
>>> integrate(x**2, (x, 0, 2))
8/3
>>> integrate(x**2, (x, 0, 2), (x, 0, 2), (y, 0, 1))
16/3
>>> float(integrate(x**2, (x, 0, 2)))
2.6666666666666665
>>> type(integrate(x**2, (x, 0, 2)))

>>> res_rational = integrate(x**2, (x, 0, 2))
>>> res_rational.evalf()
2.66666666666667
>>> res_rational.evalf(100)
2.6666666666666666666666666666666666666666666666666666666666666666666666666666666
66666666666666666667

Numeric Types

Sympy library has the numeric types Rational and Real. The Rational class represents a rational number as a pair of two integers: the numerator and the denominator, so Rational(1,4) represents 1/4, Rational(5,2) represents 5/2:
>>> a = sympy.Rational(1,10)
>>> a
1/10
>>> b = sympy.Rational(45,46)
>>> b
45/46
>>> a*b
9/92
>>> a-b
-101/115
>>> a+b
124/115
Rational class works with rational expressions exactly. This is in contrast to Python’s standard float data type which uses floating point representation to approximate numbers. We can convert the sympy.Rational type into a Python floating point variable using float or the evalf method of the Rational object. The evalf method can take an argument that specifies how many digits should be computed for the floating point approximation.
>>> c = sympy.Rational(5,6)
>>> c
5/6
>>> float(c)
0.8333333333333334
>>> c.evalf()
0.833333333333333
>>> c.evalf(20)
0.83333333333333333333

Symbol Function

In order to carry out any symbolic computation we have to define symbolic variables which means we need to invoke SYMPY library.
>>> from sympy import Symbol
>>> x = Symbol('x')
>>> type(x)

>>> y = Symbol('y')
>>> 2*x -x
x
>>> x + y + x + 10*y
2*x + 11*y
>>> y + x - y + 10
x + 10
Instead of creating variable by variable we can create multiple variables using Symbols function.

>>> import sympy
>>> x,y,z = sympy.symbols('x,y,z')
>>> x + 2 * y + 3 * z - x
2*y + 3*z
Once we’ve computed our term manipulation, it’s time to insert numbers. This can be done using ‘subs’ method.
>>> from sympy import symbols
>>> x,y,z=sympy.symbols('x,y,z')
>>> z = x + 2 * y
>>> print z
x + 2*y
>>> z.subs(x,10)
2*y + 10
>>> z.subs(x,10).subs(y,3)
16
>>> z.subs({x:10, y:3})
16
We can also substitute symbolic variable for another (we will use expression z = x + 2*y).
>>> z = x + 2*y
>>> z.subs(x,y)
3*y
>>> z.subs(x,y).subs(y,2)
6

Introduction to SymPy

In this section, we introduce some basic functionality of the SymPy (short from Symbolic Python) library. In contrast to numerical computation, in symbolic calculation we are processing and transforming generic variables.
The SymPy home page is http://sympy.org/, and provides the full (and up-to-date) documentation
for this library.
Symbolic calculation is very slow compared to oating point operation, and thus generally not for direct simulation. However, it is a powerful tool to support the preparation of code and symbolic work. Occasionally, we use symbolic operations in simulations to work out the most effcient numerical code, before that is executed.

If you’ve installed the Python (x,y) distribution you probably already have the Sympy but if not type in the following command in CMD if you’re working on the Windows operating system or Terminal if your using Linux operation system.
pip install sympy
If your using Anaconda distribution all you need to do is type in the CMD or Terminal is the following command.
conda install sympy
After installing is done open Python or Spyder which is GUI of Python and type in
import sympy
If the sympy is installed correctly then you shouldn’t get any error. 

Standard linear algebra operations

Matrix multiplication

Two arrays can be manipulated in the usual linear-algebra way using numpy.matrxmultiply. Here is an example:
>>> import numpy as N
>>> import numpy.random
>>> A = numpy.random.rand(5, 5)
>>> x = numpy.random.rand(5)
>>> b = N.dot(A, x)
>>> A
array([[ 0.78447957,  0.38316272,  0.31135453,  0.70738504,  0.07551859],
       [ 0.43140155,  0.55929425,  0.28448864,  0.38397626,  0.14849022],
       [ 0.88869344,  0.93048982,  0.70973804,  0.04341364,  0.27067106],
       [ 0.95290594,  0.17021306,  0.98377969,  0.73476449,  0.40728296],
       [ 0.05918968,  0.357716  ,  0.33474359,  0.99625301,  0.68345456]])
>>> x
array([ 0.47874927,  0.52254844,  0.53053969,  0.07469784,  0.67971556])
>>> b
array([ 0.84514734,  0.77933734,  1.47545378,  1.39880365,  0.93182836])

Solving systems of linear equations

To solve a system of equations Ax=b that is given in matrix form, we can use the linear algebra package of numpy:
>>> import numpy.linalg as LA
>>> x = LA.solve(A,b)
>>> x
array([ 0.47874927,  0.52254844,  0.53053969,  0.07469784,  0.67971556])

Computing Eigenvectors and Eigenvalues

Here is a small example that computes the Eigenvectors and Eigenvalues (eig) of the unity matrix.
>>> import numpy
>>> import numpy.linalg as LA
>>> A = numpy.eye(3)
>>> print A
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]
>>> evalues, evectors = LA.eig(A)
>>> print(evalues)
[ 1.  1.  1.]
>>> print(evectors)
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]

 Curve fitting of polynomials

Here is a small example that computes the [trivial] Eigenvectors and Eigenvalues (eig) of the unity matrix (eye)):
>>> import numpy
>>> xdata = numpy.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
>>> ydata = numpy.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
>>> z = numpy.polyfit(xdata, ydata, 3)
>>> print z
[ 0.08703704 -0.81349206  1.69312169 -0.03968254]
>>> p = numpy.poly1d(z)
>>> xs = [0.1 * i for i in range(50)]
>>> ys = [p(x) for x in xs]
>>> import pylab
>>> pylab.plot(xdata, ydata, 'o', label = 'data')
[]
>>> pylab.plot(xs, ys, label='fitted curve')
[]
>>> pylab.ylabel('y')

>>> pylab.ylabel('x')

>>> pylab.savefig('polyfit.pdf')
>>> pylab.show()


Convert from array to list or tuple

To create an array back to a list or tuple, we can use the standard python functions list(s) and tuple(s) which take sequence s as the input argument and return a list and tuple, respectively:
>>>import numpy as N
>>> a = N.array([1, 4, 10])
>>> a
array([ 1,  4, 10])
>>>list(a)
>>> tuple(a)
(1, 4, 10)

Creating Matrices using Numpy

Here are two ways of creating 2d-array:
1) By converting a list of lists (or tuples) into an array:
>>> x = N.array([[1,2,3], [4,5,6]])
>>> x
array([[1, 2, 3],
       [4, 5, 6]])
2) Using the zeros method to create a matrix with 5 rows and 5 columns
>>> x = N.zeros((5,5))
>>> print x
[[ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]]
The shape of a matrix can be queried like this.
>>> x.shape
(5, 5)
Individual elements can be accessed and set using this syntax:
>>> x = N.array([[1,2,3], [4,5,6]])
>>> x[0,0]
1
>>> x[0,1]
2
>>> x[0,2]
3
>>> x[:,0]
array([1, 4])
>>> x[0, :]
array([1, 2, 3])

Creating vectors using Numpy

The data structure that we’ll need the most often is a vector and here are few examples of how we can generate  a vector using numpy package in Python.
First we’ll create a list of numbers and then using numpy array function transform the list into a numpy array or in our case vector.
b = [0, 0.5, 1, 1.5]
print "b_type = " + str(type(b))
print "b = " + str(b)
x = np.array(b)
print "x_type = " + str(type(x))
print "x = " + str(x)
After running this code in Spyder you’ll get the following result.
b_type = 
b = [0, 0.5, 1, 1.5]
x_type = 
x = [ 0.   0.5  1.   1.5]
Now we’ll show you how to use numpy Array range function to create a vector. Type the following code.
>>> x = N.arange(0, 2, 0.5)
>>> print x
[ 0.   0.5  1.   1.5]
Function aranage() is short for array range. In our case we’ve wanted to build an array (vector) from 0 to 2 with increment 0.5. This means that our vector has 4 elements that is 0, 0.5, 1., 1.5. Well the problem is that 2 is left out and the reason for that is that upper limit of an array is not included. It’s the same procedure when you’re creating lists.
Now we’ll create a null-vector with 4 elements. So type the code bellow.
>>> x = N.zeros(4)
>>> print x
[ 0.  0.  0.  0.]
Once array is established, we can set and retrieve individual values. For example:
>>> x = N.zeros(4)
>>> x[0] = 3.4
>>> x[2] = 4
>>> print x
[ 3.4  0.   4.   0. ]
>>> print(x[0])
3.4
>>> print(x[0:-1])
[ 3.4  0.   4. ]
Note that once we have a vector we can perform calculations on every element in the vector with a single statement:
>>> x = N.arange(0, 2, 0.5)
>>> print x
[ 0.   0.5  1.   1.5]
>>> print x+10
[ 10.   10.5  11.   11.5]
>>> print x**2
[ 0.    0.25  1.    2.25]
>>> print(N.sin(x))
[ 0.          0.47942554  0.84147098  0.99749499]

Numerical Python - NUMPY

The NumPy package is short for Numerical Python. This package provides access to a new data structure called arrays which allow efficient vector and matrix operations and it also provides a number of linear algebra operations such as:
            Solving a system of linear equations,
            Computation of Eigenvectors and Eigenvalues and more.
Installing Numpy depends on the distribution of Python you have. If you installed the Python(x,y) distribution you’ll probably already have numpy installed. In order to check if you have Numpy installed all you need to do is type the following command in python console.
>>> import numpy as np 
If you don’t have Numpy installed you’ll get an ImportError exception which will tell you that there is no module named numpy.
In order to install the numpy package all you need to do is type in the; Cmd (cmd is short for Command Prompt) window if you’re working on  Windows operating system or Terminal if you’re working on Linux operating system; following command:
pip install numpy
After the installation procedure is done you can start using the numpy package. If you don’t have Pip install look here (hyperlink missing).
If you’ve installed the Anaconda distribution than you need to type the following command  (in Cmd or Terminal):
conda install numpy

And that’s it you’re good to go. 

Multiple Assignment with Dictionaries

Combining built in function items to transform dictionary into a list of tuples, tuple assignment and for loop you can create nice code pattern for traversing the keys and values of a dictionary in a single loop.
d = {'a':10, 'b':20, 'c': 30}
l = list()
for key,val in d.items():
    print key, val
    l.append((key,val))
print "l = " + str(l)
l.sort(reverse = True)
print "List of tuples after applying sort()function"
print "l = " + str(l)
Output of previous code is given below.
a 10
c 30
b 20
l = [('a', 10), ('c', 30), ('b', 20)]
List of tuples after applying sort()function
l = [('c', 30), ('b', 20), ('a', 10)]
The for loop has two iteration variable because items returns a list of tuples and key, val is a tuple assignment that successively iterates through each of the key/value pairs in the dictionary. For each iteration through the loop, both key and value are advanced to the next key/value pair in the dictionary.

Dictionaries and Tuples

Dictionaries have a method called items that returns a list of tuples, where each tuple is a key-value pair.
d = {'a':10, 'b':20, 'c': 30}
t = d.items()
print t
[('a', 10), ('c', 30), ('b', 20)]
As you can see from previous example, the items are in no particular order. So in this example we’ve started from defining a dictionary and then using built in function itmes() transformed the dictionary in into a tuple. Since the list of tuples is a list we can apply sort function on a list of tuples. Converting a dictionary to a list of tuples is a way to output the contents of a dictionary sorted by key.
d = {'a':10, 'b':20, 'c': 30}
t = d.items()
print "t = " + str(t)
t.sort()
print "list of tuples using sort() function"
print " t = " + str(t)


The output of previous example is given below.
t = [('a', 10), ('c', 30), ('b', 20)]
list of tuples using sort() function
 t = [('a', 10), ('b', 20), ('c', 30)]
The new list is sorted in ascending order by the key value. 

Tuple assignment

Tuples have unique syntactic features in Python programming language and that is the ability to have a tuple on the left hand side of an assignment statement. This allows you to assign more than one variable at a time when the left hand side is a sequence.
In this example we have a two element list (this is a sequence) and assign the first and second elements of the list to variables x and y in a single statement.
a = [‘Python’, ‘language’]
x,y = a 
print x 
print y 
Output of previous script is given below.
Python 
language 
Python translates the tuple assignment to be the following
a = [‘Python’, ‘language’]
x = a[0]
y = a[1]
print x 
print y 
Output is the same as in previous example. When we use a tuple on the left hand side of the assignment statement, we omit the parentheses, but the following is an equally valid syntax:
a = [‘Python’, ‘language’]
(x,y) = a
print x 
print y 
A particularly clever application of tuple assignment allows us to swap the values of two variables in a single statement.
x,y = y,x
In previous example both statements are tuples, but the left side is a tuple of variables while the right side is a tuple of expressions. Each value on the right side is assigned to its respective variable on the left side. All the expressions on the right side are evaluated before any of the assignments.
The number of variables on the left and the number of values on the right have to be the same otherwise you’ll get ValueError.
a, b = 1,2,3
Traceback (most recent call last):
  File "", line 1, in 
ValueError: too many values to unpack

Comparing Tuples

Comparison operator work with tuples and other sequences. For example let’s define two tuples.
>>> t1 = (0, 1, 2)
>>>t2 = (0, 4, 5)
>>> t1 < t2
True
>>>t1 > t2
False
When we type t1 < t2 the Python starts comparing the first element from each sequence. If the first element is equal in both sequences the Python is comparing next sequence elements until it finds those elements that differ. So in first statement t1 < t2 Python starts by comparing first element of t1 and first element of t2. Since they are both equal the Python is comparing the second element of t1 (1) and t2 (4). The 1 is less than 4 so that’s True. The final comparison is the third element of t1 (2) and t2 (5) and the answer is True which means that in deed the t1 is less than t2. The second statement t1 > t2 is False because elements of tuple t1 are smaller than elements in t2.
The sort function can be also applied on tuples. It sorts primarily by first element, but in the case of a tie, it sorts by second element. This feature lends itself to a pattern called DSU and DSU stands for:
DECORATE – sequence by building a list of tuples with one or more sort key preceding the elements from the sequence,
SORT – the list of tuples using Python built-in sort, and
UNDECORATE – by extracting the sorted elements of the sequence.
Try typing the following code.
txt = "Python is a programming language that lets you work quickly  and integrate systems more effectively." 
words = txt.split()
t = list()
for word in words: 
    t.append((len(word),word))
print 't = ' + str(t)
t.sort(reverse=True)
res = list()
for length, word in t: 
    res.append(word)

print res
Output is given below
t = [(6, 'Python'), (2, 'is'), (1, 'a'), (11, 'programming'), (8, 'language'), (4, 'that'), (4, 'lets'), (3, 'you'), (4, 'work'), (7, 'quickly'), (3, 'and'), (9, 'integrate'), (7, 'systems'), (4, 'more'), (12, 'effectively.')]
['effectively.', 'programming', 'integrate', 'language', 'systems', 'quickly', 'Python', 'work', 'that', 'more', 'lets', 'you', 'and', 'is', 'a']
The first loop builds a list of tuples, where each tuple is preceded by its length. Sort compare the first element, length, first and only considers the second element to break ties. The keyword argument reverse = True tells sort to go in decreasing order. 

Introduction to Tuples

Tuple is a sequence f values much like a list. The values stored in a tuple can be any type, and they are indexed by integers. The important difference is that tuples are immutable. Tuples are also comparable and hashable so we can sort lists of them and use tuples as key values in Python dictionaries.
Syntactically, a tuple is a comma separated list of values:
>>> t ='a','b', 'c', 'd','e'
>>> t
('a', 'b', 'c', 'd', 'e')
>>> type(t)
< type 'tuple'>
To create a tuple with a single element, you have to include comma after final element.
>>> t1 = ('a',)
>>> print t1
('a',)
>>> type(t1)

Without that final comma in parenthesis the Python would treat this declaration as a string.
>>> t1 = ('a')
>>> type(t1)

Another way to create a tuple is with using the built-in function tuple. With no argument, this function will create an empty tuple.
>>> t1 = tuple()
>>> type(t1)

If the argument is a sequence (string, list or tuple), the result of the call to tuple is tuple with the elements of a sequence:
>>> word = tuple('Python is awesome')
>>> print word
('P', 'y', 't', 'h', 'o', 'n', ' ', 'i', 's', ' ', 'a', 'w', 'e', 's', 'o', 'm','e')
Because tuple is a name of constructor, you should avoid using it as a variable name. Most list operators also work on tuples. The bracket operator indexes an element:
>>> z = ('a', 'b', 'c', 'd')
>>> print z[0]
a
>>> print z[1]
b
>>> i = 0
>>> while i < len(z):
...       print z[i]
...       i = i + 1
...
a
b
c
d
And the slice operator selects a range of elements:
>>> z[1:3]
('b', 'c')
If you try to modify one of the tuples elements the TypeError exception will occur.
>>> z[0] = 'a'
Traceback (most recent call last):
  File "", line 1, in 
TypeError: 'tuple' object does not support item assignment
You can modify the elements of a tuple, but you can replace one tuple with another.
>>> z = ('Z',) + z[1:]
>>> print z
('Z', 'b', 'c', 'd')

Looping and Dictionaries

If you use a dictionary as the sequence in a for statement, it traverses the key of the dictionary. The example is shown below.
guitars = {'gibson': 5050, 'fender': 500, 'jackson': 2000}
for key in guitars:
    print key, guitars[key]
After executing the code the result is.
jackson 2000
fender 500
gibson 5050
If we want to find all the entries in a dictionary above thousand, we could write the following code:
guitars = {'gibson': 5050, 'fender': 500, 'jackson': 2000}
for key in guitars:
    if guitars[key] > 1000:
        print key, guitars[key]
The result after executing previous block of code is:
jackson 2000
gibson 5050
We only see entries of the dictionary with value above 1000. If you want to print keys in alphabetical order, you first make a list of keys in the dictionary using keys method available in dictionary objects, and then sort that list and loop through the sorted list, looking up each key printing out key/value pairs in sorted order as follows:
guitars = {'gibson': 5050, 'fender': 500, 'jackson': 2000}
lst = guitars.keys()
print lst
lst.sort()
for key in lst:
    print key, guitars[key]
After executing the code we see that key elements of the dictionary are in alphabetical form.
['jackson', 'fender', 'gibson']
fender 500
gibson 5050
jackson 2000

Before entering the for loop we’ve printed out the list of keys in unsorted order that we’ve accomplished using keys method.