User:Egm6341.s11.team4/HW7

=HW 7.0 problem: redo 5.2 and 5.4 cauchy noise=

Problem Statement
Refer to lecture notes [[media:nm1.s11.mtg25.djvu|25-2]] and [[media:nm1.s11.mtg25.djvu|26-1]] for detailed description 4). Run LSSM with Cauchy Random Noise



Problem Statement
Refer to lecture slide [[media:nm1.s11.mtg29.djvu|mtg-29]] for the problem statement. c).For u=0.5 Cauchy noise and c=$$1.5*{{C}_{cr}}$$, plot $$\underline$$ For u=0.5 Cauchy noise and c=$$1.5*{{C}_{cr}}$$, plot $$\underline$$



This problem was solved by shengfeng yang

=Problem 7.1: Integrating the Logistic Equation Using Hermite-Simpson, Forward Euler, and Backward Euler=

This problem is given on Lecture Slide 40-1 with references to lecture slides 39-1 and lecture slides 40-1 and lecture slides 40-2.

Given
The logistic equation

$$ \frac{dx}{dt}=rx(1-\frac{x}{x_{max}}) $$

With the following conditions

$$ x_{max}=15 $$; $$ r=1.4 $$; $$ t=[0,20] $$

And the two initial conditions

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

Find
1. The Integral over the time span using Hermite-Simpson and the step size $$h=2^{k} \hat{h} $$ for varying values of k.

2. The Integral over the time span using Forward Euler and the step size $$h=2^{k} \hat{h} $$ for varying values of k.

3. The Integral over the time span using Backward Euler and the step size $$h=2^{k} \hat{h} $$ for varying values of k.

Solution
All Code is included at the end of this problem. The Newton-Raphson method was used to solve the nonlinear problem in Hermite-Simpson. First, we will look at the case where the initial condition is $$ x_0=3 $$ then the case where $$ x_0=9 $$. Now, we we are to develop the forward and backward Euler integration schemes using the $$ x_{bar} $$ scheme introduced in class. It should be noted that we expect the explicit forward Euler scheme to produce unstable results and the implicit backward Euler scheme to be stable.

Forward Euler Development
Starting with the general forward Euler equation we will use the bar notation introduced in class.

$$ \bar{x}:=\frac{x}{x_{max}} $$

General Forward Euler Equation

$$ x_{i+1}=x_i (1+\frac{dx_i}{dt}h) $$

Which, for our equation, becomes

$$ \bar{x}_{i+1}=x_i(1+(r\bar{x}_i(1-\bar{x}_i))h) $$

Rearranging this, the equation becomes

$$ \bar{x}_{i+1}=hr\bar{x}_i(1-\bar{x}_i) +\bar{x}_i$$

Now our x vector can be obtained from this $$ \bar{x} $$ vector simply by multiplying through by the maximum value of x.

Backward Euler Development
Starting from the general Backward Euler equation.

$$ x_{i+1} = x_i + \frac{dx_{i+1}}{dt}h $$

$$ \bar{x}_{i+1}=\bar{x}_i + hr\bar{x}_{i+1}(1-\bar{x}_{i+1}) $$

Rearranging into a form that is second order in $$ \bar{x}_{i+1} $$.

$$ hr\bar{x}_{i+1}^{2}+(1-hr)\bar{x}_{i+1}-\bar{x}_i=0 $$

Now, employing the quadratic equation.

$$ \bar{x}_{i+1} = \frac{-(1-hr) \pm \sqrt{(1-hr)^2+4hr}} {2hr}$$

Once again, the true x vector can be found by multiplying the $$ \bar{x} $$ vector by the maximum value of x.

Both of these equations will be used when employing the Forward and Backward Euler integration.

Exact Solution
Starting with the logistic equation and separating the variables, it becomes


 * $$\displaystyle \frac{1}{x(1-\frac{x}{x_{max}})}\;dx=r\;dt$$

Now, multiply the left-hand side of the equation by $$\frac{x_{max}}{x_{max}}$$  to simplify the denominator


 * $$\displaystyle \frac{x_{max}}{x(x_{max}-x)}\;dx=r\;dt$$

