User:Egm6341.s10.Team4/HW4

= Problem 1: Higher Order Corrected Trapezoidal Rule =

Given
Refer Lecture slide 20-1 for problem statement.

Error for Higher Order Error for Trapezoidal rule on slide 18-2, is given by:


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

$$ $$
 * $$\displaystyle E_n^1 := I - I_n
 * $$\displaystyle E_n^1 := I - I_n
 * $$\displaystyle \longrightarrow (1)
 * }.
 * }.

where $$\displaystyle I_n = T(n) $$ is the Composite Trapezoidal approximation of the Integral $$\displaystyle I = \int_{a}^{b} f(x)dx $$ with $$\displaystyle n $$ intervals and $$\displaystyle n+1 $$ points $$\displaystyle \forall x \in [a,b] $$.


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

T(n)=\frac{h}{2}(f_0+2f_1+2f_2+ \cdots +2f_{n-1}+f_n) $$ $$
 * $$ \displaystyle\
 * $$ \displaystyle\
 * $$\displaystyle \longrightarrow (2)
 * }.
 * }.

where $$\displaystyle\ f_i=f(x_i)$$.

Equation (1) is a Euler-Maclaurin series given by Refer 18-3:


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

$$ $$
 * $$\displaystyle E_n^1 = \sum_{q=1}^{\infty} a_q h^{2q}
 * $$\displaystyle E_n^1 = \sum_{q=1}^{\infty} a_q h^{2q}
 * $$\displaystyle \longrightarrow (3)
 * }.
 * }.

where $$\displaystyle h = \frac{b-a}{n} $$.

Find
1. To develop a general higher order corrected Trapezoidal Rule.

2. To find $$\displaystyle n, I_n, E_n^1 $$ using $$\displaystyle CT_1, CT_2, CT_3 $$ for which $$\displaystyle I - I_n $$ is of the order $$\displaystyle \Theta(10^{-6}) $$ for $$\displaystyle I = \int_{0}^{1} \frac{e^x-1}{x}dx $$.

Solution
1. To develop a general Higher Order Corrected Trapezoidal Rule:

Equation (3),


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

\begin{align} E_n^1 &= \sum_{q=1}^{\infty} a_q h^{2q} \\ \Rightarrow     E_n^1 &= I - T(n) = \sum_{q=1}^{\infty} a_q h^{2q} &= I - T(n) = a_1h^{2} + \sum_{q=2}^{\infty} a_q h^{2q} \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.


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

\begin{align} \Rightarrow     E_n^1 &= I - \left[\underbrace{T(n) - a_1h^{2}}_{CT_1}\right] = \sum_{q=2}^{\infty} a_q h^{2q} \\ &= I - CT_1(n) = a_2h^{4} + \sum_{q=3}^{\infty} a_q h^{2q} \\ \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.


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

\begin{align} \Rightarrow     E_n^1 &= I - \left[\underbrace{CT_1 + a_2h^{4}}_{CT_2}\right] = \sum_{q=3}^{\infty} a_q h^{2q} \\ &= I - CT_2(n) = a_3h^{6} + \sum_{q=4}^{\infty} a_q h^{2q} \\ \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.


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

\begin{align} \Rightarrow     E_n^1 &= I - \left[\underbrace{CT_2 + a_3h^{6}}_{CT_3}\right] = \sum_{q=4}^{\infty} a_q h^{2q} \\ &= I - CT_3(n) = a_4h^{8} + \sum_{q=5}^{\infty} a_q h^{2q} \\ &.\\ &.\\ &.\\ &.\\ &.\\ \Rightarrow     E_n^1 &= I - \left[\underbrace{CT_{(k-1)} + a_kh^{2k}}_{CT_k}\right] = \sum_{q=(k+1)}^{\infty} a_q h^{2q} \\ &= I - \underbrace{CT_k(n)}_{I_n} = a_{(k+1)}h^{2(k+1)} + \sum_{q=k+2}^{\infty} a_q h^{2q} \\
 * $$\displaystyle
 * $$\displaystyle

