User:Egm6341.s10.team3.sa/HW7

HW Problem Set 7 = Linearization about $$ x_{max} $$ = Ref: Lecture Notes [[media:Egm6341.s10.mtg38.djvu|p.38-4]]

Problem Statement

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

$$
 * $$\displaystyle x=\hat{x}+ y
 * $$\displaystyle x=\hat{x}+ y
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle f(x)=rx \left(1- \frac{x}{x_{max}} \right)
 * $$\displaystyle f(x)=rx \left(1- \frac{x}{x_{max}} \right)
 * }
 * }

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

$$ Linearization about equil. pt $$\hat {x} $$
 * $$\displaystyle \hat{x}=x_{max}
 * $$\displaystyle \hat{x}=x_{max}
 * }
 * }

Solution

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

$$
 * $$\displaystyle At \hat{x}=x_{max}
 * $$\displaystyle At \hat{x}=x_{max}
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

\begin{align} f(x)&=f(x_{max}+y)\\ &=r(x_{max}+y)(1- \frac{x_{max}+y}{x_{max}})\\ &=-ry \left(1+ \frac{y}{x_{max}} \right) \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle where yy<<x_{max} thus 1+ \frac{y}{x_{max}} \approx 1
 * $$\displaystyle where yy<<x_{max} thus 1+ \frac{y}{x_{max}} \approx 1
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle f(x)\approx -ry
 * $$\displaystyle f(x)\approx -ry
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle \frac{dx}{dt}= \frac{d(x_{max}+y)}{dt}= \frac{dy}{dt}
 * $$\displaystyle \frac{dx}{dt}= \frac{d(x_{max}+y)}{dt}= \frac{dy}{dt}
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle \Rightarrow \frac{dy}{dt}=-ry \Rightarrow y=y_0e^{-rt}
 * $$\displaystyle \Rightarrow \frac{dy}{dt}=-ry \Rightarrow y=y_0e^{-rt}
 * }
 * }

-- By Min Zhong 12:00, 22 Apr 2010 (UTC)

= Verify and solve the differential equation = Ref: Lecture Notes [[media:Egm6341.s10.mtg39.djvu|p.39-1]]

Problem Statement
The exact expression of x(t) is :
 * {| style="width:100%" border="0" align="left"

$$ i) Verify that:
 * $$\displaystyle x(t)= \frac{x_0x_{max}e^{rt}}{x_{max}+x_0(e^{rt}-1)}
 * $$\displaystyle x(t)= \frac{x_0x_{max}e^{rt}}{x_{max}+x_0(e^{rt}-1)}
 * }
 * }
 * {| style="width:100%" border="0" align="left"

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

ii) solve the above differential equation with hint on Lecture Notes [[media:Egm6341.s10.mtg40.djvu|p.40-1]]

Solution
i) Verify the differential equation
 * {| style="width:100%" border="0" align="left"

\begin{align} \frac{dx(t)}{dt}&= \frac {x_0x_{max}e^{rt}r}{x_{max}+x_0(e^{rt}-1)}+x_0x_{max}e^{rt}(-1) \left( \frac{1}{x_{max}+x_0(e^{rt}-1)} \right) ^2x_0e^{rt}r \\ &= \frac{x_ox_{max}e^{rt}}{x_{max}+x_0 \left( e^{rt}-1\right)}r+ \frac{r}{x_{max}}  \left( \frac{x_0x_{max}e^{rt}}{x_{max}+x_0 \left( e^{rt}-1\right)} \right)^2\\ &=xr-\frac{r}{x_{max}}x^2\\ &=rx \left(1- \frac{x}{x_{max}} \right) \end{align} $$ End the verification.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

ii)Solve the differential equation Partial fraction:
 * {| style="width:100%" border="0" align="left"

\begin{align} \frac{1}{x \left(1- \frac{x}{x_{max}} \right)}&= \frac{a}{x}+ \frac{b}{1- \frac{x}{x_{max}}}\\ &= \frac{a \left(1- \frac{x}{x_{max}}\right)+bx}{x \left(1- \frac{x}{x_{max}} \right)}\\ \end{align} $$ Let
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle a\left(1- \frac{x}{x_{max}}\right)+bx=1 \Rightarrow a-x \left(b- \frac{a}{x_{max}} \right)=1
 * $$\displaystyle a\left(1- \frac{x}{x_{max}}\right)+bx=1 \Rightarrow a-x \left(b- \frac{a}{x_{max}} \right)=1
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

\Rightarrow \begin{cases} & a= 1 \\ & b=  \frac{1}{x_{max}} \end{cases} $$ Separate the variables of the differential equation:
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle \int_{x_0}^{x}  \frac{dx}{x \left(1- \frac{x}{x_{max}} \right)}= \int_{0}^{t}rdt
 * $$\displaystyle \int_{x_0}^{x}  \frac{dx}{x \left(1- \frac{x}{x_{max}} \right)}= \int_{0}^{t}rdt
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle  \Rightarrow \int_{x_0}^{x}  \frac{1}{x} +  \frac{ \frac{1}{x_{max}}}{1- \frac{x}{x_{max}}}= \int_{0}^{t}rdt
 * $$\displaystyle  \Rightarrow \int_{x_0}^{x}  \frac{1}{x} +  \frac{ \frac{1}{x_{max}}}{1- \frac{x}{x_{max}}}= \int_{0}^{t}rdt
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

$$ :{| style="width:100%" border="0" align="left" $$ :{| style="width:100%" border="0" align="left" $$
 * $$\displaystyle  \Rightarrow  \left[ lnx-ln \left(1- \frac{x}{x_{max}} \right) \right]_{x_0}^x =  \left[rt \right]_{0}^t
 * $$\displaystyle  \Rightarrow  \left[ lnx-ln \left(1- \frac{x}{x_{max}} \right) \right]_{x_0}^x =  \left[rt \right]_{0}^t
 * }
 * }
 * $$\displaystyle  \Rightarrow  \left[ lnx-ln \left(1- \frac{x}{x_{max}} \right) \right]_{x_0}^x =  rt
 * $$\displaystyle  \Rightarrow  \left[ lnx-ln \left(1- \frac{x}{x_{max}} \right) \right]_{x_0}^x =  rt
 * }
 * }
 * $$\displaystyle  \Rightarrow ln \left( \frac{x}{x_0} \right) \left( \frac{1- \frac{x_0}{x_{max}}}{1- \frac{x}{x_{max}}} \right) =  rt
 * $$\displaystyle  \Rightarrow ln \left( \frac{x}{x_0} \right) \left( \frac{1- \frac{x_0}{x_{max}}}{1- \frac{x}{x_{max}}} \right) =  rt
 * }
 * }


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

