# University of Florida/Egm6341/s10.Team2/HW4
## Problem 1: Development of higher order Corrected Trapezoidal rules and application
#### Statement

P. 20-1
Develop higher-order corrected trapezoidal rules for $CT_{1}(n)\quad CT_{2}(n)\quad CT_{3}(n)$
Then use these new corrected trapezoidal rules to find:
$I=\int _{0}^{1}{\frac {\exp ^{x}-1}{x}}dx$
until the error is of order of $10^{-6}\quad$ using $\;n=2,4,8,16,32...$

#### Solution

Higher Order corrected trapezoidal rules are of the form:
$\;CT_{k}(n)=CT_{k-1}(n)+a_{k}h^{2*k}$
where the a coefficient is defined as follows:
$a_{i}=d_{i}\left[f^{(2i-1)}(b)-f^{(2i-1)}(a)\right]$
and $di={\frac {-B_{2i}}{(2i)!}}$
B are Bernoulli numbers. The first values for $d_{i}\quad$ are:
$\;d_{1}={\frac {-1}{12}}\qquad d_{1}={\frac {1}{720}}\qquad d_{3}={\frac {-1}{30240}}$
These corrected trapezoidal rules define the integral as follows:
$\;I=CT_{k}(n)+\Theta (h^{(2*(k+1)})$
This means that higher order corrected trapezoidal rules can be developed as long as the function is continuous and its higher order derivatives can be calculated.
The first three corrected trapezoidal rules are found as follows:
$\;CT_{1}(n)=CT_{0}(n)+a_{1}h^{2}$
$\;CT_{2}(n)=CT_{1}(n)+a_{2}h^{4}$
$\;CT_{3}(n)=CT_{2}(n)+a_{2}h^{6}$
Using matlab to calculate the integrals, yields the results as follows:
It can be noticed that the corrected trapezoidal rules yield a solution to the integral much faster than the uncorrected Composite Trapezoidal Rule.

#### Matlab Code Used for Problem

MATLAB Code: Finding the Corrected Trapezoidal Rules
MATLAB Code: Calculating the Values of the function
MATLAB Code: Calculating the Composite Trapezoidal Rule
solved by--- Egm6341.s10.Team2.GV 21:56, 3 March 2010 (UTC)
checked by---Egm6341.s10.team2.patodon 20:53, 3 March 2010 (UTC)

## Problem 2: Higher order error analysis of Trapezoidal rule for periodic functions

    For which n does $\;E_{n}^{1}\simeq 0$

#### Statement

To prove that the value of the trapezoidal rule error: $\;E_{n}^{1}\simeq 0$ for a periodic function and find the effect of $\;n$(number of intervals or panels in trapezoidal rule) on the convergence of the function.

#### Solution

Let $f\left(x\right)$ be a periodic function in the interval $\left[a,b\right]$
We know that for a periodic function with periodicity of $\;b-a,$ $\;f\left(a\right)=f\left(a+\left(b-a\right)\right)=f\left(b\right)$
Differentiating $\;f\left(x\right)$ w.r.t $\;x$ :
$f'\left(a\right)=f'\left(a+\left(b-a\right)\right)=f'\left(b\right)$
similarly,
$f^{k}\left(a\right)=f^{k}\left(a+\left(b-a\right)\right)=f^{k}\left(b\right)$
Therefore the coefficients of the expansion of $\;E_{n}^{1}$ i.e. $\;a_{1},a_{2},a_{3},a_{4}$ becomes zero due the following reason:
We know that:
$a_{i}=d_{i}\times \left[f^{\left(2i-1\right)}\left(b\right)-f^{\left(2i-1\right)}\left(a\right)\right]$ where $\;d_{i}={\frac {-B_{2i}}{\left(2i\right)!}}$
$f^{\left(2i-1\right)}\left(b\right)\simeq f^{\left(2i-1\right)}\left(a\right)$ when periodicity of the function is $\;b-a$
Therefore, in the case of periodic function $f\left(x\right)$ with periodicity $\;b-a$, $\;a_{i}\simeq 0\;\forall \;i=1,2,3...$
Hence the error expressed as:
$\;E_{n}^{1}=\sum _{i=1}^{\infty }a_{i}h^{2i}\simeq 0$ as $\;a_{i}\simeq 0\;\forall \;i=1,2,3...$
The same result is tested by using a Matlab program for trapezoidal rule and it proves to be true that the function converges rapidly for periodic functions using trapezoidal rule as the number of intervals 'n' increases and thus the error term $\;E_{n}^{1}$ tends to zero.
solved by-----Srikanth Madala
checked by---Egm6341.s10.team2.patodon 20:55, 3 March 2010 (UTC)

## Problem 3: Discussions of pros and cons of various quadrature methods
#### Statement

Discuss the pros and cons of the following quadrature methods:
* Taylor Series
* Composite Trapezoidal rule
* Composite Simpson's rule
* Romberg table (including Richardson's extrapolation)
* Corrected Trapezoidal rule

#### Solution

Taylor Series:
Pros:
* 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)$
Cons:
* The function $\;f(x)$ should be differentiable in the given limits 'a' and 'b' in which it is approximated as a polynomial.
* The product of two sub-functions in a given function $\;f(x)$ like for example $e^{x}.\cos \left(x\right)$ would make it more difficult and more costly (in terms of computation) to evaluate successive derivatives.
Composite Trapezoidal rule:
Pros:
* Very intuitive and easy to use
* The numerical integral results obtained for periodic functions like trigonometric functions (sinx, cosx, tanx etc) is very accurate owing to the fact that the error is dependent on the difference between the odd derivatives of the limiting points - 'a' and 'b'.
* It yields better accuracy than simple trapezoidal rule in calculating the numerical integral of a function.
Pros:
* Needs many subintervals to converge and the convergence rate is low for non-periodic functions.
* Runge Phenomena.
Composite Simpson's rule:
Pros:
* It is still intuitive and easy to use
* The convergence rate is OK
* Exact for cubic polynomials
Cons:
* 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...)
* Runge Phenomena
Romberg table (including Richardson's extrapolation):
Pros:
* The successive computation of $\;I(2n)$ is cheaper and faster when $\;I(n)$ is already computed.
* 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.
Cons:
* 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 The function $\;x\to x^{\frac {1}{3}}$ is not differentiable at $\;x=0$, so the required conditions are not satisfied for any extrapolation.
solved by --- --Srikanth Madala
checked by --- --Egm6341.s10.Team2.GV 22:00, 3 March 2010 (UTC)

## Problem 4: Transformation of variables in the error expression of HO trapezoidal rule
#### Statement

Given the expression
    $\displaystyle {\boldsymbol {E_{n}^{1}=\sum _{k=0}^{n-1}[\int _{x_{k}}^{x_{k+1}}f(x)dx-{\frac {h}{2}}[f(x_{k})+f(x_{k+1})]]}}$ (1)
where ${\boldsymbol {x_{k}=a+kh}}$ and ${\boldsymbol {h={\frac {(b-a)}{n}}}}$,
Show that by transformation of variables we get
    $\displaystyle {\boldsymbol {E_{n}^{1}={\frac {h}{2}}\sum _{k=0}^{n-1}[\int _{-1}^{1}g_{k}(t)dt-[g_{k}(-1)+g_{k}({+1})]]}}$ (2)
where ${\boldsymbol {g_{k}(t)=f(x(t));x\epsilon [x_{k},x_{k+1}]}}$

#### Solution

Given the expression
    $\displaystyle {\boldsymbol {E_{n}^{1}=\sum _{k=0}^{n-1}[\int _{x_{k}}^{x_{k+1}}f(x)dx-{\frac {h}{2}}[f(x_{k})+f(x_{k+1})]]}}$ (1)
we use the transformation $\displaystyle {x_{k}=a+kh}$ to transform the interval from $\displaystyle {[x_{k},x_{k+1}]}$ to $\displaystyle {[-1,1]}$ by introducing $\displaystyle {[x_{k},x_{k+1}]}$.We get
$\displaystyle {x(-1)=x_{k}}$
$x(0)={\frac {x_{k}+x_{k+1}}{2}}$
$\displaystyle {x(1)=x_{k+1}}$
substituting in (1) we get
    $\displaystyle E_{n}^{1}=\sum _{k=0}^{n-1}[\int _{-1}^{+1}f(x(t))dt-{\frac {h}{2}}[f(x(-1))+f(x(+1)]]$
    $\displaystyle E_{n}^{1}=\sum _{k=0}^{n-1}[\int _{-1}^{+1}f(t{\frac {h}{2}}+0)dt-{\frac {h}{2}}[f(x(-1))+f(x(+1)]]$
    $\displaystyle E_{n}^{1}=\sum _{k=0}^{n-1}[\int _{-1}^{+1}f(t{\frac {h}{2}}+0)dt-{\frac {h}{2}}[f(x(-1))+f(x(+1)]]$
taking $\displaystyle {\frac {h}{2}}$ out we get
    $\displaystyle {\boldsymbol {E_{n}^{1}={\frac {h}{2}}\sum _{k=0}^{n-1}[\int _{-1}^{1}g_{k}(t)dt-[g_{k}(-1)+g_{k}({+1})]]}}$ (2)
solved by---Egm6341.s10.team2.niki 17:46, 27 February 2010 (UTC)
checked by--- --Srikanth Madala

## Problem 5: Integration by parts
#### Statement

Show the following relation via integration by parts
$\int _{-1}^{1}(-t)\cdot g^{(1)}(t)=\int _{-1}^{1}g(t)dt-\left[g(-1)+g(1)\right]$

#### Solution

Integration by Parts
${\begin{aligned}\int _{-1}^{1}g(t)dt&=\int _{-1}^{1}g(t)\cdot (1)dt\\&=\left[g(t)\cdot (t)\right]_{-1}^{1}-\int _{-1}^{1}t\cdot g^{(1)}(t)\\&=\left[g(-1)+g(1)\right]+\int _{-1}^{1}(-t)\cdot g^{(1)}(t)\\\end{aligned}}$
Rearrangement yields
$\int _{-1}^{1}(-t)\cdot g^{(1)}(t)=\int _{-1}^{1}g(t)dt-\left[g(-1)+g(1)\right]$
solved by---Egm6341.s10.team2.lee 01:45, 2 March 2010 (UTC)
checked by---Egm6341.s10.team2.niki 05:23, 2 March 2010 (UTC)
checked by---Egm6341.s10.team2.patodon 20:48, 3 March 2010 (UTC)

## Problem 6: Derivative of $g_{k}(t)$ as it relates to the proof to the theorem of higher order trapezoidal rules
#### Statement

The following is defined:
$g_{k}(t)=f(x(t))$
where,
$x(t)=t{\frac {h}{2}}+{\frac {x_{k}+x_{k+1}}{2}}$
Proof the following:
$g_{k}^{i}=\left({\frac {h}{2}}\right)^{i}f^{i}(x(t))$

#### Solution

The first step is to state the chain rule for a composite function which is the case in this situation,
$Y=F(G(x))$
${\frac {dY}{dx}}=F^{'}(G(x))G^{'}(x)$
The first derivative is found as follows:
$g_{k}(t)=f(x(t))$
$g_{k}^{'}(t)=f^{'}(x(t))x^{'}(t)$
$x^{'}(t)={\frac {h}{2}}$
By this the first derivative is defined as:
$g_{k}^{(1)}(t)={\frac {h}{2}}f^{'}(x(t))$
Using the chain rule again to find the second derivative as follows:
$g_{k}^{(2)}(t)={\frac {h}{2}}\left[{\frac {d}{dt}}f^{1}(x(t))\right]$
$g_{k}^{(2)}(t)=\left({\frac {h}{2}}\right)^{2}f^{(2)}(x(t))$
Apply the chain rule once more to find the third derivative as follows:
$g_{k}^{(3)}(t)=\left({\frac {h}{2}}\right)^{2}{\frac {d}{dt}}f^{(2)}(x(t))$
$g_{k}^{(3)}(t)=\left({\frac {h}{2}}\right)^{3}f^{(3)}(x(t))$
One notices that the value of ${\frac {h}{2}}$ comes out of the derivative everytime it is applied.
One is able to then conclude a general expression for the derivative of the function as follows:
$g_{k}^{(i)}(t)=\left({\frac {h}{2}}\right)^{i}f^{(i)}(x(t))$
solved by---Egm6341.s10.Team2.GV 05:03, 27 February 2010 (UTC)
checked by---Egm6341.s10.team2.patodon 21:10, 3 March 2010 (UTC)

## Problem 7: Verifying the Error for Higher Order Trapezoidal Rule
#### Statement

Show the following relation via integration by parts
$\int _{-1}^{1}g^{(1)}(t)\cdot P_{1}(t)dt=\left[P_{2(t)}\cdot g^{(1)}(t)\right]_{-1}^{1}-\int _{-1}^{1}P_{2}(t)\cdot g^{(2)}(t)dt$

#### Solution

Integration by Parts
$\int _{-1}^{1}g^{(1)}(t)\cdot P_{1}(t)dt=\left[P_{2(t)}\cdot g^{(1)}(t)\right]_{-1}^{1}-\int _{-1}^{1}P_{2}(t)\cdot g^{(2)}(t)dt$
${\mbox{where,}}\qquad P_{1}(t)=-t\qquad {\mbox{and}}\qquad P_{2}(t)=\int P_{1}(t)dt=-{\frac {t^{2}}{2}}+c$
solved by---Egm6341.s10.team2.patodon 20:59, 3 March 2010 (UTC)
checked by---Egm6341.s10.team2.lee 02:15, 2 March 2010 (UTC)

## Problem 8: Installation of chebfun program
#### Statement

Install Chebfun program (Tefethen et al. 2004)

#### Solution

Status - Successful
solved by---
checked by---Egm6341.s10.team2.patodon 21:27, 3 March 2010 (UTC)

## Problem 9: Comparison of codes
#### Statement

For the integral $\displaystyle {I=\int _{0}^{1}{\frac {e^{x}-1}{x}}dx}$
* Optimise the code for Romberg table
* Compare second column of Romberg table with the results of Simpson's rule and derive any relationship.
* Compute $\displaystyle {I_{n}}$ to the order of $\displaystyle {10^{-10}}$ and find the computational time using
* Code for Composite Trapezoidal Rule
* Code for Composite Simpson's Rule
* Code for Romberg table
* Chebfun "sum" command
* Built-in MATLAB function "trapz"
* Built-in MATLAB function "quad"
* Built-in MATLAB function "clencurt"

#### Solution

Exact value of the Integral from Mathematica is calculated by teh command "Integrate[((exp(x)-1)/x),+{x,+0,+1}]" the value is
1.3179021514544038948600088442492318379749012457927839928404611969976461077561394826119536468343922075 [[((exp(x)-1)/x),+{x,+0,+1}]]

##### Part1:Corrected code for Romberg table

MATLAB Code: For Efficient Production of the Romberg Table

##### Part2:Comaprison of Simpson's rule and Romberg table

The below table shows the values for the given integral computed using Composite Simpsons rule and the second column of the Romberg table for the same. By comparing the values of the Integral (or of the corresponding errors)) for any n, we see that the value from the composite Simpson's rule is the same as the value for n/2 (previous row) of $\displaystyle {T_{1}(n)}$.
For e.g: Denoting Simpson's as $\displaystyle {S(n)}$ we see that
    $\displaystyle \displaystyle {S(16)=T_{1}(8)}$
    $\displaystyle \displaystyle {S(32)=T_{1}(16)}$
    $\displaystyle \displaystyle {S(64)=T_{1}(32)}$
etc. In general we can write
    $\displaystyle \displaystyle {S(2^{i})=T_{1}(2^{i-1})}$
wkt
    $\displaystyle T_{1}(2^{i})={\frac {4T(2^{i})-T(2^{i-1})}{3}}$
    $\displaystyle \therefore S(2^{i})={\frac {4T(2^{i})-T(2^{i-1})}{3}}$ (1)

##### Part3:Codes

Using Code for Composite trapezoidal rule
Gives n= 32768 ,In =1.31790215149321 and time 4.973 sec
Using Code for Composite Simpson's Rule
Gives n= 128 ,In= 1.31790215146089 and time=0.017sec
Using code for Romberg table
Gives time =0.04sec
Using chebfun - sum:
I = 1.317902151454404 and time = 0.009175 sec.
Using Trapz function:
n=32768 ,In=1.31790215149321 and time = 0.006879 sec
Using Quad function
Gives n=2 ,In= 1.31790215145441 and time = 0.004 sec
Using clencurt function
Gives II=1.317902151454404 and time = 0.000308 sec

##### Part3:Computational times

solved by---Egm6341.s10.team2.niki 18:20, 28 February 2010 (UTC)
checked by---Egm6341.s10.team2.lee 22:20, 2 March 2010 (UTC)

## Problem 10: Another comparison of codes
#### Statement

Repeat Problem 9 for integral $\displaystyle {I=\int _{-5}^{5}{\frac {1}{1+x^{2}}}dx}$

#### Solution

The exact value of the definite integral is atan(5)-atan(-5)=2.746801533890034...

##### own code - composite trapezoidal rule

Function Definition
Call Function
Result

##### own code - composite simpson's rule

Function Definition
Call Function
Result

##### Romberg quadrature rule

Function Definition
Call Function
Result

##### matlab - trapz and quad

Call Function
Result

##### matlab - clencurt

Call Function
Result

##### "sum" command of chebfun matlab toolbox

Call Function
Result
solved by---Egm6341.s10.team2.lee 22:22, 2 March 2010 (UTC)
checked by-----Egm6341.s10.team2.niki 20:45, 6 March 2010 (UTC)

## Problem 11: Finding The Arc Length of an Ellipse using Numerical Integration
#### Statement

Find the Arc Length of the ellipse where:
$r(\theta )={\frac {1-e^{2}}{1-e\cos(\theta )}}$
and
$eccentricity=e=\sin({\frac {\pi }{12}})$
Using:
1)Composite Trapezoidal Rule
2)Composite Simpson's Rule
3)Romberg Table
4)"Sum" Command of the chebfun toolbox
5)"Trapz" Command of Matlab
6)"Quad" Command of Matlab
7)"Clencurt" Command of Matlab
Until the error is less than $10^{-6}$

