User:Egm6341.s11.team5/HW7

=Problem 7.0 - Redo HW 5.2 and 5.4= From [[media:Nm1.s11.mtg39.djvu|Mtg 39-2]]

Problem 5.2 - Linear State Space Model with Random Noise
From [[media:Nm1.s11.mtg25.djvu|Mtg 25-2]]

Given

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

$$ $$ Where
 * $$ \mathbf{x_{k+1}}= \mathbf{F} \mathbf{x_{k}}+\mathbf{G} \mathbf{w_{k+1}}
 * $$ \mathbf{x_{k+1}}= \mathbf{F} \mathbf{x_{k}}+\mathbf{G} \mathbf{w_{k+1}}
 * $$\displaystyle (5-1)
 * }
 * }
 * {| style="width:100%" border="0" align="left"

\mathbf{F}&= \left[ \mathbf{I} +\Delta \mathbf{A} \right]\\ &= \left( {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right) + .02\left( {\begin{array}{*{20}{c}} {-0.2}&1\\ { - 1}&{-0.2} \end{array}} \right)\end{align} $$ $$
 * $$\begin{align}
 * $$\begin{align}
 * $$\displaystyle (5-2)
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

1\\ 1\end{array}} \right) $$ $$
 * $$G =\alpha \left( {\begin{array}{*{20}{c}}
 * $$G =\alpha \left( {\begin{array}{*{20}{c}}
 * $$\displaystyle (5-3)
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

{x_{k}^1}\\ {x_{k}^2}\end{array}} \right), \mathbf{x_{0}} =  \left( {\begin{array}{*{20}{c}} {3}\\ {-2}\end{array}} \right) $$ $$
 * $$\mathbf{x_{k}} = \left( {\begin{array}{*{20}{c}}
 * $$\mathbf{x_{k}} = \left( {\begin{array}{*{20}{c}}
 * $$\displaystyle (5-4)
 * }
 * }

And$$ \mathbf{w_{k}}$$ is a series of random values according to some given distribution

Find
Find and Plot $$\mathbf{x_{k}}$$ for $$k=1,2,3...N$$ with $$\alpha=0.5,1,2$$. With no noise, Gaussian random noise, and Cauchy random noise. Highlighting the initial condition and equlibrium point.

Solution




Matlab code used to generate the data and plots.



The Cauchy random variable generating equation was found using the Mathematica command:

In: InverseCDF[CauchyDistribution[m,b],y]

Out: $$b \tan \left(\pi \left(y-\frac{1}{2}\right)\right)+m$$

where m is the mean, b is the half height, and y is the random variable.

The Equlibrium point was found by:

$$\lim_{k\to \infty } \, F^k.x = \left( \begin{array}{c} 0. \\ 0. \end{array} \right)$$

Author and proof-reader
[Author] J Davis 20:17, 20 April 2011 (UTC) [Proof-reader]

Problem 5.4 - Mass, Spring, Damper Analysis (Cauchy Analysis)
From [[media:Nm1.s11.mtg29.djvu|Mtg 29-1]] And User:Egm6341.s11.team5/HW5

Given
Diagram of the mass, spring, damper under consideration From [[media:Nm1.s11.mtg29.djvu|Mtg 29-1]]

Also the step size, h, equals 0.02

Find
For u=0.5 cauchy noise and $$c=\frac{3}{2}{c_{cr}}$$, plot $$\mathbf{x_{k}}$$

Solution
Using Cauchy noise as the forcing function for this system and again using the over damped value for c, the plot of x becomes



The heavy tails cause the system to react to much greater changes in input, which causes the large magnitude spikes seen in figure 3. As we can see the system never settles down to an equilibrium. The drastic spikes in the noise cause wild oscillations. Because the Cauchy distribution does not have a constant mean value, each time the matlab code was run a different graph was produced.

Author and proof-reader
[Author]

[Proof-reader]

=Problem 7.1 - Integrating the Verhulst Equation= From [[media:Nm1.s11.mtg40.djvu|Mtg 40-1]]

Given
The Verhulst equation as follows: $$\frac{dx}{dt}= rx\left( 1 - \frac{x}{x_{max}} \right)$$

Find
1. Use HS to integrate the Verhulst equation, st $$AbsTol=O(10^{-6})$$, using: $$x_{max}=15$$ and $$r=1.4$$. Case 1: $$x_{0}=3$$ Case 2: $$x_{0}=9$$

2. Develop Forward Euler method for step size $$h = 2^{k}h$$.

3. Develop Backwards Euler method for step size $$h = 2^{k}h$$.

Solution
1. HS Integration The logistics equation is integrated using the octave code at the end. The graphs are as follows: At $$x_{0}=3$$:  At $$x_{0}=9$$:

2. Forward Euler Method, with $$h=2^{k}\hat{h}$$. $$\frac{dx}{dt}=r\bar{x}(1-\bar{x})$$ $$\frac{\bar{x}_{i+1}-\bar{x}_{i}}{h}=r\bar{x}(1-\bar{x})$$ $$\bar{x}_{i+1}-\bar{x}_{i}=hr\bar{x}(1-\bar{x})$$ $$\bar{x}_{i+1}-\bar{x}_{i}=2^{k}\hat{h}r\bar{x}(1-\bar{x})$$ $$\bar{x}_{i+1}=2^{k}\hat{h}r\bar{x}(1-\bar{x})+\bar{x}_{i}$$ Applying the octave code, we get the following graph: Case 1: $$x_{0}=3$$, at $$k=1$$. $$x_{0}=3$$, at $$k=2$$. $$x_{0}=3$$, at $$k=3$$. Case 2: $$x_{0}=9$$, at $$k=1$$. $$x_{0}=9$$, at $$k=$$. $$x_{0}=9$$, at $$k=3$$. 3. Backwards Euler Method The Backwards Euler Method is found by solving for $$ \bar{x}_{i+1}$$ as follows: $$ \frac{d\bar{x}}{dt}=r\bar{x}_{i+1}(1-\bar{x}_{i+1})$$ $$ \frac{\bar{x}_{i+1}-\bar{x}_{i}}{h}=r\bar{x}_{i+1}(1-\bar{x}_{i+1})$$ $$ \bar{x}_{i+1}-\bar{x}_{i}=hr\bar{x}_{i+1}(1-\bar{x}_{i+1})$$ $$ \bar{x}_{i+1}-\bar{x}_{i}=hr\bar{x}_{i+1}-hr\bar{x}_{i+1}^{2}$$ $$ \bar{x}_{i+1}-\bar{x}_{i}-hr\bar{x}_{i+1}+hr\bar{x}_{i+1}^{2}=0$$ $$ \bar{x}_{i+1}-\bar{x}_{i}-hr\bar{x}_{i+1}+hr\bar{x}_{i+1}^{2}=0$$ $$hr\bar{x}_{i+1}^{2}+\bar{x}_{i+1}(1-h)-\bar{x}_{i}=0$$ $$\bar{x}_{i+1}=\frac{(hr-1)+ \sqrt{(hr-1)^{2}+4hr\bar{x}_{i}}}{2hr}$$ $$\bar{x}_{i+1}=\frac{(hr-1)- \sqrt{(hr-1)^{2}+4hr\bar{x}_{i}}}{2hr}$$ Applying the octave code, we get the following graph: Case 1: $$x_{0}=3$$, at $$k=1$$. $$x_{0}=3$$, at $$k=2$$. $$x_{0}=3$$, at $$k=3$$. Case 2: $$x_{0}=9$$, at $$k=1$$. $$x_{0}=9$$, at $$k=$$. $$x_{0}=9$$, at $$k=3$$.

Author and proof-reader
[Author]Cavalcanti

[Proof-reader]

=Problem 7.2 - Hermite-Simpson Algorithm = From [[media:Nm1.s11.mtg40.djvu|Mtg 40-1]]

Given
Use the Hermite-Simpson algorithm to determine the trajectory of an interceptor given the following inputs.


 * $$T(t)=\left\{\begin{matrix}

6000 N, & 0\leq t< 27, & 33\leq t\leq 40\\ 1000 N, & 27\leq t< 33 & \end{matrix}\right.$$

For part 1
 * $$\alpha(t)=\left\{\begin{matrix}

0.03 rad & 0\leq t< 21\\ \frac{0.13-0.09}{27-21}(t-21)+0.09 rad & 21\leq t< 27\\ \frac{-0.2+0.13}{33-27}(t-27)-0.13 rad & 27\leq t< 33\\ \frac{-0.13+0.2}{40-33}(t-33)-0.2 rad & 33\leq t\leq 40 \end{matrix}\right.$$

For part 2
 * $$\alpha(t)=\left\{\begin{matrix}

0.03 \; rad, & 0\leq t < 21 \\ \frac{0.13-0.03}{27-21}(t-21) + 0.03 \; rad, & 21\leq t < 27 \\ \frac{-0.2+0.13}{33-27}(t-27) - 0.13 \; rad, & 27\leq t < 33 \\ - 0.2\; rad, & 33 \leq t \leq 40 \end{matrix}\right.$$


 * $$\mathbf{f}=

\begin{bmatrix} \dot{\gamma}\\ \dot{v}\\ \dot{x}\\ \dot{y} \end{bmatrix}=\begin{bmatrix}

\frac{T-D}{mv}sin(\alpha)+\frac{L}{mv}cos(\alpha)-\frac{g}{v}cos(\gamma)\\ \frac{T-D}{m}cos(\alpha)+\frac{L}{m}sin(\alpha)-gsin(\gamma)\\ vcos(\gamma)\\ vsin(\gamma) \end{bmatrix}$$


 * $$\mathbf{z}=\begin{bmatrix}

\gamma\\ v\\ x\\ y \end{bmatrix}$$


 * $$DL=\left\{\begin{matrix}

D(\alpha,v,y)=0.5C_d\rho v^2S_{ref}\\ L(\alpha,v,y)=0.5C_c\rho v^2S_{ref} \end{matrix}\right. $$


 * $$\left\{\begin{matrix}

C_d(\alpha)=A_1\alpha^2+A_2\alpha+A_3\\ C_l(\alpha)=B_1\alpha+B_2\\ \rho(y)=C_1y^2+C_2y+C_3 \end{matrix}\right.$$


 * $$\left\{\begin{matrix}

m=1005 kg\\ g=9.81 \frac{m}{s^2}\\ S_{ref}=0.3376 m^2\\ A_1=-1.9431\\ A_2=-0.1499\\ A_3=0.2359\\ B_1=21.9\\ B_2=0\\ C_1=3.312e^{-9} \frac{kg}{m^5}\\ C_2=-1.142e^{-4} \frac{kg}{m^4}\\ C_3=1.224 \frac{kg}{m^3} \end{matrix}\right.$$


 * $$t\in [0,40]$$ seconds
 * $$h=1$$ seconds

Find
1.Reproduce results of spring 2010 using the same $$T(t)$$ and $$\alpha(t)$$.Compute $$J_1=\int_{0}^{t_f}y(t)dt$$ and show that $$y(t_f)=0$$. 2.Use the same $$T(t)$$ and $$\alpha(t)$$ to plot all quantities(state variables) versus time(t).Compute $$J_2=\int_{0}^{t_f}y(t)dt$$ and show that $$y(t_f)=0$$.Also compare $$J_2$$ and $$J_1$$

Solution
The solution will be very similar to to the ones used in last term.There is however a difference,being input of alpha function.

From equation(1) on Lecture-[[media:nm1.s11.mtg38.djvu|38-3]] we get:

$$\displaystyle \mathbf{z}_{i+1} = \mathbf{z}_{i} + \frac{h/2}{3}[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})]$$ where, $$\displaystyle \mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}) = \mathbf{z}_{i+1/2} = \frac{1}{2}\left(\mathbf{z}_{i} + \mathbf{z}_{i+1}\right) + \frac{h}{8}(\mathbf{f}(\mathbf{z}_{i}) - \mathbf{f}(\mathbf{z}_{i+1})$$

F(zi+1) can then be defined: $$
 * $$\displaystyle \mathbf{F}(\mathbf{z}_{i+1}) = \mathbf{z}_{i+1} - \mathbf{z}_{i} + \frac{h/2}{3}[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})] = \mathbf{0}

From Lecture [[media:nm1.s11.mtg39.djvu|39-1]] we get: $$\displaystyle \mathbf{z}_{i+1}^{(k+1)} = \mathbf{z}_{i+1}^{(k)} - [\frac{\partial \mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{\partial \mathbf{z}}]^{-1}\mathbf{F}(\mathbf{z}_{i+1}^{(k)}) $$

The Jacobian is as follows: $$\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{z}} = [\frac{\partial F_i}{\partial z_j}] $$

Matlab code for the Jacobian is as follows:

For initial values of alpha Computing $$\ J_1$$ $$ \ J_1 = 1.68*10^4 $$(for s10 values, obtained from code with alpha=1)

For Modified values of alpha Computing$$\ J_2$$. $$ \ J_2 = 9.0246*10^3 $$(for s11 values, obtained from code with alpha=2) $$ \ J_1 > J_2 $$ So, better results are obtained by using s11 values. The codes for sub functions:

Author and proof-reader
[Author] | Raghunathan

[Proof-reader]

=Problem 7.3 - Solve logistic equation using Inconsistent Trapezoidal-Simpson algorithm= From [[media:Nm1.s11.mtg41.djvu|Mtg 41-2]]

Given
The logistic equation for population dynamics is defined as follows.


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

$$ \dot{x}(t) = f(x) = rx \left(1-\frac{x}{x_{\text{max}}}\right) $$     (7-1)$$ ,where $$\displaystyle x_{\text{max}}=15;\ r = 1.4,\ $$ and $$\displaystyle t \in [0,20] $$ The analytical solution to the equation is given as follows.
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle
 * }


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

$$ x(t) = \frac{x_0x_{\text{max}}e^{rt}}{x_{\text{max}} + x_0(e^{rt}-1)} $$     (7-2)$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle
 * }

Find
Solve the logistic equation defined in Eq (7-1) for two initial conditions, $$\displaystyle x_0=3$$ and $$\displaystyle x_0=9 $$ using Inconsistent Trapezoidal-Simpson algorithm defined in Eq.(7-3).
 * {| style="width:100%" border="0"

$$ z_{i+1} = z_{i} + \frac{h/2}{3} \left[f_i + 4f_{(i+1/2)} +f_{i+1}\right] $$     (7-3)$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle
 * }
 * {| style="width:100%" border="0"

$$ z_{i+1/2} = \frac{1}{2} [z_{i} + z_{i+1}] $$     (7-4)$$
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle
 * }

Solution
Consider the Eq.(7-1) in which $$\displaystyle f_{i+1/2}$$ is a function of $$\displaystyle z_{i}$$ and $$\displaystyle z_{i+1}$$, where $$\displaystyle z_i, z_{i+1} \rightarrow$$ values of $$\displaystyle x$$ at $$\displaystyle t = t_i$$ and $$\displaystyle t = t_{i+1}$$.

Since, this is an initial value problem, $$\displaystyle z_i$$ is known. So, Eq. (7-3) as a whole becomes a function of $$\displaystyle z_{i+1}$$ given by,
 * {| style="width:100%" border="0"

$$ \displaystyle F(z_{i+1})=0 $$     (7-5)$$ The solution to Eq. (7-5) is found using Newton-Raphson methon.
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle
 * }
 * {| style="width:100%" border="0"

$$ z_{i+1}^{(k+1)} = z_{i+1}^{(k)} - \left[\frac{dF(z_{i+1}^{(k)})}{dz}\right]^{-1} F(z_{i+1}^{(k)}) $$     (7-6)$$ The iteration in Newton-Raphson algorithm starts with an initial guess as $$\displaystyle z_{i+1}^{0} = z_{i}$$, is stopped when the absolute tolerance condition is satisfied as expressed in Eq. (7-7).
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle
 * }


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

$$ \displaystyle \text{AbsTol} = ||z_{i+1}^{(k+1)} - z_{i+1}^{(k)}|| \le 10^{-6} $$     (7-7)$$  FIRST CASE : $$\displaystyle x_0 = 3$$ Matlab code:
 * style="width:95%" |
 * style="width:95%" |
 * $$\displaystyle
 * }
 * {| style="width:100%" border="0" align="left"

Plot:
 * style="width:50%; padding:10px; border:1px solid #888888" align="center" |
 * style="width:50%; padding:10px; border:1px solid #888888" align="center" |
 * }
 * }
 * }



SECOND CASE : $$\displaystyle x_0 = 9$$ Matlab code:
 * {| style="width:100%" border="0" align="left"

Plot:
 * style="width:50%; padding:10px; border:1px solid #888888" align="center" |
 * style="width:50%; padding:10px; border:1px solid #888888" align="center" |
 * }
 * }
 * }

Author and proof-reader
[Author] Shin [Proof-reader]

=References=

=Contributing members=