$$
 * $$\displaystyle  \Rightarrow x= \frac{x_0x_{max}e^{rt}}{x_{max}+x_0(e^{rt}-1)}
 * $$\displaystyle  \Rightarrow x= \frac{x_0x_{max}e^{rt}}{x_{max}+x_0(e^{rt}-1)}
 * }
 * }

-- By Min Zhong 12:00, 22 Apr 2010 (UTC)

= Hermite-Simpson Algorithm to integrate Verhulst equation = Ref: Lecture notes [[media:Egm6341.s10.mtg39.djvu|p.39-1]]

Problem Statement

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


 * Given that,
 * }
 * }
 * }


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

\frac{dx}{dt} = rx \big(1-\frac{x}{xmax}\big) $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }


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

with $$\displaystyle xmax = 10 \;and\; r=1.2 $$ and the absolute tolerance being $$\displaystyle \theta(10^{-6})$$ for 2 cases of initial conditions Case(i) $$\displaystyle x_0 = 2$$ Case(ii) $$\displaystyle x_0 = 7$$
 * Use Hermite-Simpson rule to integrate the above given Verhulst equation
 * Use Hermite-Simpson rule to integrate the above given Verhulst equation
 * }
 * }

Solution

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

Hermite-Simpson algorithm.
 * A MATLAB code was written to solve the above differential equation using
 * A MATLAB code was written to solve the above differential equation using
 * }
 * }


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


 * MATLAB CODE:
 * MATLAB CODE:
 * }
 * }

Case(1) : Initial condition: $$ \displaystyle x_0 = 2.0 $$

Case(2) : Initial condition: $$ \displaystyle x_0 = 7.0 $$

Subramanian Annamalai 06:16, 19 April 2010 (UTC)

= Chaotic Systems (Chaos) =

Problem Statement
$$ \displaystyle Reproduce\, Fig.\, 15.6\, and\, 15.7\, in\, King\, et.\, al.\, 2003.\,\, (pp.\, 455-456)\, $$

Ref: Differential Equations: Amazon

$$ \displaystyle Logistic\,\, map\, $$

$$ \displaystyle x_{n+1}=r*x_n*(1-x_n) $$

$$ \displaystyle where\,\,\, r=constant,\, 0<r\leq 4\,\, and\,\, n=integer\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg40.djvu|p.40-1]]

Solution
$$ \displaystyle i)\, Reproduce\, the\, period\, doubling\, process\, as\, a\, bifurcation\, diagram\, using\, Matlab\, $$

$$ \displaystyle ii)\, Plot\, the\, bifurcation\, diagram\, $$



$$ \displaystyle iii)\, Reproduce\, a\, sequence\, of\, 100\, iterates\, of\, the\, logistic\, map\, with\, r=4,\, x_0=0.1\, and\, x_0=0.1+10^{-16}\, $$

$$ \displaystyle iiii)\, Plot\, the\, sequence\, $$



Heejun Chung 02:28, 22 April 2010 (UTC)

=Use the code of HS algorithm and run the code for $$ h = 2^k h $$ and also develop and run the code for Forward and backward Eular methods for the same step size =

Ref Lecture [[media:Egm6341.s10.mtg40.djvu|p.40-2]]

Problem Statement
Integrate the logistic equation and use the code of question 3 for HS and run the code for $$ h = 2^k h $$ and also develop Forward Eular and Back Eular methods for the same step sizes.

Solution
Matlab code with change in step size is given below.

MATLAB CODE:

For k = 1



For k = 2



For k=3



We can see that as we increase the value of k due to decrease in step size our curve shifted to the left side.

Below is the code for Forward Eular method

MATLAB CODE:

