User:Egm6341.s11.team1.arm/HW5

Problem 1: Integrate using Composite Trapezoidal Rule and Composite Simpson's
 Referenced the open-source code gaussleg on the MathWorks File Exchange 

Find: Plots and error associated with the different discretizations
A.) For Eq. 1.1, plot $$ f_n^T(x) \ $$ (piecewise linear) and $$ f_n^S(x) \ $$ (piecewise quadratic) on $$ f(x) \ $$ B.) For Eq. 1.2, apply Trapezoidal and Simpson's Rule using uniform discretization C.) For Eq. 1.2, apply Trapezoidal and Simpson's Rule using non-uniform discretization D.) For Eq. 1.2, apply Trapezoidal and Simpson's Rule using Gauss-Legendre quadrature E.) For Eq. 1.3, apply Trapezoidal and Simpson's Rule using uniform discretization

Plot Trapezoidal Rule and Simpson's Rule
The following MATLAB (v2009b) code was developed to display area differences of integrand between Composite Trapezoidal rule and Composite Simpson's rule. Note that Langrange polynomials with n=2 were used to fit parabolas to $$f(x) \ $$.

Uniform discretization
The exact value of the integral in Eq. 1.2 is:

The following MATLAB code then developed to numerically integrate Eq 1.2 using uniform discretization and produce the following table.

Non-uniform discretization
In order to make h small close to zero, the following method of choosing x values was used:

The following MATLAB code then developed to numerically integrate Eq 1.2 using non-uniform discretization and produce the following table.

Gauss-Legendre Quadrature
The following MATLAB code was developed to evaluate the integral in Eq. 1.2.

NOTE: The user defined function legendrepoly was adapted from the work done in the open-source function gaussleg, available on the MathWorks File Exchange.

Periodic functions
Wolfram Alpha was used to calculate the exact value of the integral in Eq. 1.3.

The following MATLAB code was developed to numerically integrate Eq. 1.3.

Problem 2: Linear Space Model
 Solved without assistance 

Run linear space model and plot $$ f \underline{x} _j $$, for $$ j=1,2,... \ $$ in space $$ (x^1,x^2)=(x,y) \ $$
Given the following:

Solution

The following MATLAB code was developed to evaluate the state space model until an equilibrium point was reached. The results are shown in Figure 5.2.1.


 * Nm1.s11.team1.HW4.fig5.2.1.png

====Find equilibrium point as $$ \lim_{k \to \infty} \underline{x} _{k+1} = \lim_{k \to \infty} \underline{F} ^{k+1} \underline{x} _0 =: \underline{\hat{a} } $$ ==== Plot $$ \underline{\hat{a}} $$ as big red dot and $$ x_0 \ $$ as big blue dot in the same plot for $$ \left \lbrace \underline{x} _j \right \rbrace \ $$, where $$ j=1,2,... \ $$

Solution

In order to determine the equilibrium point, the Euclidean distance between the next $$x_{n+1} \ $$ and $$ x_{n} \ $$ was calculated for each step in the model. In this case, the Euclidean distance is defined as:

When this distance became less than $$ 10^{-6} \ $$, the equilibrium point was taken as "found".

The MATLAB code presented above displays the equilibrium point, which is also plotted in Figure 5.2.1 above.

Gaussian Random Noise
Let $$ \underline{G} = \begin{Bmatrix} 1 \\ 1 \end{Bmatrix} \alpha\ $$, thus $$ \underline{w} -{kH} = \begin{Bmatrix} w_{kH} \end{Bmatrix} $$. Use MATLAB randn to gererate $$ \begin{Bmatrix} w_j \end{Bmatrix} $$, for $$ j=0,1,2,... \ $$ and plot $$ \begin{Bmatrix} \underline{x} _j \end{Bmatrix} $$, for $$ j=0,1,2,... \ $$. Use a big blue dot for $$ j=0 \ $$, use small dots for the other points. All for $$ \alpha\ = 0.5 $$, $$ \alpha\ = 1 $$ and  $$ \alpha\ = 2 $$

Solution

The following MATLAB code was developed to evaluate the state space model for 1000 iterations, since an equilibrium point was not reached. The results are shown in Figure 5.2.2.


 * Nm1.s11.team1.HW4.fig5.2.2.png

Cauchy Random Noise
Same as 2.3 except with Cauchy. HINT: Find a MATLAB command to generate $$ \begin{Bmatrix} \theta\ _j \end{Bmatrix} $$, for $$ j=0,1,2,... \ $$ in a single slit diffraction.

Solution

A set of MATLAB functions freely available on the MathWorks File Exchange called cauchy was used to generate random values from the Cauchy distribution.

The following MATLAB code was developed to evaluate the state space model for 1000 iterations, since an equilibrium point was not reached. The results are shown in Figure 5.2.2.


 * Nm1.s11.team1.HW4.fig5.2.3.png

Problem 5: Modify your MATLAB code to make it more efficient
 Solved without assistance 

Refer to Homework 2.4 for original code. Improved code is shown below, along with the resulting table similar to the one generated from the original code. The original code took 0.016680 seconds to run; the modified code below took 0.004149 seconds to run, meaning the new code runs in approximately 1/4th of the time.

Do Romberg Table; Compare to previous results in HW2.4
Using equation 3 on [[media:Nm1.s11.mtg29.djvu|Lecture page 7-5]] for $$T_k \ $$ and the MATLAB code written below, the following Romberg Integration Table was generated.

Note that $$T_4(4) \ $$ results in an error of 1.3323e-015. Compared to $$T_0(64) \ $$, which has an error of 3.742e-006, $$T_4(4) \ $$ is approximately $$10^9$$ times better than the uncorrected trapezoidal rule at n=64.

Problem 6: Compute Integral until error order of magnitude reaches 10^-6
 Solved without assistance 

Reference [[media:Nm1.s11.mtg7.djvu|HW2.4 from Lecture 7-3]] and [[media:Nm1.s11.mtg30.djvu|HW5.4 from Lecture 29-6]] Compute $$ I_n \ $$ using $$ CT_k(n) \ $$, for $$ k=1,2,3 \ $$ for $$ n=2,4,8,16... \ $$ until error $$ I-I_n= O(10^{-6}) \ $$

Solution

From [[media:Nm1.s11.mtg30.djvu|Lecture 30-1]], the corrected trapezoidal rule states that

where

The following MATLAB code was developed to evaluate the corrected trapezoidal rule using $$ k=1,2,3 \ $$ until the error $$ I - I_n^T \leq O(10^{-6}) \ $$