User:EGM6341.s11.TEAM1.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 State Space Model
Referenced the open-source code cauchyrnd on the MathWorks File Exchange

Given: Linear State Space Model
Given the following:

Find: Model plot with various random noise
A.) Run model and plot $$\underline{x} _0 \ $$ B.) Find the 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} } $$ C.) Use MATLAB randn to generate $$ \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 $$ D.) Run model with Cauchy random noise and plot $$\underline{x} _0 \ $$

Run linear state space model
The following MATLAB (v2009b) 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
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
The following MATLAB (v2009b) 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
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 (v2009b) 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 3: Cauchy Heavy Tails
 Solved without assistance 

Given Quartile points: $$ Q_1 \ $$, $$ Q_2 \ $$ and $$ Q_3 \ $$, where

Find the first and third Quartile points for Cauchy distribution with parameter ( x0, gamma )
Note: Culmulative distribution function(CDF) of Cauchy distribution:

Plugging Q1 into CDF function and value of 1/4 yields,

Plugging Q3 into CDF function and value of 3/4 yields,


 * Nm1 s11 team1 HW5 p3 fig1.png

Find the first and third Quartile points for Normal distribution with parameter ( u, sigma )
Note: Cumulative distribution function (CDF):

Plugging x=Q1 into the CDF function,

Multiplying 2 on both side of the eqn,

Subtracting 1 from the both side,

Tanking inverse erf function on both side of the eqn,

Multiplying sqrt(2) and sigma from the both side,

Adding mu on the both side yields,

Similarlly,

Multiplying 2 on both side of the eqn,

Subtracting 1 from the both side,

Tanking inverse erf function on both side of the eqn,

Multiplying sqrt(2) and sigma from the both side,

Adding mu on the both side yields,


 * Nm1 s11 team1 HW5 p3 fig1.png

Find appropriate sigma^1 such that gamma^G=1 and plot C(0,1) and N(0, sigma^1)
Given:

Where $$ \gamma\ ^c = $$ half width of $$ C \left ( x_0, \gamma\ ^c \right ) $$ and $$ \gamma\ ^G = $$ half width of $$ N \left ( \mu\, \sigma\ \right ) $$

Result Matlab


 * Nm1 s11 team1 HW5 p3 fig2.png

==== Find $$ \left \{ Q_1^c, Q_3^c \right \}  $$ for $$ C \left ( x_0, \gamma\ \right ) $$ and $$ C \left ( 0, 1 \right ) $$ also find $$ \left \{ Q_1^G,Q_3^G  \right \}  $$ for $$ N \left ( \mu\, \sigma\ \right ) $$ and $$ N \left ( 0, \sigma\ ^1 \right ) $$ ====

Note:  General form of $$\displaystyle \{\ Q_1^c, Q_3^c \}\ $$ for $$\displaystyle C \left( x_0, \gamma \right) $$ is already derived in the Part.A as follows,

Simply, plugging the given parameters(0,1) into the general form is the best way to find $$\{\ Q_1^c, Q_3^c \}\ $$ for $$\displaystyle C \left( 0, 1 \right) $$

Note:  General form of $$\displaystyle \{\ Q_1^G, Q_3^G \}\ $$ for $$\displaystyle N \left( x_0, \gamma \right) $$ is already derived in the Part.B as follows,

Plugging the parameters N(0,sigma) we got in part.3 into the general form,

Plot $$ \left \{ Q_1^c, Q_3^c \right \}  $$ for  $$ C \left ( 0, 1 \right ) $$ also plot $$ \left \{ Q_1^G,Q_3^G  \right \}  $$ for  $$ N \left ( 0, \sigma\ ^1 \right ) $$ and comment on results

Matlab


 * Nm1 s11 team1 HW5 p3 fig2.png
 * Nm1 s11 team1 HW5 p3 fig cdf.png

Problem 4: Spring and Dampener system
 Solved without assistance 

Find: Equation of the motion and control scheme
A.) Equation of motion B.) State space model of system C.) $$c_{cr}$$ D.) Plot $$x_k \ $$ for $$c=\frac{1}{2}c_{cr},c_{cr},\frac{3}{2}c_{cr} $$, and with Gaussian and Cauchy noise.

Equation of motion
The system(Spring-mass-damper) equation can be defined as,

Each of the elements are

Plugging each of the elements in to the eq4.1 becomes EOM(Equation of the motion) as

Define state space model
State space

Applying the given expression on p29-7, Eq(4):

Rewriting to discrete system equation, as

where I is the identity matrix, and A, B, and h are defined as above

Critical damping
Frequency and damping ratio are given as

Thus,

State Space Models with Random Noise
When external force(Control force) u=0, plot $$\displaystyle \mathbf$$ $$\displaystyle c=0.5C_{cr}$$, $$\displaystyle c=C_{cr}$$, $$\displaystyle c=1.5C_{cr}$$


 * Nm1 s11 team1 HW5 p4 fig a.png