Forward Eular breaks down at higher values of step size h .This shows that it is unstable in nature. Below is the code for Backword Eular Method

MATLAB CODE:

For k = 1



For k=2



Abhishekksingh 10:29, 23 April 2010 (UTC)

= Solution of optimal control of bunt maneuver of aircraft using Hermite-Simpson (HS) algorithm and Newton-Raphson (NR) method for the given control inputs = Ref: Lecture notes [[media:Egm6341.s10.mtg40.djvu|p.40-3]] and [[media:Egm6341.s10.mtg41.djvu|p.41-1]]

Problem Statement


The given equations of motion of aircraft are


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

\mathbf{\dot{z}} = \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 - g\sin\gamma \\ v\cos\gamma \\ v\sin\gamma \end{Bmatrix} = \begin{Bmatrix} f_{1}(\mathbf{z}) \\ f_{2}(\mathbf{z}) \\ f_{3}(\mathbf{z}) \\ f_{4}(\mathbf{z}) \end{Bmatrix} = \mathbf{f}(\mathbf{z}) $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 1)
 * }
 * }

Here, $$\displaystyle D$$ and $$\displaystyle L$$ are
 * {| style="width:100%" border="0" align="left"

\begin{cases} D(y,v,\alpha) = \frac{1}{2}C_{d}\rho v^{2}S_{ref} \\ L(y,v,\alpha) = \frac{1}{2}C_{l}\rho v^{2}S_{ref} \end{cases} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

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

\begin{cases} C_{d}(\alpha) = A_{1}\alpha^2 + A_{2}\alpha + A_{3} \\ C_{l}(\alpha) = B_{1}\alpha + B_{2} \\ \rho(y) = C_{1}y^2 + C_{2}y + C_{3} \end{cases} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Solve the above equations of motion of aircraft when the below control inputs $$\displaystyle T(t)$$ and $$\displaystyle \alpha(t)$$ are given.


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

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} $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 2)
 * }
 * }


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

\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-21}(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} $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 3)
 * }
 * }

The solution should be obtained using Hermite-Simpson (HS) algorithm and Newton-Raphson (NR) method.

Solve the same equations of motion using "ode45" command in Matlab to verify the solution by HS and NR is correct.

Solution
All required parameters can be obtained from the Subchan and Zbikowski paper. (Subchan and Zbikowski, Optim. Control Appl. Meth. 2007; 28:311-353)
 * {| style="width:100%" border="0" align="left"

\begin{cases} 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{cases} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

The below Matlab function returns the given control input thrust $$\displaystyle T(t)$$.



The below Matlab function returns the given control input angle of attack $$\displaystyle \alpha(t)$$.



Let us use the notation $$\displaystyle \mathbf{z}(t_{i}) = \mathbf{z}_{i}$$ for simplification. In class meeting, we already obtained the below HS algorithm. (see [[media:Egm6341.s10.mtg37.djvu|p.37-2]])
 * {| style="width:100%" border="0" align="left"

\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
 * $$\displaystyle
 * $$\displaystyle (Eq. 4)
 * }
 * }

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

\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}\left(\mathbf{f}(\mathbf{z}_{i}) - \mathbf{f}(\mathbf{z}_{i+1})\right) $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 5)
 * }
 * }

The below Matlab function returns vector $$\displaystyle \mathbf{f}(\mathbf{z})$$ (see (Eq. 1)) for a given vector $$\displaystyle \mathbf{z}(t)$$.

By rearrange terms of (Eq.4), let us define $$\displaystyle \mathbf{F}(\mathbf{z}_{i+1})$$ as below.
 * {| style="width:100%" border="0" align="left"

\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] = 0 $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 6)
 * }
 * }

We can find the value of $$\displaystyle \mathbf{z}_{i+1}$$ which satisfies (Eq.6) using NR method. (see [[media:Egm6341.s10.mtg38.djvu|p.38-1]])
 * {| style="width:100%" border="0" align="left"

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

Since time step is small, $$\displaystyle \mathbf{z}_{i}$$ can be used as the initial value of $$\displaystyle \mathbf{z}_{i+1}^{(k)}$$, i.e. $$\displaystyle \mathbf{z}_{i+1}^{(0)}=\mathbf{z}_{i}$$.

Now, note that we need to calculate $$\displaystyle \frac{d\mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{d\mathbf{z}}$$ to use NR method.
 * {| style="width:100%" border="0" align="left"

\left[\frac{d\mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{d\mathbf{z}}\right] = \mathbf{I} - \frac{h/2}{3}\left[4\frac{d\mathbf{f}(\mathbf{g})}{d\mathbf{g}}\frac{\mathbf{g}(\mathbf{z}_{i+1})}{d\mathbf{z}_{i+1}} + \frac{d\mathbf{f}(\mathbf{z}_{i+1})}{d\mathbf{z}_{i+1}}\right] $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 8)
 * }
 * }

where $$\displaystyle \mathbf{I}$$ is a identity matrix. From (Eq. 5), we can calculate $$\displaystyle \frac{\mathbf{g}(\mathbf{z}_{i+1})}{d\mathbf{z}_{i+1}}$$ as below.
 * {| style="width:100%" border="0" align="left"

\frac{\mathbf{g}(\mathbf{z}_{i+1})}{d\mathbf{z}_{i+1}} = \frac{1}{2}\mathbf{I} - \frac{h}{8}\frac{d\mathbf{f}(\mathbf{z}_{i+1})}{d\mathbf{z}_{i+1}} $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 9)
 * }
 * }