\end{align} $$
 * }.
 * }.

Hence the Higher Order Corrected Trapezoidal rule is given by,


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

$$\displaystyle
 * style="width:20%; padding:10px; border:2px solid #8888aa" |
 * style="width:20%; padding:10px; border:2px solid #8888aa" |

\Rightarrow     I_n = CT_k(n) = CT_{(k-1)}(n) + a_kh^{2k}

$$ $$
 * style= |
 * $$\displaystyle \longrightarrow (4)


 * }.
 * }.

2. To find $$\displaystyle I_n $$ using $$\displaystyle CT_1, CT_2, CT_3 $$:

From Equation (4), $$\displaystyle I_n $$ using $$\displaystyle CT_1, CT_2, CT_3 $$ can be given by:


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



\begin{align} \Rightarrow     I_n &= CT_3 \\ &= CT_2 + a_3h^{6}\\ &= CT_1 + a_2h^{4} + a_3h^{6}\\ &= T(n) + a_1h^{2} + a_2h^{4} + a_3h^{6} \end{align} $$
 * $$\displaystyle


 * }.
 * }.

Here,


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



\begin{align} a_i &= d_i \left[ f^{(2i-1)}(b) - f^{(2i-1)}(a)\right]\\ d_i &= - \frac{B_{2i}}{2i}; \ \ \ B_{2i} \rightarrow Bernoulli\ Numbers;\ \ \ d_1 = \frac{-1}{12}; \ \ \ d_2 = \frac{1}{720}; \ \ \ d_3 = \frac{-1}{30240} \end{align}
 * $$\displaystyle

$$


 * }.
 * }.

At $$\displaystyle n = 8$$, the Error is reduced to an order of $$\displaystyle \Theta(10^{-6}) $$.



 Matlab Code: 

Author
Solved and typed by - Egm6341.s10.Team4.andy 08:44, 3 March 2010 (UTC)

= Problem 2: Applying periodic function for higher-order of Trap. rule =

Given
Envisage the higher-order error for Trap. rule which has been introduced as following [[media:Egm6341.s10.mtg19.djvu|19-3]], for a function $$\displaystyle f(x) $$ which is periodic on $$\displaystyle [a,b] $$


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle E_n=a_1h^{2\times 1}+a_2h^{2\times 2}+...+a_nh^{2\times n}=\sum_{i=1}^\infty a_ih^{2i}


 * }
 * }


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle for k=0,1,2,...        f^{(k)}(a)=f^{(k)}(b)


 * }
 * }


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle a_i=d_i[f^{(2i-1)}(b)-f^{(2i-1)}(a)]


 * }
 * }

Find
Demonstrate that for which $$\displaystyle n $$, the $$\displaystyle E_n=0 $$ would be achieved?

Solution
Since $$\displaystyle f(x) $$ is a periodic function on the domain, we can use the following property of periodic functions (refer to first equation on the [[media:Egm6341.s10.mtg20.djvu|20-2]]:

If $$\displaystyle f(x) $$ is periodic on $$\displaystyle [a,b] $$, the odd derivatives of it will be equal at the same interval:


 * $$\displaystyle

f^{(k)}(a)=f^{(k)}(b) $$



k=2i-1; i\in \mathbb{N} $$



\Rightarrow f^{(2i-1)}(a)=f^{(2i-1)}(b) $$



\Rightarrow a_i=d(di=-\dfrac {B_{2i}}{(2i!)}\cancelto{0}{[f^{(2i-1)}(b)-f^{(2i-1)}(a)]}=0 $$

Thus, recalling the higher-order error of Trapezoidal rule;



E_n=\sum_{i=1}^n a_ih^{2i}=0 $$

Author
Solved and typed by - Egm6341.s10.Team4.nimaa&amp;m 21:51, 27 February 2010 (UTC) .

= Problem 3: Pros and cons of different numerical methods =

Find

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




 * Discuss the pros and cons of different numerical methods.
 * 1) Taylor's series.
 * 2) Composite Trapezoidal rule.
 * 3) Composite Simpson's rule.
 * 4) Romberg Table (Richardson Extrapolation).
 * 5) Corrected Trapezoidal rule.
 * }.
 * }.

