User:EML4500.f08.JAMAMA/NM/ren

=HW 5.4=

Given: An ideal mass-spring-damper system
An ideal mass-spring-damper system with mass m, spring constant k and viscous damper of damping coefficient c is subject to a force F=u.

Find: Obtain the governing equation and plot system respond
1) Derive equation of motion in terms of $$d,c,k,m,u$$.

2) Let $$ x=\left\{ {\begin{array}{*{20}{c}} d \\  \dot d  \\ \end{array}} \right\} =\left\{ {\begin{array}{*{20}{c}} x_1  \\  x_2  \\ \end{array}} \right\}  \ $$, find $$\left(F,G\right)$$

3) Find $$c_{cr}$$ in terms of $$k,m$$ such that system is critically damped.

4) Let $$k=1,/ m=1/2,/ x_0=[0.8,-0.4]^T$$

a) For $$u=0$$, plot $$x_k$$ for $$c=\frac{1}{2}c_{cr}, c_{cr}, \frac{3}{2}c_{cr}$$.  b) For $$u=0.5$$ Gaussian noise, and $$c=\frac{3}{2}c_{cr}$$, plot $$x_k$$ c) For $$u=0.5$$ Cauchy noise, and $$c=\frac{3}{2}c_{cr}$$, plot $$x_k$$

Solution
1)