Based on (Eq. 8) and (Eq. 9), we can calculate $$\displaystyle \frac{d\mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{d\mathbf{z}}$$ once we know $$\displaystyle \frac{d\mathbf{f}(\mathbf{z}_{i+1})}{d\mathbf{z}_{i+1}}$$. Therefore, let us calculate $$\displaystyle \frac{d\mathbf{f}(\mathbf{z})}{d\mathbf{z}}$$ using (Eq. 1).
 * {| style="width:100%" border="0" align="left"

\begin{align} \frac{d\mathbf{f}(\mathbf{z})}{d\mathbf{z}} & =
 * $$\displaystyle
 * $$\displaystyle

\begin{bmatrix} \frac{df_{1}(\mathbf{z})}{d\gamma} & \frac{df_{1}(\mathbf{z})}{dv} & \frac{df_{1}(\mathbf{z})}{dx} & \frac{df_{1}(\mathbf{z})}{dy}    \\ \frac{df_{2}(\mathbf{z})}{d\gamma} & \frac{df_{2}(\mathbf{z})}{dv} & \frac{df_{2}(\mathbf{z})}{dx} & \frac{df_{2}(\mathbf{z})}{dy}    \\ \frac{df_{3}(\mathbf{z})}{d\gamma} & \frac{df_{3}(\mathbf{z})}{dv} & \frac{df_{3}(\mathbf{z})}{dx} & \frac{df_{3}(\mathbf{z})}{dy}    \\ \frac{df_{4}(\mathbf{z})}{d\gamma} & \frac{df_{4}(\mathbf{z})}{dv} & \frac{df_{4}(\mathbf{z})}{dx} & \frac{df_{4}(\mathbf{z})}{dy} \end{bmatrix} \\ & = \begin{bmatrix} \frac{g}{v}\sin\gamma & -\frac{T}{mv^2}\sin\alpha - \frac{C_{d}\rho S_{ref}}{2m}\sin\alpha + \frac{C_{l}\rho S_{ref}}{2m}\cos\alpha + \frac{\rho}{v^2}\cos\gamma & 0 & -\frac{C_{d}(2C_{1}y+C_{2})vS_{ref}}{2m}\sin\alpha + \frac{C_{l}(2C_{1}y+C_{2})vS_{ref}}{2m}\cos\alpha    \\ -g\cos\gamma & -\frac{C_{d}\rho vS_{ref}}{m}\cos\alpha - \frac{C_{d}\rho vS_{ref}}{m}\sin\alpha & 0 & -\frac{C_{d}(2C_{1}y+C_{2})v^2S_{ref}}{2m}\cos\alpha - \frac{C_{l}(2C_{1}y+C_{2})v^2S_{ref}}{2m}\sin\alpha    \\ -v\sin\gamma & \cos\gamma & 0 & 0    \\ v\cos\gamma & \sin\gamma & 0 & 0 \end{bmatrix} \\ \end{align} $$ $$
 * $$\displaystyle (Eq. 10)
 * }
 * }

The below Matlab function returns vector $$\displaystyle \frac{d\mathbf{f}(\mathbf{z})}{d\mathbf{z}}$$ (see (Eq. 10)) for a given vector $$\displaystyle \mathbf{z}(t)$$.

Now, we can calculate (Eq. 6), (Eq. 8) and (Eq. 9) using the provided Matlab functions. Since we have all required information in order to perform NR method given in (Eq. 7), the solution of the equations of motion given in (Eq. 1) can be obtained.

The below Matlab function returns vector $$\displaystyle \mathbf{z}(t)$$ for the given time range $$\displaystyle t \in [0, 40]$$.

To verify the obtained solution from our own code, let us solve the same equations of motion using "ode45" command in Matlab.

Now, let us compare our results to the results from "ode45" command.

Plot $$\displaystyle \gamma(t) $$ vs $$\displaystyle t $$



Plot $$\displaystyle v(t) $$ vs $$\displaystyle t $$



Plot $$\displaystyle x(t) $$ vs $$\displaystyle t $$



Plot $$\displaystyle y(t) $$ vs $$\displaystyle t $$



Plot $$\displaystyle y(x) $$ vs $$\displaystyle x $$



As shown above figures, the solution from our own code is really consistent with the result from Matlab ode45 command.

Yong Nam Ahn 03:55, 22 April 2010 (UTC)

= ''' Inconsistent Trap. Simpson algorithm '''=

Problem Statement
$$ \displaystyle i)\, Solve\, logistic\, equation\, using\, the\, inconsistent\, trap.\, simpson\, algorithm\, $$

$$ \displaystyle z_{i+1} = z_i + \frac{(\frac{h}{2})}{3}*[f_i + 4f_{i+\frac{1}{2}}+f_{i+1}]\, $$

$$ \displaystyle with\,\, z_{i+1} = \frac{1}{2}*(z_i + z_{i+1})\, $$