Solution
The pros and cons of the above mentioned methods are discussed below.

Taylors Series:

Advantages:
 * 1) The solution is highly accurate for a small order compared to Trapezoidal and Simpson's rule.
 * 2) It is one step and explicit.

Disadvantages:
 * 1) It needs the explicit form of derivatives of the function.

Composite Trapezoidal rule:

Advantages:
 * 1) The method is quite simple in implementing, compared to other methods.
 * 2) It takes a piecewise linear approximation of the function and hence execution of the solution is faster.
 * 3) The solution is very rapidly convergent.
 * 4) It consists of the fact that weighting coefficients are nearly equal to each other.

Disadvantages:
 * 1) The error is much higher compared to the other methods.
 * 2) This method is less accurate for non linear functions since it has straight line approximation.

Composite Simpson's rule:

Advantages:
 * 1) It assumes piecewise quadratic interpolation and hence the accuracy is much higher compared to trapezoidal rule.
 * 2) It can integrate polynomials up to third order accurately.
 * 3) The error is less compared to that of trapezoidal rule.
 * 4) Weighting coefficients are simple and do not fluctuate in terms of magnitude.

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

Romberg Table (Richardson Extrapolation):

Advantages:
 * 1) Error is reduced subsequently as we go from TK(n) to TK+1(n).
 * 2) It is efficient in the sense that TK+1(n) is calculated from the already computed values of TK(n) and TK(2n).
 * 3) The above condition reduces the time and space requirement. Also, the table format is good for pursuing the history of calculations and making some comparisons.

Disadvantages:
 * 1) In Richardson extrapolation, We are neglecting the higher orders of aih2i.

Corrected Trapezoidal rule:

Advantages:
 * 1) When integrals of periodic functions are approximated numerically, Corrected Trapezoidal functions are the best choice. Even we can reach to the error equal to zero in this condition.

Disadvantages:
 * 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.

Author
Solved and typed by - Egm6341.s10.team4.anandankala 13:37, 3 March 2010 (UTC) Reviewed by - Egm6341.s10.Team4.nimaa&amp;m 13:57, 3 March 2010 (UTC) .

= Problem 4: Proof of Higher Order Error for Trap. rule Theorem =

Given
Higher Order Error for Trap. rule equation (1) on slide [[media:Egm6341.s10.mtg21.djvu|21-1]], which is:


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle Eq. (1) \ E_n^1= \sum_{k=0}^{n-1}[ \int_{x_k}^{x_{k+1}}f(x)dx- \frac{h}{2}(f_k(x_k)+f(x_{k+1}))]
 * }
 * }

Show
That eq. (1) can be transformed to the following: equation (5) on slide [[media:Egm6341.s10.mtg21.djvu|21-2]]


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle E_n^1=\dfrac{h}{2} \sum_{k=0}^{n-1}[ \int_{-1}^{1}g_k(t)dt-(g_k(-1)+g_k(+1))]

where, equation (2) on slide [[media:Egm6341.s10.mtg21.djvu|21-1]]

h= \frac{(b-a)}{n} \ \ \ x_i=a+ih $$ and, equation (3) on slide [[media:Egm6341.s10.mtg21.djvu|21-1]]

x(t)= t*\frac{h}{2} + \frac{x_i+x_{i+1}}{2} \ \ \ t\in[-1,1] $$
 * style= |
 * }.
 * }.

Solution
Defining: equation (1) on slide [[media:Egm6341.s10.mtg21.djvu|21-2]]

g_i(t):=f(x(t)) \ where \ x\in[x_i,x_{i+1}] $$

For one Panel, we have:



x_i=-1 \ \ and \ \ x_{i+1}=1 $$



f(x)_i=g(-1) \ \ \ f(x)_{i+1}=g(1) $$

Also,



\frac{dx(t)}{dt}=\frac{h}{2} $$