The left-hand side can now be deconstructed into the partial fractions


 * $$\displaystyle \left[\frac{1}{x}+\frac{1}{x_{max}-x}\right]dx=r\;dt$$

Integrating both sides yields


 * $$\displaystyle ln|x|-ln|x_{max}-x|-ln|x_0|+ln|x_{max}-x_0|=r(t-t_0)$$

Taking the inverse natural log of both sides


 * $$\displaystyle \frac{x(x_{max}-x_0)}{x_0(x_{max}-x)}=e^{r(t-t_0)}$$

Now algebraically solving for x


 * $$\displaystyle x(x_{max}-x_0)=x_0(x_{max}-x)e^{r(t-t_0)}$$


 * $$\displaystyle x(x_{max}-x_0)=x_0x_{max}e^{r(t-t_0)}-x_0xe^{r(t-t_0)}$$


 * $$\displaystyle x(x_{max}-x_0+x_0e^{r(t-t_0)})=x_0x_{max}e^{r(t-t_0)}$$


 * $$\displaystyle x = \frac{x_0x_{max}e^{r(t-t_0)}}{x_{max}-x_0+x_0e^{(t-t_0)}}$$

Initial Condition X=3
Looking at figure 1 in the set of plots below, we see that the Hermite-Simpson integration method produces highly accurate results when plotted against the analytic solution. This accuracy maintains even as we increase the step size in accordance with the $$ h=2^k \hat{h} $$ rule. Similarly, as expected, the Forward Euler Integration scheme blows up as we increase the step size, ultimately providing results that are useless. The Backward Euler Integration scheme appears to be stable and provides results that rival the Hermite-Simpson integration scheme.

Initial Condition X=9
Looking at figure 2 in the set of plots below, we see similar results as obtained when the initial condition X=3 was used. First, we see that the Hermite-Simpson integration method once again produces highly accurate results when plotted against the analytic solution. This accuracy maintains even as we increase the step size in accordance with the $$ h=2^k \hat{h} $$ rule. Similarly the Forward Euler Integration scheme blows up as we increase the step size. Lastly, the Backward Euler Integration scheme is stable and provides results that rival the Hermite-Simpson integration scheme.

MATLAB Code
This problem was solved by Brendan Mahon

Given
The parameters for this problem are given in both the problem statement in the lectures as well as Subchan and Zbikowski's paper Computational optimal control of the terminal bunt manoeuvre—Part 1: Minimum altitude case.

From Subchan and Zbikowski's paper we have the following:

Equations of Motion:


 * {| style="width:70%" border="0" align="center"

(2.1)
 * $$\displaystyle \dot \gamma = \frac{T - D}{mV}sin(\alpha) + \frac{L}{mV}cos(\alpha) - \frac{g cos(\gamma)}{V}$$
 * 
 * 
 * }
 * }


 * {| style="width:70%" border="0" align="center"

(2.2)
 * $$\displaystyle \dot V = \frac{T - D}{m}cos(\alpha) - \frac{L}{m}sin(\alpha) - g sin(\gamma)$$
 * 
 * 
 * }
 * }


 * {| style="width:70%" border="0" align="center"

(2.3)
 * $$\displaystyle \dot x = V cos(\gamma)$$
 * 
 * 
 * }
 * }


 * {| style="width:70%" border="0" align="center"

(2.4)
 * $$\displaystyle \dot y = V sin(\gamma)$$
 * 
 * 
 * }
 * }

And,


 * {| style="width:70%" border="0" align="center"

\begin{bmatrix} \gamma \\ v \\ x \\ y \end{bmatrix} $$
 * $$\mathbf{z} =
 * $$\mathbf{z} =
 * }
 * }

Where &gamma;, V, x and y are state variables defined such that &gamma; is the flight path angle, V is the speed, x is the horizontal position and y is the altitude.