$$ \displaystyle ii)\, Compare\, with\, the\, results\, of\, the\, Hermite-Simpson\, Algorithm\, using\, same\, values\, of\, h\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg41.djvu|p.41-2]]

Solution
$$ \displaystyle i)\, Matlab\, Code-\,\,\, only\, difference\, between\, the\, Inconsistent\, Trap.\, Simpson\, algorithm\, and\, the\, Hermite-Simpson\, Algorithm\, for\, Matlab\, is\, a\, definition\, of\, H_{half}\, $$

$$ \displaystyle a)\, Inconsistent\, Trap.\, Simpson\, algorithm\, $$

$$ \displaystyle x_{half} = 0.5*(x(i)+x(i+1));\, $$

$$ \displaystyle b)\, Hermite-Simpson\, Algorithm\, $$

$$ \displaystyle x_{half} = 0.5*(x(i)+x(i+1)) + (h/8)*(f1-f2);\, $$

$$ \displaystyle ii)\, Case(1):\, Initial\, condition:\, x_0 = 2.0\, $$



$$ \displaystyle iii)\, Case(2):\, Initial\, condition:\, x_0 = 7.0\, $$



$$ \displaystyle iiii)\, Comparison\, with the\, Hermite-Simpson\, Algorithm\, $$

$$ \displaystyle a)\, Case(1):\, Initial\, condition:\, x_0 = 2.0\, $$



$$ \displaystyle b)\, Case(2):\, Initial\, condition:\, x_0 = 7.0\, $$



$$ \displaystyle Therefore,\, the\, Inconsistent\, Trap.\, Simpson\, algorithm\, and\, the\, Hermite-Simpson\, Algorithm\, have\, the\, same\, results\, $$

Heejun Chung 16:00, 23 April 2010 (UTC)

= Linear Momentum Term in the x' Direction with tan(dr) =

Problem Statement
$$ \displaystyle Show\, \overline{CD} = \overline{AB} + hot\,\, for\, small\, dr\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg42.djvu|p.42-1]]

Solution
$$ \displaystyle i)\, Let's\, consider\, x'\, and\, y'\, coordinate\, system.\, $$



$$ \displaystyle ii)\, Define\, length\, of\, \overline{CD}\, and\, \overline{AB}\, with\, tan(dr)\, in\, x'\, and\, y'\, coordinate\, system.\, $$

$$ \displaystyle \overline{AB} = tan(dr) * V\, $$

$$ \displaystyle \overline{CD} = tan(dr) * (V+dV)\, $$

$$ \displaystyle iii)\, Identify\, tan(dr)\, Using\, Taylor\, Series\, expansion\, to\, approximate\, the\, tangent\, of\, the\, small\, angle.\, $$

$$ \displaystyle tan(dr) \cong dr\, $$

$$ \displaystyle iiii) Rewrite\, \overline{CD}\, and\, \overline{AB}\, with\, the\, small\, angle\, approximation\, of\, tangent.\, $$

$$ \displaystyle \overline{AB} = dr * V\, $$

$$ \displaystyle \overline{CD} = dr * (V+dV)\, = dr*V + \cancelto{}{dr*dV}\, $$

$$ \displaystyle The\, higher\, order\, term (dr*dV)\, can\, be\, negligible.\, $$

$$ \displaystyle Therefore\, \overline{CD} = \overline{AB}\, $$

$$ \displaystyle iiiii)\, Verification\, of\, the\, small\, angle\, approximation\, of\, tangent.\, $$

$$ \displaystyle From\, the\, above\, results,\, the\, small\, angle\, approximation\, of\, tangent\, may\, work\, pretty\, well\, out\, up\, to\, about\, 10\, degrees.\, $$

$$ \displaystyle Therefore\, we\, can\, conclude\, that\, \overline{CD} = \overline{AB}\, with\, small\, dr.\, $$

Heejun Chung 05:25, 23 April 2010 (UTC)

=  Parametrization of ellipse =

Problem Statement
$$ \displaystyle Show\, \frac{x^2}{a^2}+\frac{y^2}{b^2}=1\, $$

$$ \displaystyle where\,\,\, x=a*cos(t)\,\, and\,\, y=b*sin(t)\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg42.djvu|p.42-1]]

Solution
$$ \displaystyle i)\, For L.H.S,\,\, \frac{x^2}{a^2}+\frac{y^2}{b^2} = (\frac{x}{a})^2+(\frac{y}{b})^2\,\,\,\,\,\,\,\, (1)\, $$

$$ \displaystyle ii)\, x=a*cos(t)\, \rightarrow\, (\frac{x}{a})=cos(t)\,\, and\,\, y=b*cos(t)\, \rightarrow\, (\frac{y}{b})=sin(t)\,\,\,\,\,\,\,\, (2)\, $$

$$ \displaystyle (1) + (2) is given by $$

$$ \displaystyle (cos(t))^2 + (sin(t))^2\,\,\,\,\,\,\,\,(3) $$

$$ \displaystyle iii) (3) is always 1 because of Pythagorean Identities.\, $$

$$ \displaystyle Thus,\, \frac{x^2}{a^2}+\frac{y^2}{b^2} = 1\, $$

Heejun Chung 20:47, 22 April 2010 (UTC)