\Rightarrow E_n^1= \sum_{k=0}^{n-1} \int_{-1}^{1}g_k(t)\dfrac{h}{2}dt-\dfrac{h}{2}(g_k(-1)+g_k(+1)) $$

Giving,


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

$$\displaystyle E_n^1=\dfrac{h}{2} \sum_{k=0}^{n-1}[ \int_{-1}^{1}g_k(t)dt-(g_k(-1)+g_k(+1))] $$
 * style="width:20%; padding:10px; border:2px solid #8888aa" |
 * style="width:20%; padding:10px; border:2px solid #8888aa" |


 * style= |


 * }
 * }

Author
Solved and typed by - --Egm6341.s10.Team4.roni 07:39, 6 March 2010 (UTC) Reviewed by - --Egm6341.s10.Team4.roni 07:39, 6 March 2010 (UTC) .

= Problem 5: Theorem of Higher Order Error for Trap. rule Step1 =

Given
Higher Order Error for Trap. rule equation (5) on slide [[media:Egm6341.s10.mtg21.djvu|21-1]], which is:


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle Eq. (5)\ E_n^1=\dfrac{h}{2} \sum_{k=0}^{n-1}[ \int_{-1}^{1}g_k(t)dt-(g_k(-1)+g_k(+1))]
 * }
 * }

Find
That using eq. (1) on slide [[media:Egm6341.s10.mtg21.djvu|21-2]],


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle g_k(t):=f(x(t)) \  such \  that \  x \   \epsilon  \left[ x_k,x_{k+1}\right]

For one panel Step 1 :



Eq.(2) \int_{-1}^{1}(-t) {g}^{(1)}(t)dt = \int_{-1}^{1}g_k(t)dt-[(g_k(-1)+g_k(+1))] $$


 * style= |
 * }
 * }

Solution
Using integration by Parts :



\int u {v}^{'}= uv - \int v {u}^{'}\ \ where \ \ u=t \ \ {v}^{'}=g(t) $$

Writing the left side of Eq. (2):



\int_{-1}^{1}(-t) {g}^{(1)}(t)dt = -\int_{-1}^{1}(t) {g}^{(1)}(t)dt $$



= -\Bigg[(t \, g(t)) \Bigg|^{1}_{-1} -\int_{-1}^{1}g(t)dt \Bigg] $$



= -(1*g(1)-(-1)g(-1)-\int_{-1}^{1}g(t)dt) $$



= \int_{-1}^{1}g(t)dt)-(g(-1)+g(1)) $$


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

$$\displaystyle \int_{-1}^{1}(-t) {g}^{(1)}(t)dt = \int_{-1}^{1}g(t)dt)-(g(-1)+g(1)) $$
 * style="width:20%; padding:10px; border:2px solid #8888aa" |
 * style="width:20%; padding:10px; border:2px solid #8888aa" |


 * style= |


 * }
 * }

Author
Solved and typed by - --Egm6341.s10.Team4.roni 07:40, 6 March 2010 (UTC) Reviewed by - --Egm6341.s10.Team4.roni 07:40, 6 March 2010 (UTC)

= Problem 6: i-derivative of g function in higher-order of Trap. proof =

Given
Consider the following definitions for these two functions such that $$\displaystyle x\in [x_k,x_{k+1}] $$ interval:


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle x(t):=t\dfrac {h}{2}+\dfrac {x_k+x_{k+1}}{2}
 * }
 * }


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle g_k(t):=f(x(t))


 * }
 * }

Find
Show that $$\displaystyle g_k^{(i)}(t)=(\dfrac {h}{2})^i.f^{(i)}(x(t)) $$:

Solution
$$\displaystyle g_k(t)=f(x(t))=f(t\dfrac {h}{2}+\dfrac {x_k+x_{k+1}}{2}) $$

$$\displaystyle \Rightarrow g_k^{(1)}(t)=\dfrac {d[f(\dfrac {th}{2}+\dfrac {x_k+x_{k+1}}{2})]}{dt} $$

Using the Chain rule for derivation:

$$\displaystyle \Rightarrow g_k^{(1)}(t)=\dfrac {d[f(x(t))]}{dx}\times\dfrac {h}{2} $$