For the spring k 
 * {| style="width:100%" border="0"

$$ f_k=kd \ $$ $$ (Eq.5.4.1) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

For the damping c 
 * {| style="width:100%" border="0"

$$ f_c=c\cdot \frac{d}{dt}d(t)=c\dot d \ $$ $$ (Eq.5.4.2) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

For the mass m 
 * {| style="width:100%" border="0"

$$ m\frac{d^2}{dt^2}d(t) =u-f_k-f_c  \ $$ $$ (Eq.5.4.3) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Combine (Eq.5.4.1-3), we obtain the governing equation of the system


 * {| style="width:100%" border="0"

$$ m\ddot d = u - c\dot d-kd \quad \Rightarrow\quad m\ddot d + c\dot d+kd= u  \ $$ $$ (Eq.5.4.4) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

2)


 * {| style="width:100%" border="0"

$$ \underline{x}=\left\{ {\begin{array}{*{20}{c}} d \\  \dot d  \\ \end{array}} \right\} =\left\{ {\begin{array}{*{20}{c}} x_1  \\  x_2  \\ \end{array}} \right\}  \ $$ $$ (Eq.5.4.5) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

From (Eq.5.4.4), we have


 * {| style="width:100%" border="0"

$$ \left\{ {\begin{array}{*{20}{c}} \dot x_1 = x_2 \\  \dot x_2 = u - \frac{c}{m}x_2-\frac{k}{m}x_1\\ \end{array}} \right. $$   $$ (Eq.5.4.6) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

that is


 * {| style="width:100%" border="0"

$$ \left\{ {\begin{array}{*{20}{c}} \dot x_1 \\ \dot x_2 \\ \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} 0 & 1 \\  -\frac{k}{m} & -\frac{c}{m} \\ \end{array}} \right] \left\{ {\begin{array}{*{20}{c}} x_1 \\ x_2 \\ \end{array}} \right\} + \left[ {\begin{array}{*{20}{c}} 0 & 0 \\  1 & 0 \\ \end{array}} \right] \left\{ {\begin{array}{*{20}{c}} u \\ 0 \\ \end{array}} \right\} $$   $$ (Eq.5.4.7) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }


 * {| style="width:100%" border="0"

$$ \Rightarrow \quad \underline{A} = \left[ {\begin{array}{*{20}{c}} 0 & 1 \\  -\frac{k}{m} & -\frac{c}{m} \\ \end{array}} \right] \quad \underline{B}= \left[ {\begin{array}{*{20}{c}} 0 & 0 \\  1 & 0 \\ \end{array}} \right]$$ $$ (Eq.5.4.8) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

According to the lecture notes, we have

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \underline{F}=\underline{I}+\Delta t \underline{A} = \left[ {\begin{array}{*{20}{c}} 1 & \Delta t \\ -\frac{k\Delta t}{m} & 1-\frac{c\Delta t}{m} \\ \end{array}} \right] \quad \underline{G}=\Delta t \underline{B}= \left[ {\begin{array}{*{20}{c}} 0 & 0 \\  \Delta t & 0 \\ \end{array}} \right]$$ $$ (Eq.5.4.9) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

3) For critical damping system, the governing equation can be written as <span id="(1)">
 * {| style="width:100%" border="0"

$$ \ddot d + 2\omega_0\dot d+\omega^2 d=const. \ $$   $$ (Eq.5.4.10) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Comparing (Eq.5.4.10) with (Eq.5.4.4), we have

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{k}{m} = \left(\frac{c_{cr}}{2m}\right)^2 \quad\Rightarrow\quad c_{cr}=2\sqrt{mk}\ $$ $$ (Eq.5.4.11) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

4) The governing equation (Eq.5.4.7) can be integrated according to

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \underline{x}_{k+1}=\underline{F}\underline{x}_k+\underline{G}\underline{u}_{k+1} \quad\ $$ $$ (Eq.5.4.12) \ $$ with  $$\quad x_0=[0.8,-0.4]^T \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Matlab code is written to solve the problem and attached below:

Find: understand Kessler's code line by line
Using $$(p_{2i},\ p_{2i+1}),\ i=1,2,3$$ to understand Kessler's code line by line, starting with the best of Spring 2010

Solution
The Kessler's code is to compute the coefficients $$p_{2n}$$ and the constants $$c_{n}$$ where $$n=1,2,3,...,n$$ ,using the following algorithm:

1) $$i=0$$: <span id="(1)">
 * {| style="width:100%" border="0"

$$p_1(t)=c_1t,\ c_1=-1 $$ $$ (Eq.6.6.1) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

2) $$i=1,2,3,...$$: <span id="(1)">
 * {| style="width:100%" border="0"

$$c_{2i+1}=-\sum^{i-1}_{j=0}\frac{c_{2j+1}}{[2(i-j)+1]!}$$ $$ (Eq.6.6.2) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$p_{2i}(t)=\sum^{i}_{j=0}c_{2j+1}\frac{t^{2(i-j)}}{[2(i-j)]!}$$ $$ (Eq.6.6.3) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * function traperror

The input of this function $$n$$ is a integer that ranges from 1 to $$\infty$$. The function returns two arrays $$c$$ and $$p$$, corresponding to the constants $$c_{n}$$ and the coefficients $$p_{2n}$$.

Set initial values for f, g, cn and cd, where the initial value of cn and cd are the numerator and denominator of $$c_1=-1$$, respectively.

Start iteration from $$k=1$$ to $$k=n$$

Update vector f according to $$f^k_i=f^{k-1}_ig^{k-1}_i(g^{k-1}_i + 1)$$, where $$f^k_i$$ is the denominator of the factors of $$c_i$$ in (Eqn 6.6.2) as $$i=k$$

Calculate $$c_{2k-1}$$ according to (Eqn 6.6.2)

Update vector f and g, where $$g_i$$ and $$f_i$$ are respectively the numerator and denominator of the factors of $$c_{2i+1}$$ in (Eqn 6.6.2) as $$i=k$$

Calculate $$p_{2k}$$ according to (Eqn 6.6.3)

Add $$c_{2k+1}$$ and $$p_{2k}$$ to the output arrays c and p


 * function fracsum

This function is to calculate the summation of $$n_i/d_i$$ from vectors $$n$$ and $$d$$. The outputs nsum and dsum returns the numerator and denominator,respectively. $$\frac{nsum}{dsum}=\sum_{i}\frac{n_i}{d_i}$$

The function round is used to round off a number to its nearest integer. The function gcd returns an array containing the greatest common divisors of the corresponding elements of the two input arrays.

This is to remove the common divisors from $$n_i$$ and $$d_i$$, while the ratio $$n_i/d_i$$ is preserved.

The function lcm returns the least common multiple of corresponding elements of two input arrays

Given: inconsistent Trapezoidal Simpson algorithm

 * inconsistent Trapezoidal Simpson algorithm

<span id="(1)">
 * {| style="width:100%" border="0"

$$ z_{i+\frac{1}{2}}=\frac{1}{2}(z_{i}+z_{i+1}) \ $$ $$ (Eq.7.3.1) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * Hermite - Simpson algorithm

<span id="(1)">
 * {| style="width:100%" border="0"

$$ z_{i+\frac{1}{2}}=\frac{1}{2}(z_{i}+z_{i+1})+\frac{h}{8}(f_{i}-f_{i+1}) \ $$ $$ (Eq.7.3.2) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Find: Solve the logistic equation
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{dz}{dt}=rx(1-\frac{x}{x_{max}}) \ $$ $$ (Eq.7.3.3) \ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

with $$x_{max}=15$$, $$r=1.4$$, and two initial conditions

<span id="(1)">
 * {| style="width:100%" border="0"

$$x_0=3<\frac{1}{2}x_{max}$$
 * style="width:95%" |
 * style="width:95%" |

$$x_0=9<\frac{1}{2}x_{max}$$

$$t\in [0,20]$$
 * }
 * }

Solution
Then complete Hermite-Simpson time-stepping algorithm is given in. The only difference of the inconsistent Trapezoidal Simpson algorithm from the Hermite-Simpson algorithm is that the value of $$z_{i+1/2}$$ is evaluated as the mean value of $$z_{i}$$ and $$z_{i+1}$$, as shown in (Eq.7.3.2).



From observation of the above figures, the solutions given by the two algorithms almost coincide, regardless of the different initial conditions. Yet, the Hermite-Simpson algorithm is superior in the sense that it gives better estimation to the function value at the middle point.

128.227.113.35 19:25, 21 April 2011 (UTC)