=  Arc Length =

Problem Statement
$$ \displaystyle Show\, dl=[dx^2+dy^2]^{\frac{1}{2}}\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg42.djvu|p.42-2]]

Solution
$$ \displaystyle i)\, Let's\, consider\, the\, Pythagorean\, theorem\, written\, by\, $$

$$ \displaystyle a^2 + b^2 = c^2\, $$

$$ \displaystyle where\, c\, represents\, the\, length\, of\, the\, hypotenuse,\, and\, a\, and\, b\, represent\, the\, lengths\, of\, the\, other\, two\, sides.\, $$

$$ \displaystyle ii)\, If\, the\, angle(dr)\, for\, dl\, is\, small,\, we\, can\, apply\, the\, Pythagorean\, theorem\, to\, the\, arc\, length\, formula\, like\, the\, below\, figure.\, $$



$$ \displaystyle iii) From\, the\, above\, figure,\, dl\, is\, given\, by\, $$

$$ \displaystyle dl^2 \cong dx^2 + dy^2\, $$

$$ \displaystyle Thus,\, dl = [dx^2 + dy^2]^{\frac{1}{2}}\, $$

Heejun Chung 22:31, 22 April 2010 (UTC)

= Circumference of Ellipse = Ref: Lecture notes [[media:Egm6341.s10.mtg42.djvu|p.42-2]]

Problem Statement

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

$$\displaystyle C = \int dl = a \int_{0}^{2\pi} [1-e^2cos^2t]^{1/2}\,dt$$ where, Eccentricity, $$ \displaystyle e = \Bigg[1-\frac{b^2}{a^2}\Bigg]^{1/2} $$
 * Show that the circumference of an ellipse is given by,
 * Show that the circumference of an ellipse is given by,
 * }
 * }

Solution

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

x = a cost and y = b sint
 * The parametric representation of an ellipse is given by,
 * The parametric representation of an ellipse is given by,
 * }
 * }


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

$$\displaystyle \frac {dy}{dt} = b cos t $$
 * $$\displaystyle \frac {dx}{dt} = -a sin t $$
 * $$\displaystyle \frac {dx}{dt} = -a sin t $$
 * }
 * }

From Problem 10 of this HW, we know,
 * {| style="width:100%" border="0" align="left"


 * $$ \displaystyle dl = [dx^2 + dy^2 ]^{1/2} $$
 * }
 * }
 * }


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

$$ \displaystyle dl = \Bigg[\bigg(\frac{dx}{dt}\bigg)^2 + \bigg(\frac{dy}{dt}\bigg)^2 \Bigg]^{1/2} dt $$
 * Multiplying and dividing by dt,
 * Multiplying and dividing by dt,
 * }
 * }


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

$$ \displaystyle = \Bigg[(-a sin t)^2 + (b cos t)^2\Bigg]^{1/2} dt$$
 * }
 * }


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

$$ \displaystyle = \Bigg[(a^2 sin^2 t) + (b^2 cos^2 t)\Bigg]^{1/2} dt$$
 * }
 * }


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

$$ \displaystyle = a\Bigg[(sin^2 t) + \bigg(\frac{b^2}{a^2}\bigg) cos^2 t)\Bigg]^{1/2} dt$$
 * }
 * }


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

$$ \displaystyle cos^2 t $$ inside the root
 * Now adding and subtracting ,
 * Now adding and subtracting ,
 * }
 * }


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


 * $$ \displaystyle dl = a\Bigg[(sin^2 t) + cos^2 t - cos^2 t +\bigg(\frac{b^2}{a^2}\bigg) cos^2 t)\Bigg]^{1/2} dt$$
 * }
 * }
 * }


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

$$ \displaystyle cos^2 t + sin^2 t =1 $$
 * Since,
 * Since,
 * }
 * }


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


 * $$ \displaystyle dl = a\Bigg[1- \Bigg(1-\bigg(\frac{b^2}{a^2}\bigg)\Bigg) cos^2 t\Bigg]^{1/2} dt$$
 * }
 * }
 * }


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

$$ \displaystyle e = \Bigg[1-\frac{b^2}{a^2}\Bigg]^{1/2} $$
 * Since,
 * Since,
 * }


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


 * $$ \displaystyle dl = a [1-e^2cos^2t]^{1/2} dt$$
 * }
 * }
 * }


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


 * And therefore,
 * }
 * }
 * }


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

$$\displaystyle C = \int dl = a \int_{0}^{2\pi} [1-e^2cos^2t]^{1/2}\,dt$$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = |
 * }
 * }

Subramanian Annamalai 06:16, 19 April 2010 (UTC)

=Convert circumference of ellipse formula from t to $$\alpha$$ variable =

Ref Lecture [[media:Egm6341.s10.mtg40.djvu|p.42-2]]

Problem Statement
Prove that $$ a \int\limits_{t=0}^{2\pi}[1-e^2 cos^2 t]^{1/2}dt = 4a \int\limits_{\alpha=0}^{\pi/2} [1-e^2 sin^2 \alpha]^{1/2}d\alpha   $$

Solution
$$ a \int\limits_{t=0}^{2\pi}[1-e^2 cos^2 t]^{1/2}dt $$