For second derivative of $$\displaystyle g_k^{(2)} $$, we have:

$$\displaystyle g_k^{(2)}(t)=\frac {d[g_k^{(1)}(t)]}{dt}=\frac {d[\frac {d[f(x(t))]}{dx}\times\frac {h}{2}]}{dt}=\frac {d[\frac {d[f(\frac{th}{2}+\frac {x_k+x_{k+1}}{2})]}{dx}\times\frac {h}{2}]}{dt} $$

Using the Chain rule for derivation:

$$\displaystyle \Rightarrow g_k^{(2)}(t)=\frac {d^2(f(x(t))}{dx^2}\times (\dfrac {h}{2})\times(\dfrac {h}{2})=\frac {d^2(f(x(t))}{dx^2}\times (\dfrac {h}{2})^2 $$

Similarly, for third derivative of $$\displaystyle g_k^{(3)} $$ we will find:

$$\displaystyle \Rightarrow g_k^{(3)}(t)=\frac {d^3(f(x(t))}{dx^3}\times (\dfrac {h}{2})\times (\dfrac {h}{2})\times (\dfrac {h}{2})=\frac {d^3(f(x(t))}{dx^3}\times (\dfrac {h}{2})^3 $$

Thus, for $$\displaystyle i^{th} $$, $$\displaystyle g_k^{(i)} $$ can be predicted as:

$$\displaystyle g_k^{(i)}(t)=(\dfrac {h}{2})^i.f^{(i)}(x(t)) $$

Author
Solved and typed by - Egm6341.s10.Team4.nimaa&amp;m 00:08, 2 March 2010 (UTC) .

= Problem 7: Theorem of Higher Order Error for Trap. rule [Step:2a - Part2] =

Given
Refer Lecture Slide Refer 21-2 for problem statement.

Higher Order Error for Trap. rule is given by [ Refer 21-1 ]:


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

E_n^1=\frac{h}{2} \sum_{k=0}^{n-1} \left[ \underbrace{\int_{-1}^{1}g_k(t)dt-(g_k(-1)+g_k(+1))}_{E}\right] $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle \longrightarrow (1)
 * }.
 * }.

where $$\displaystyle E $$ can be written as,


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

E =\int_{-1}^{1}(-t)g^{(1)}(t)dt $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle \longrightarrow (2)
 * }.
 * }.

Find
Show that,


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

$$\displaystyle E =\int_{-1}^{1}\underbrace{(-t)}_{p_{1}(t)}g^{(1)}(t)dt = \left[p_2(t)g^{(1)}(t) \right]_{-1}^{+1} - \int_{-1}^{+1}p_2(t)g^{(2)}(t)dt $$

$$
 * $$\displaystyle \longrightarrow (3)
 * }.
 * }.

where $$\displaystyle p_2(t) = \int p_1(t)dt$$.

Solution
Integrating by Parts,


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

\int u dv = uv - \int v du $$
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.

where,


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

$$\displaystyle \begin{align} u &= g^{(1)}(t) \\ \Rightarrow du &= g^{(2)}(t)dt \end{align} $$


 * }.
 * }.

and,


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

$$\displaystyle \begin{align} dv &= p_1(t)dt \\ &= -tdt \\ & \\ \Rightarrow v &= \int p_1(t)dt &=p_2(t) \end{align} $$


 * }.
 * }.

Equation (3) is then,


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

$$\displaystyle \begin{align} \Rightarrow E &=\int_{-1}^{1}\underbrace{(-t)}_{p_{1}(t)}g^{(1)}(t)dt \\ &= \left[p_2(t)g^{(1)}(t) \right]_{-1}^{+1} - \int_{-1}^{+1}p_2(t)g^{(2)}(t)dt \end{align} $$


 * }.
 * }.

Author
Solved and typed by - Egm6341.s10.Team4.andy 08:44, 3 March 2010 (UTC)

= Problem 8: Chebfun coding for Clenshaw-Curtis quad.=