Part.2 with Gaussian noise


 * Nm1 s11 team1 HW5 p4 fig b.png

Part3. With Cauchy noise


 * Nm1 s11 team1 HW5 p4 fig c.png

Problem 5: Modify Trapezoidal Rule MATLAB code and create Romberg Table
 Solved without assistance 

Find: Improved Trapezoidal Rule and Romberg Table
A.) Modify your Trapezoidal Rule MATLAB code to make it more efficient B.) Create Romberg Table and compare to previous results in HW2.4

More Efficient Trapezoidal Rule Code
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.

Create Romberg Table
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}) \ $$

Problem 7: Discuss the pros and cons of the following integration methods
 Solved with references listed here:  http://en.wikipedia.org/wiki/Taylor_series

http://en.wikipedia.org/wiki/Maclaurin%27s_series

http://en.wikipedia.org/wiki/Taylor_approximation

http://www.answers.com/topic/trapezium-rule

http://en.wikipedia.org/wiki/Trapezium_rule

http://en.wikiversity.org/w/index.php?title=User:Egm6341.s2010.Team1/HW4&oldid=542621

http://en.wikiversity.org/w/index.php?title=Egm6341.s10.Team2/HW4&oldid=541906

http://en.wikiversity.org/w/index.php?title=User:Egm6341.s10.team3.sa/HW4&oldid=541673

http://en.wikiversity.org/w/index.php?title=User:Egm6341.s10.Team4/HW4&oldid=541867

Taylor Series
 Positives:  1) The partial sums of the series can be used as approximations of the entire function.

2) Differentiation and integration of power series can be performed term by term and is hence particularly easy.

3) The solution is highly accurate for a small order compared to Trapezoidal and Simpson's rule.

4) It is one step and explicit.

5) Accuracy and speediness.

6) Powerful tool that helps to approximate any type of simple or complex functions into a polynomial form $$ P_n(x) \ $$ with an error in the form of remainder $$ R_{n+1}(x) \ $$

 Negatives:  1) The partial sums of the series are good if sufficiently many terms are included.

2) Problems arise with natural logarithm functions, such as log(1 + x), and some of its Taylor polynomials around a = 0. These approximations converge to the function only in the region −1 < x ≤ 1; outside of this region the higher-degree Taylor polynomials are worse approximations for the function. This is similar to Runge's phenomenon.

3) It needs the explicit form of derivatives of the function.

4) Number of operations can be large.

5) Need to know about smoothness of solution.It is senseless to use high order method when the solution has unbounded high order derivatives.

6) For oscillatory functions, a high order Taylor seems a good choice, but can lead to problems of loss of significance in the computation if not programmed carefully.

7) The function $$ F(x) \ $$ should be differentiable in the given limits 'a' and 'b' in which it is approximated as a polynomial.

8) The product of two sub-functions in a given function $$ F(x) \ $$ like for example $$ e^xcos(x) \ $$ would make it more difficult and more costly (in terms of computation) to evaluate successive derivatives.

9) It is only accurate in the area of the point that the Taylor series is calculated for (for example, $$ x=0 \ $$)

Composite Trapezoidal
 Positives:  1) The error in using the trapezium rule is approximately proportional to $$ \frac{1}{n^2} $$, so that if the number of sub-intervals is doubled, the error is reduced by a factor of 4.

2) Often converges very quickly for periodic functions

3) The method is quite simple in implementing, compared to other methods.

4) It takes a piecewise linear approximation of the function and hence execution of the solution is faster.

5) The solution is very rapidly convergent.

6) It consists of the fact that weighting coefficients are nearly equal to each other.

7) Approximate derivative by constant

8) Average (not discriminate) endpoints.

9) Very intuitive and easy to use

10) It yields better accuracy than simple trapezoidal rule in calculating the numerical integral of a function.

 Negatives:  1) Error is difficult to identify, can be eaily over or under estimated

2) The error is much higher compared to the other methods.

3) This method is less accurate for non linear functions since it has straight line approximation.

4) For various classes of functions that are not twice-differentiable,the trapezoidal rule has sharper bounds than Simpson's rule.

5) Needs many subintervals to converge and the convergence rate is low for non-periodic functions.

6) Runge Phenomena.

Composite Simpson
 Positives:  1) If the function being integrated is relatively smooth over the interval [a,b], a smooth quadratic interpolant like the one used in Simpson's rule will give good results. This is achieved with composite Simpson's by breaking up larger interval into smaller intervals.

2) It assumes piecewise quadratic interpolation and hence the accuracy is much higher compared to trapezoidal rule.

3) Exact for cubic polynomials

4) The error is less compared to that of trapezoidal rule.

5) Weighting coefficients are simple and do not fluctuate in terms of magnitude.