$$ = 4a \int\limits_{t=0}^{\pi/2}[1-e^2 cos^2 t]^{1/2}dt $$

Let $$ cos t = Sin \alpha $$

$$ - sin t dt = cos \alpha d\alpha $$

Using above values we obtain

$$ = 4a (-) \int\limits_{\alpha=0}^{\pi/2} [1-e^2 sin^2 \alpha]^{1/2} (-) \frac {cos \alpha}{sin t} d\alpha  $$

$$ = 4a (-) \int\limits_{\alpha=0}^{\pi/2} [1-e^2 sin^2 \alpha]^{1/2} (-) \frac {cos \alpha}{sin t} d\alpha   $$

$$ = 4a \int\limits_{\alpha=0}^{\pi/2} [1-e^2 sin^2 \alpha]^{1/2} \frac {cos \alpha}{(1- sin^2 \alpha)^(1/2)} d\alpha         $$

$$ = 4a \int\limits_{\alpha=0}^{\pi/2} [1-e^2 sin^2 \alpha]^{1/2} \frac {cos \alpha}{cos \alpha} d\alpha        $$

$$ = 4a \int\limits_{\alpha=0}^{\pi/2} [1-e^2 sin^2 \alpha]^{1/2}  d\alpha        $$

Hence Proved

Abhishekksingh 10:26, 23 April 2010 (UTC)

= Change of variable   = Ref: Lecture Notes [[media:Egm6341.s10.mtg42.djvu|p.42-2]]

Problem Statement
By Changing the variable of integration,
 * {| style="width:100%" border="0" align="left"

$$ show that:
 * $$\displaystyle  x=cos \theta
 * $$\displaystyle  x=cos \theta
 * }
 * }
 * {| style="width:100%" border="0" align="left"


 * $$\displaystyle  I= \int_{-1}^{+1}f(x)dx= \int_{0}^{\pi}f(cos  \theta)sin \theta d\theta
 * $$\displaystyle  I= \int_{-1}^{+1}f(x)dx= \int_{0}^{\pi}f(cos  \theta)sin \theta d\theta

$$
 * }
 * }

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

$$ and
 * $$\displaystyle  x=cos \theta \Rightarrow dx=-sin \theta d \theta
 * $$\displaystyle  x=cos \theta \Rightarrow dx=-sin \theta d \theta
 * }
 * }
 * {| style="width:100%" border="0" align="left"

\begin{cases} & x=-1;  \theta =\pi \\ & x=+1;\theta=0 \end{cases} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"


 * $$\displaystyle  I= \int_{-1}^{+1}f(x)dx= \int_{\pi}^{0}f(cos\theta)(-1)sin \theta d\theta=\int_{0}^{\pi}f(cos\theta)sin \theta d\theta
 * $$\displaystyle  I= \int_{-1}^{+1}f(x)dx= \int_{\pi}^{0}f(cos\theta)(-1)sin \theta d\theta=\int_{0}^{\pi}f(cos\theta)sin \theta d\theta

$$ End the proof -- By Min Zhong 12:00, 22 Apr 2010 (UTC)
 * }
 * }

= Find the item $$a_k$$ in cosine series of $$f(cos\theta)$$  = Ref: Lecture Notes [[media:Egm6341.s10.mtg42.djvu|p.42-3]]

Problem Statement
The known cosine series of f(x) is as follows:
 * {| style="width:100%" border="0" align="left"

$$ Find out the expression for $$a_k$$
 * $$\displaystyle  f(cos \theta)= \frac{a_0}{2}+ \sum_{k=1}^{\infty}a_kcos \left( k\theta\right)
 * $$\displaystyle  f(cos \theta)= \frac{a_0}{2}+ \sum_{k=1}^{\infty}a_kcos \left( k\theta\right)
 * }
 * }
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle a_k= \frac{2}{\pi} \int_{0}^{\pi}f(cos\theta)cos(k\theta)d\theta
 * $$\displaystyle a_k= \frac{2}{\pi} \int_{0}^{\pi}f(cos\theta)cos(k\theta)d\theta
 * }
 * }

Solution
The proof starts by multiplying $$cos(k\theta)$$ on both sides of cosine series of $$ f(cos\theta)$$
 * {| style="width:100%" border="0" align="left"

f(cos \theta) = \frac{a_0}{2}+ \sum_{k=1}^{\infty}a_kcos \left( k\theta\right) $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

f(cos \theta)cos(k\theta) = \frac{a_0}{2}cos(k\theta)+ \sum_{k=1}^{\infty}a_kcos \left( k\theta\right)^2 $$ Do integration from $$ [0,\pi] $$ on both sides of the above equation:
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }
 * {| style="width:100%" border="0" align="left"

\int_{0}^{\pi}f(cos \theta)cos(k\theta) = \int_{0}^{\pi}\frac{a_0}{2}cos(k\theta)+  \int_{0}^{\pi}\sum_{k=1}^{\infty}a_kcos \left( k\theta\right)^2 $$ where
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }
 * {| style="width:100%" border="0" align="left"

\begin{align} \int_{0}^{\pi}\frac{a_0}{2}cos(k\theta)&= \frac{a_0}{2k} \left[sin(k\theta) \right]_0^\pi\\ &= \frac{a_0}{2k} \left( sin(k\pi)-sin(0)\right)\\ &=0\\ \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

\begin{align} \int_{0}^{\pi}\sum_{k=1}^{\infty}a_kcos \left( k\theta\right)^2&= \sum_{k=1}^{\infty}a_k \int_{0}^{\pi}cos^2(k\theta)d\theta\\ &=\sum_{k=1}^{\infty}a_k\int_{0}^{\pi} \frac{cos(2\theta+1)}{2}\\ &= \frac{1}{2} \sum_{k=1}^{\infty}a_k \left[ \frac{1}{2}sin(2\theta)+\theta\right]_0^\pi\\ &= \frac{1}{2} \sum_{k=1}^{\infty}a_k \left(  \frac{1}{2} sin(2\pi)+\pi- \frac{1}{2}sin0-0\right)\\ &= \frac{\pi}{2} \sum_{k=1}^{\infty}a_k\\
 * $$\displaystyle
 * $$\displaystyle

\end{align} $$
 * }
 * }

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