Given
Codes of the Chebfun project are available on the lecture plan link of the course web site as a zip. file to apply in programming softwares. These codes are written for estimating the Clenshaw-Curtis quadrature by (Trefethen et al.)with a command in which started with (sum).

Find
Install Chebfun codes to see how this set of codes can estimate Clenshaw-Curtis quadrature.

Solution
Installation instructions: Chebfun Version 3 is compatible with MATLAB 7.4 (2007a) and more recent versions. First, contents of the zip file to a directory. The name "chebfun" for this directory is suggested. Then, In MATLAB, the chebfun directory is added to the path. This can be done by selecting "File/Set Path..." from the main or Command window menus. It is recommended to select the "Save" button on this dialog so that chebfuns are on the path automatically in future MATLAB sessions. To check whether the system appears to be working or not, "chebtest" should be typed at the command window prompt. Now, $$\displaystyle >> f=chebfun ('f(x)') $$ can be used for estimating the Clenshaw-Curtis quadrature.

The curve between the data points is the polynomial interpolant, which is evaluated by the barycentric formula introduced by Salzer [Berrut & Trefethen 2004, Salzer 1972]. This method of evaluating polynomial interpolants is stable and efficient even if the degree is in the millions [Higham 2004].

As an example if we want to know what is the integral of f from a to b? we should use this command:



>> sum ('f(x)') $$

This command can be applied for integration of $$\displaystyle f(x) $$ by Clenshaw-Curtis quadrature, as an alternative method for composite trapezoidal rule, composite simpson's rule, some Matlab commands and Richardson extrapolation (Romberg table).

Author
Solved and typed by - Egm6341.s10.Team4.nimaa&amp;m 21:59, 27 February 2010 (UTC) .

= Problem 9 : A Comparison of Integration Methods =

Given
We have been given various integration methods thus far throughout the class, we are now asked to compare their efficiency and accuracy.

Find
We are asked to

1) correct matlab code for Romberg efficiently

2) compare the second column of the Rhomberg Table (T1(n)) to the Composite Simpson's Rule, and derive any relationship that exists.

3) integrate the function $$ \displaystyle\ f(x)=\frac{e^x-1}{x}$$ over the interval $$ \displaystyle\ x \epsilon [0,1]$$

using the following methods and to compare the computational cost as a function of run-time -
 * Composite Trapezoidal Rule(optimized)
 * Composite Simpson's Rule
 * Rhomberg Interpolation
 * chebfun sum command - a function that uses Chebyshev polynomials to integrate
 * the matlab functions trapz, quad, and clencurt.

Solution
1) correct matlab code for Romberg efficiently

The composit Trapozoida rule for $$\displaystyle n=2^j$$ and Romberg table was computed efficiently using this equation. (the modified matlab codes are in below)


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

T_0(2^{j})= \frac{1}{2}T_o(2^{j-1})+ h \sum_{i=1}^{2^{j-1}}f(x_{2i-1}) $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle \longrightarrow (1)

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

h = \frac{b-a}{2^{j}} $$
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.