6) It is still intuitive and easy to use

 Negatives:  1) A large number of ordinates are needed in between the interval.

2) for various classes of rougher functions (ones with weaker smoothness conditions), the trapezoidal rule has faster convergence in general than Simpson's rule.

3) This rule can be applied to only even number of finite intervals between the limits- 'a' and 'b'. So the value of 'n' in Simpson's rule is always an even number. i.e. n=2i; where (i=1,2,3...)

4) Runge Phenomena

Romberg (Richardson)
 Positives:  1) Evaluates the integrand at equally-spaced points. The integrand must have continuous derivatives, though fairly good results may be obtained if only a few derivatives exist.

2) Newton Cotes methods fail to converge for many integrals, while Romberg integration is more stable.

3) Error is reduced subsequently as we go from $$ T_K(n) \ $$ to $$ T_{K+1}(n) \ $$.

4) It is efficient in the sense that $$ T_{K+1}(n) $$ is calculated from the already computed values of $$ T_K(n) \ $$ and $$ T_K(2n) \ $$

5) The above condition reduces the time and space requirement.

6) By using Richardson extrapolation repeatedly on the trap. rule is more efficient.

7) Converge to high accuracy.

8) The successive computation of $$ I(2n) \ $$ is cheaper and faster when $$ I(2n) \ $$ is already computed.

9) The accuracy of the numerical integration results obtained are better than other methods. The higher the power of 'h' term in the trapezoidal error, the better the accuracy.

 Negatives:  1) In Richardson extrapolation, We are neglecting the higher orders of $$ a_ih^{2i} \ $$.

2) The success of Romberg integration is only justified if the integrand-'f' satisfies the hypotheses of the Euler–Maclaurin Theorem.For example, as illustrated in the An Introduction to Numerical Analysis by Suli and Mayers Text Pg.218 (http://www.netlibrary.com/Search/SearchResults.aspx?t1=Mayers%2c+D.+F.&tt1=Author&ql=ENG%7C) The function $$ x \ $$ going to $$ x^{\frac{1}{3}}\ $$ is not differentiable at $$ x=0 \ $$, so the required conditions are not satisfied for any extrapolation.

Corrected Trapezoidal Rules: $$ cT_k(n) \ $$
 Positives:  1) When integrals of periodic functions are approximated numerically, Corrected Trapezoidal functions are the best choice.

2) More accurate than composite trapezoidal rule and Romberg table.

3) Converge more quickly than trapezoidal rule.

 Negatives:  1) This method demands derivatives of the function at the end points of the interval.

2) Bernoulli's number should be calculated for each order of derivatives.

Problem 8: The Theorem of Higher Order Trapezoidal Rule Error(HOTRE)
 Solved without assistance 

Given
,where $$\begin{align} x_k=a+kh, \; \ h=(b-a)/n, \; \ x(t)= t*\frac{h}{2} + \frac{x_i+x_{i+1}}{2}, \; \ t \epsilon [-1,1] \end{align}$$ ,where $$\begin{align} g_k(t):=f(x(t)) \; \ x\in[x_k,x_{k+1}] \end{align}$$

Show eqn.8.1 is equal to eqn.8.2
from eqn.8.3 we have $$\begin{align} x_k=-1, \; x_{k+1}=1 \end{align}$$ $$\begin{align} \Rightarrow f(x_k)=g(-1)\end{align}$$ $$\begin{align}, \; f(x_{k+1})=g(1)\end{align}$$ $$\begin{align}, \; \frac{dx(t)}{dt}=\frac{h}{2} \end{align}$$

Plugging above terms and eqn7.3 into eqn.7.1

Show HOTRE for Trapezoidal rule Step 1 is true
From lecture note [[media:Nm1.s11.mtg30.djvu|30-3]], First step of HOTRE for Trap. rule is as following; Using Integration by parts, we can rewrite LHS of eqn.8.5; $$\displaystyle\begin{align} \int u {v}^{'}= uv - \int v {u}^{'} \end{align}$$ Let $$\displaystyle\begin{align} u=t \ \ {v}^{'}=g(t) \end{align}$$

$$\displaystyle \begin{align} = \int_{-1}^{1}g(t)dt-(g(-1)+g(1)) \\ \end{align}$$

So, HOTRE for Trap. rule Step 1(eqn.8.5) is true

Show HOTRE for Trap. rule Step 2a is true
From lecture note [[media:Nm1.s11.mtg30.djvu|30-3]], Step 2a of HOTRE for Trap. rule is as following; where $$\displaystyle \begin{align} \; g_k^{(i)}(t)=\frac{d^i}{dt^i} g_k^{(t)}, \;x \epsilon[x_k, x_{k+1}] \end{align}$$ From eqn.8.1 and eqn.8.2,

=Contributing Members & Referenced Lecture=