Also from Subchan and Zbikowski we have:


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle D(h,V,\alpha) = \frac{1}{2}C_d \rho V^2 S_{ref}, \,\,\, C_d = A_1 \alpha^2 + A_2\alpha + A_3$$
 * }
 * }
 * }


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle L(h,V,\alpha) = \frac{1}{2}C_l \rho V^2 S_{ref}, \,\,\, C_l = B_1\alpha + B_2$$
 * }
 * }
 * }


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle \rho = C_1y^2 + C_2y + C_3$$
 * }
 * }
 * }

Where,


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

\begin{align} & m = 1005 \; kg \\ & g = 9.81 \; 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.312 \times 10^{-9} \; kg/m^5 \\ & C_{2} = -1.142 \times 10^{-4} \; kg/m^4 \\ & C_{3} = 1.224 \; kg/m^3 \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

The initial conditions are given as:


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

\begin{align} & \gamma_0 = 0 \, rad\\ & V_0 = 272 \, m/s \\ & x_0 = 0 \, m \\ & y_0 = 30 \, m \\ \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

From the problem statement in the lecture notes we get the profiles for the two control inputs thrust, T(t) and angle of attack &alpha;(t):


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

T(t) = \begin{cases} 6000 \; N, & \mbox{for }0 \leq t < 27 \; \mbox{and} \; 33 \leq t \leq 40 \\ 1000 \; N, & \mbox{for }27 \leq t < 33 \end{cases} $$     (2.5)
 * $$\displaystyle
 * $$\displaystyle
 * 


 * }
 * }

For part 1:


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

\alpha(t) = \begin{cases} 0.03 \; rad, & \mbox{for }0 \leq t < 21 \\ \frac{0.13-0.09}{27-21}(t-21) + 0.09 \; rad, & \mbox{for }21 \leq t < 27 \\ \frac{-0.2+0.13}{33-27}(t-27) - 0.13 \; rad, & \mbox{for }27 \leq t < 33 \\ \frac{-0.13+0.2}{40-33}(t-33) - 0.2\; rad, & \mbox{for }33 \leq t \leq 40 \end{cases} $$     (2.6)
 * $$\displaystyle
 * $$\displaystyle
 * 
 * }
 * }

For part 2:


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

\alpha(t) = \begin{cases} 0.03 \; rad, & \mbox{for }0 \leq t < 21 \\ \frac{0.13-0.03}{27-21}(t-21) + 0.03 \; rad, & \mbox{for }21 \leq t < 27 \\ \frac{-0.2+0.13}{33-27}(t-27) - 0.13 \; rad, & \mbox{for }27 \leq t < 33 \\ - 0.2\; rad, & \mbox{for }33 \leq t \leq 40 \end{cases} $$     (2.7)
 * $$\displaystyle
 * $$\displaystyle
 * 
 * }
 * }

Objective
a.) Solve for z(t), for t contained on [0,40] such that successive numerical solutions are of the order 10-6, then plot the state variables as a function of time and compare them to the solution using the MatLab function "ode45".  Also compute


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle J_1 = \int_0^{t_f} y(t)dt \,\,\, s.t. \,\, y(t_f) = 0$$
 * }
 * }
 * }

Use given alpha for part 1.

b.) Solve for z(t), for t contained on [0,40] such that successive numerical solutions are of the order 10-6, then plot the state variables as a function of time and compare them to the solution using the MatLab function "ode45".


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle J_2 = \int_0^{t_f} y(t)dt \,\,\, s.t. \,\, y(t_f) = 0$$
 * }
 * }
 * }

Use given alpha for part 2.

Solution
Both parts of the problem will be solved using the same technique. The only difference will be the alpha function input into the Matlab code.

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


 * {| style="width:70%" border="0" align="center"

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

Where:


 * {| style="width:70%" border="0" align="center"

= \frac{1}{2}\left(\mathbf{z}_{i} + \mathbf{z}_{i+1}\right) + \frac{h}{8}\left(\mathbf{f}(\mathbf{z}_{i}) - \mathbf{f}(\mathbf{z}_{i+1})\right) $$     (2.9)
 * $$\displaystyle \mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}) = \mathbf{z}_{i+1/2}
 * $$\displaystyle \mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}) = \mathbf{z}_{i+1/2}
 * 
 * }
 * }