[reference] derivation of equation(1)
 * {| style="width:100%" border="0" align="left"


 * style="width:20%; padding:10px; border:2px solid #8888aa" |
 * {| style="width:100%" border="0" align="left"
 * {| style="width:100%" border="0" align="left"

T_o(2^{j})= \frac{b-a}{2^{j}} [ \frac{1}{2}f(x_0)+f(x_1)+f(x_2)+f(x_3)+f(x_4)+......+ \frac{1}{2} f(x_{2^{j}})] $$ For example
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.


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

T_0(64)= \frac{b-a}{64} [ \frac{1}{2}f(x_0)+f(x_1)+f(x_2)+f(x_3)+f(x_4)+......f(x_{63}) + \frac{1}{2} f(x_{64})] $$ the rear part can be arranged like that (decouple to the even number & odd number)
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.
 * {| style="width:100%" border="0" align="left"

T_0(64)= \frac{b-a}{64} [ \underbrace{\frac{1}{2}f(x_0)+f(x_2)+f(x_4)+... \frac{1}{2} f(x_{64})}_{even} ~ + \underbrace{f(x_1) +f(x_3)+f(x_5)+......f(x_{63})}_{odd}] $$ The even number part is exactly same with $$\displaystyle \frac{1}{2}T_0(32) $$
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.


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

\frac{b-a}{64}[ \frac{1}{2}f(x_0)+f(x_2)+f(x_4)+... \frac{1}{2} f(x_{64})] = \frac{1}{2} * \frac{b-a}{32}[ \frac{1}{2}f(x_0)+f(x_1)+f(x_2)+... \frac{1}{2} f(x_{32})] = \frac{1}{2}T_o(32)$$ So, the orginal $$\displaystyle T_0(64) $$ is like that
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.


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

T_0(64)= \frac{1}{2}T_o(32) + \frac{b-a}{64} [f(x_1) +f(x_3)+f(x_5)+......f(x_{63})] = \frac{1}{2}T_o(32) + \frac{b-a}{64} \sum_{i=1}^{32}f(x_{2i-1}) $$ if this is generalized,
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.


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

T_0(2^{j})= \frac{1}{2}T_o(2^{j-1})+ h \sum_{i=1}^{2^{j-1}}f(x_{2i-1}) $$ where
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.
 * {| style="width:100%" border="0" align="left"

h = \frac{b-a}{2^{j}} $$
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * style= |
 * }.
 * }.
 * }.

2) compare the second column of the Rhomberg Table (T1(n)) to the Composite Simpson's Rule, and derive any relationship that exists.

For the first three methods, pre-existing codes written by the authors were used in order to evaluate the computational time required to perform the integration. For the Rhomberg Extrapolation, the computation time was determine by the time required to compute T0(n) plus the time to perform the extrapolation.

When the Simpson's Rule and T1 (Richardsen Extrapolation) are compared to each other, the results are identical out to all of the significant digits calculated, thus we believe that the two are equivalents processes. Furthermore,

$$ \displaystyle T_1(2n)=\frac{4T_0(2n)-T_0(n)}{4-1}$$

Breaking this into $$ \displaystyle T_0(n)=h(\frac{1}{2}f_0+f_1+f_2+...\frac{1}{2}f_n$$ and $$ \displaystyle T_0(n/2)=2h(\frac{1}{2}f_0+f_2+f_4+...\frac{1}{2}f_n$$ we can see that $$ \displaystyle T_1(n)=\frac{h}{3}(f_0+4f_1+2f_2+4f_3+...f_n)$$, which is equal to the Simpson's Rule for n points (where h is twice the spacing between two points).

For the Trapezoidal Rule and Rhomberg Extrapolation

For the Simpson's Rule

For Matlab functions and chebshev

For ClenCurt Integration

Author
Egm6341.s10.Team4.riherd 14:29, 3 March 2010 (UTC) --Egm6341.s10.Team4.roni 07:41, 6 March 2010 (UTC) Proof read by Egm6341.s10.Team4.yunseok 22:45, 7 March 2010 (UTC)

= Problem 10 - Another comparison of integration methods =

Find
We are asked to integrate the function $$ \displaystyle\ f(x)=\frac{1}{1+x^2}$$ over the interval $$ \displaystyle\ x \epsilon [-5,5]$$ using the following methods and to compare the computational cost as a function of run-time -
 * Composite Trapezoidal (optimized)
 * Composite Simpson's Rule
 * Rhomberg Interpolation
 * chebfun sum command - a function that uses Chebyshev polynomials to integrate
 * the matlab functions trapz, quad, and clencurt.

Solution
Operating in the same manner as above, only for a different $$ \displaystyle f(x)$$ and bounds of integration.

For the Trapezoidal and Rhomberg Extrapolation

For the Composite Simpson's Rule

For the Matlab functions

For ClenCurt Integration

Author
Solved and Written by Egm6341.s10.Team4.riherd 16:06, 3 March 2010 (UTC)

ClenCurt Solved and Written by --Egm6341.s10.Team4.roni 07:56, 6 March 2010 (UTC)

Proof read by Egm6341.s10.Team4.yunseok 22:44, 7 March 2010 (UTC)

= Problem 11: Application of Engineering Orbital Mech. =

Given
Refer Lecture slide 25-2 for problem statement

Polar form relative to focus
 * {| style="width:100%" border="0" align="left"

r(\theta) = \frac{1-e^2}{1-e \, cos(\theta)} $$
 * [[Image:Ellipse Polar.svg|200px]]
 * $$\displaystyle
 * $$\displaystyle
 * }.
 * }.

where, eccentricity $$\displaystyle e = \sin(\frac{\pi}{12}) $$

Find
Find Arc length of ellipse (orbital)using below integration method.

compute $$\displaystyle I_n$$ to the error $$\displaystyle 10^{-10}$$ and

compute time using tic/toc matlab commend

1) Composit Trapozoidal rule

2) Composit Simpson's rule