#### Solution
##### Composite Simpson's Rule

MATLAB Code:

##### Romberg Table

Romberg Table for n=2
Table for Error of Romberg Table
Romberg Table for n=4
Error for Romberg Table n=4
Romberg Table for n=8
Table for Error when n=8

##### "Sum" Command of Chebfun Matlab Toolbox

$In=6.069090959564776$
$Error=I-In=-8.457248767967940e-008$

##### "Trapz" Command of Matlab

'''Trapz

##### "Quad" Command of Matlab

The value found using the "Quad" Matlab Command was:
$I=6.069090875$
with an error smaller than $10^{-6}$
solved by--- --Egm6341.s10.Team2.GV 22:04, 3 March 2010 (UTC)
checked by---Egm6341.s10.team2.patodon 21:15, 3 March 2010 (UTC)

## Contributing Authors and Assignments
#### Signatures

solved problems 4 and 9, proofread problems 5 and 10:--Egm6341.s10.team2.niki 18:22, 28 February 2010 (UTC)
solved problems 5 and 10, proofread problems 7 and 9:--Egm6341.s10.team2.lee 22:26, 2 March 2010 (UTC)
solved problems 7, proofread problems 1, 2, 5, 6, 8 and 11 Egm6341.s10.team2.patodon 23:01, 2 March 2010 (UTC)
solved problems 1, 6 and 11. Proofread 3, 9. --Egm6341.s10.Team2.GV 22:06, 3 March 2010 (UTC)
solved problems 3 and 8. Proofread 4, 9--Srikanth Madala