F(zi+1) can then be defined by subtracting the right side of eqn. (2.8) from the left side of eqn. (2.8).


 * {| style="width:70%" border="0" align="center"

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

From Lecture Slide [[media:nm1.s11.mtg39.djvu|39-1]] we get the following:


 * {| style="width:70%" border="0" align="center"

$$     (2.11)
 * $$\displaystyle \mathbf{z}_{i+1}^{(k+1)} = \mathbf{z}_{i+1}^{(k)} - \left[\frac{\partial \mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{\partial \mathbf{z}}\right]^{-1}\mathbf{F}(\mathbf{z}_{i+1}^{(k)})
 * $$\displaystyle \mathbf{z}_{i+1}^{(k+1)} = \mathbf{z}_{i+1}^{(k)} - \left[\frac{\partial \mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{\partial \mathbf{z}}\right]^{-1}\mathbf{F}(\mathbf{z}_{i+1}^{(k)})
 * 
 * }
 * }

In order to solve the problem we need to calculate the Jacobian as follows:


 * {| style="width:70%" border="0" align="center"

$$     (2.12) The values in the Jacobian can be seen in the MatLab code below.  "Main Script" 
 * $$\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{z}} = \left[\frac{\partial F_i}{\partial z_j}\right]
 * $$\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{z}} = \left[\frac{\partial F_i}{\partial z_j}\right]
 * 
 * }
 * }

 Function "zdot" 

 Function "dfdz" 

Function "nonlinode1" and "nonlinode2" to be used with the ode45 function.  Function "nonlinode1"   Function "nonlinode2" 

a.) Part 1: Initial values for alpha










b.) Part 2: Modified values for alpha










c.) Compute $$\ J_1$$, $$\ J_2$$ and compare.
For two different input values, calculate $$ J = \int_{t=0}^{t_f}y(t) dt, y(t_f) = 0 $$ The integration is calculated numerically by using "trapz" function of matlab, and the algorithm is in the above "main script" m-file. For input value of s10. $$ \ J_1 = 16803.19 $$ For input value of s11. $$ \ J_2 = 9024.61 $$ $$ \ J_1 > J_2 $$ So, using modified input makes better results.

Given
Refer to Lecture note S10 41-2 for more detailed problem statement

The logistic equation for population dynamics is given in S10 38-3,


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

$$\displaystyle \dot{x}(t) = f(x) = rx \left(1-\frac{x}{x_{max}}\right) $$ $$
 * <p style="text-align:right;">$$\displaystyle
 * }.
 * }.

Where, $$\displaystyle x_{max}=10;\ r = 1.2,\ $$ and $$\displaystyle t \in [0,10] $$

And the analytical solution is given by,


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

$$\displaystyle x(t) = \frac{x_0x_{max}e^{rt}}{x_{max} + x_0(e^{rt}-1)} $$ $$
 * <p style="text-align:right;">$$\displaystyle
 * }.
 * }.

Find
Integrate the population function using the inconsistent Trapezoidal-Simpson Algorithm given by,


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

$$\displaystyle z_{i+1/2} = \frac{1}{2} [z_{i} + z_{i+1}] $$
 * }.
 * }.

Two initial conditions, $$\displaystyle x_0=2\ and\ x_0 = 7,\ $$.

Solution
The Hermite-Simpson algorithm behaves better than inconsistent Trapezoidal-Simpson algorithm, that the error is much smaller comparing to inconsistent Trapezoidal-Simpson algorithm. Below are the plots for comparing Hermite-Simpson algorithm with Trapezoidal-Simpson algorithm.

Reference
S10.Team4

=Signatures=

Brendan Mahon 16:53, 17 April 2011 (UTC)

shengfeng yang 19:14 20 April 2011 (UTC)

William Kurth 03:18, 21 April 2011 (UTC) Taeho Kim 00:45, 21 April 2011 (UTC)

Erle Fields 05:05, 21 April 2011 (UTC)

Hailong Chen 10:58 21 April 2011 (EST)