3) Romberg Table

4) Chebfun, sum commend

5) matlab commend : Trapz, quad, clencurt

Solution
At first, we tried to calculate the arc length of ellipse using the $$\displaystyle r(\theta) = \frac{1-e^2}{1-e \, cos(\theta)}$$ however, we couldn't find any equation or relationship between arc length of ellipse and $$\displaystyle r(\theta) $$.

so, we found another method to compute the circumference of ellipse in wikipidia.

The circumference $$C$$ of an ellipse is $$\displaystyle 4 a E(e)$$, where the function $$E$$ is the complete elliptic integral of the second kind. Here, a=1:

The complete elliptic integral of the second kind E is defined as


 * $$E(e) = \int_0^{\pi/2}\sqrt {1-e^2 \sin^2\theta}\ d\theta\!$$

So,


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

I= 4\int_0^{\pi/2}\sqrt {1-e^2 \sin^2\theta}\ d\theta $$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle \longrightarrow (1)
 * }.
 * }.

using the equation (1), the arc length of ellipse was calculated by several numerical method.

The result is summarized.

 

Ramanujan approximation is 's:
 * $$C \approx \pi \left[3(a+b) - \sqrt{(3a+b)(a+3b)}\right]= \pi(3(a+b)-\sqrt{10ab+3(a^2+b^2)})$$

where, a = longest radius, b= shortest radius in ellipse

the resuts of our integration method and Ramanujan approximation were same. $$\displaystyle (c \simeq 6.1766)$$

(1) Composit Trapoziodal rule

 Matlab Code: 

 subfunction ff(x) 

(2) Composit Simpson's rule

 Matlab Code: 

(3) Romberg Table

 Matlab Code: 

(4) Chebfun

 Matlab Code: 

(5) Trapz (matlab commend)

 Matlab Code: 

(6) quad (matlab commend)

 Matlab Code: 

(7) ClenCurt Integration)

 Matlab Code: 

ClenCurt Integration Results -

.

Author
Solved and typed by - Egm6341.s10.Team4.yunseok 16:38, 3 March 2010 (UTC)

ClenCurt Solved and typed by--Egm6341.s10.Team4.roni 07:58, 6 March 2010 (UTC) .

Proof read by - Egm6341.s10.Team4.riherd 17:51, 3 March 2010 (UTC)

And - Egm6341.s10.Team4.nimaa&amp;m 20:52, 7 March 2010 (UTC)

= Contributing Members = Egm6341.s10.Team4.andy 08:55, 3 March 2010 (UTC)- Author of Problems:1, 7 ; Proof Read: 2,4,5,6

Egm6341.s10.Team4.riherd 16:07, 3 March 2010 (UTC) - Author of Problems 9,10; Proof Read:11

Egm6341.s10.team4.anandankala 13:37, 3 March 2010 (UTC)- Author of problems 3; Proof Read:1,7

Egm6341.s10.Team4.yunseok 16:38, 3 March 2010 (UTC)

Egm6341.s10.Team4.nimaa&amp;m 22:15, 3 March 2010 (UTC)- Author of problems 2, 6, 8; Proof Read:3,11

Egm6341.s10.Team4.roni 07:58, 6 March 2010 (UTC)