\int_{0}^{\pi}f(cos \theta)cos(k\theta) = \frac{\pi}{2} \sum_{k=1}^{\infty}a_k $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }


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

\Rightarrow a_k= \frac{2}{\pi} \int_{0}^{\pi}f(cos\theta)cos(k\theta)d\theta $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

-- By Min Zhong 12:00, 23 Apr 2010 (UTC)

= Find the expression of integration by fourier transfer  = Ref: Lecture Notes [[media:Egm6341.s10.mtg42.djvu|p.42-3]]

Problem Statement
The following is known:
 * {| style="width:100%" border="0" align="left"

I= \int_{-1}^{+1}f(x)dx= \int_{0}^{\pi}f(cos\theta)sin\theta d\theta $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

The known cosine series of f(x) is as follows:
 * {| style="width:100%" border="0" align="left"

$$ Find the expression of Integration as:
 * $$\displaystyle  f(cos \theta)= \frac{a_0}{2}+ \sum_{k=1}^{\infty}a_kcos \left( k\theta\right)
 * $$\displaystyle  f(cos \theta)= \frac{a_0}{2}+ \sum_{k=1}^{\infty}a_kcos \left( k\theta\right)
 * }
 * }
 * {| style="width:100%" border="0" align="left"

$$
 * $$\displaystyle I=a_0+ \sum_{j=1}^{\infty} \frac{2a_{2j}}{1- \left( 2j\right)^2}
 * $$\displaystyle I=a_0+ \sum_{j=1}^{\infty} \frac{2a_{2j}}{1- \left( 2j\right)^2}
 * }
 * }

Solution

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

\begin{align} I&= \int_{0}^{\pi}f(cos\theta)sin\theta d\theta\\ &= \int_{0}^{\pi} \left[  \frac{a_0}{2}+ \sum_{j=1}^{\infty}a_{2j}cos(2j\theta)\right]sin\theta d\theta\\ &= \frac{a_0}{2}\int_{0}^{\pi}sin\theta d \theta +  \sum_{j=1}^{\infty} a_{2j}\int_{0}^{\pi}sin\theta cos(2j\theta) d \theta\\ &= - \frac{a_0}{2} \left[ cos\theta\right]_0^\pi+ \sum_{j=1}^{\infty}a_{2j} \int_{0}^{\pi} \frac{1}{2} \left [sin(1+2j) \theta +sin (1-2j)\theta \right] d\theta\\ &= a_0+ \sum_{j=1}^{\infty}a_{2j} \frac{1}{2} \left[\left( \frac{-1}{1+2j} cos(1+2j) \theta \right)_0^\pi+ \left(  \frac{-1}{1-2j}cos(1-2j) \theta\right) _0^\pi \right]\\ &=a_0+ \sum_{j=1}^{\infty}a_{2j} \frac{1}{2} \left[ \frac{2}{1+2j} + \frac{2}{1-2j}\right]\\ &=a_0+ \sum_{j=1}^{\infty}a_{2j}  \left( \frac{1}{2}\right) \frac{4}{1-(2j)^2}\\ &=a_0+ \sum_{j=1}^{\infty} \frac{2a_{2j}}{1-(2j)^2}\\ \end{align} $$ The end of proof
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

-- By Min Zhong 12:00, 23 Apr 2010 (UTC)

=Contributing Team Members=

1.Subramanian Annamalai 06:18, 19 April 2010 (UTC) Authored 3 and 11 Proof read 6 2. Yong Nam Ahn 03:55, 22 April 2010 (UTC) Authored 6 Proof read 4,7,8,9 and 10 3. Abhishekksingh 10:32, 23 April 2010 (UTC) Authored 5 and 12 Proof read 1,2,13,14 and 15 4. Heejun Chung 16:01, 23 April 2010 (UTC) Authored 4,7,8,9 and 10 Proof read 3 and 11 5. Min Zhong 16:10, 23 April 2010 (UTC) Authored 1,2,13,14,15 Proof read 5 and 12

=References= 1. Introduction to Numerical Methods,2nd Edition by Kendall E Atkinson Second Edition 2. Numerical Methods for Engineers,5th Edition by Steven C Chapra and Raymond P Canale