User:Egm6341.s10.team3.sa/HW4

HW Poblem Set 4

= Proof: Trapezoidal rule converges rapidly for periodic functions = Ref: Lecture Notes [[media:Egm6341.s10.mtg20.djvu|p.20-2]]

Problem Statement
The error in the higher order trapezoidal rule is given by,
 * {| style="width:100%" border="0" align="left"

E_n^{1} = \sum_{i=1}^{\infty}a_i \cdot h^{2i} $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$\displaystyle (Eq. 1)$$
 * }

For which value of n is $$\displaystyle E_n^{1} = 0$$

Solution
In Eq.(1) the coefficients $$a_i$$ are given by,
 * {| style="width:100%" border="0" align="left"

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

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

$$ \displaystyle f\,\,\,\,Periodic\,\,\,\,function $$ $$ \displaystyle a \,\,and\,\,b\,\,\,\,end\,\,\, points\,\,\, of\,\,\, the\,\,\, function\,\,\,f.$$
 * $$ \displaystyle d_i\,\,\,constant $$
 * $$ \displaystyle d_i\,\,\,constant $$
 * }

When f is periodic in [a b] $$ \displaystyle f^{k}(a) = f^{k}(b) \,\,\,\, ,k=0,1,2........ $$

Hence , $$ \displaystyle f^{2i-1}(a) = f^{2i-1}(b) \,\,\,\, \forall \,\, i$$. $$\therefore \displaystyle a_i = 0 \forall \,\,\, i $$ $$\therefore \displaystyle E_n^1 = 0 $$ irrespective of the value of n.

This being Trapezoidal rule $$\displaystyle n\geq 2. $$

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

$$ \displaystyle E_n^1 = 0 \,\,\, \forall \,\,\, n\geq 2 $$
 * style="width:25%; padding:10px; border:2px solid #8888aa" |
 * style="width:25%; padding:10px; border:2px solid #8888aa" |
 * style = |
 * }

Subramanian Annamalai 18:21, 1 March 2010 (UTC)

= Discuss Pros and cons of different Numerical Methods =

Ref. Lecture notes [[media:Egm6341.s10.mtg20.djvu|p.20-2]]

Problem Statement
Discuss Pros and cons of below methods.

1. Taylor Series 2. Composite Trapezoidal rule 3. Composite Simpson's rule 4. Romberg table (including Richardson's extrapolation)

Solution
 Taylor Series : 

Pros

1. Powerful tool that helps to approximate any type of simple or complex functions into a polynomial form with an error in the form of remainder.

2. 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.

Cons

1. The function should be differentiable in the given limits in which it is approximated as a polynomial.

2. Number of operations can be large. This is not always a problem these days considering how cheap and fast computers are at  present. In general the step size can be made larger the higher the order of the Taylor scheme. But a computational count will tell you whether it is more effective to compute a lot at each step, and take larger steps, or compute little at each step and make due with smaller step sizes. Usually it is more advantageous to go with good lower order method and small steps, but this depends on the problem.

Composite Trapezoidal rule:

Pros:

1. This methods works well for periodic functions like trigonometric functions (sinx, cosx) and it is very accurate because the error is dependent on the difference between the odd derivatives of the limits.

2. Better accuracy then simple trapezoidal rule and also easy to use.

Cons :

1. It's convergence rate is low and we need to have more number of sub intervals.

Composite Simpson's rule:

Pros:

1. Exact for cubic (less then that) polynomials.

2. More accurate then composite trapezoidal rule.

Cons:

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

 Romberg table (including Richardson's extrapolation): 

Pros:

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

2. Better accuracy then other methods.

Cons:

1. There is one limitation with this rule.Romberg integration is successful only when integrand satisfies the hypotheses of the Euler–Maclaurin.

Abhishekksingh 07:55, 3 March 2010 (UTC)

= Proof: Higher order error Trapezoidal rule : Integration by parts = Ref: Lecture Notes [[media:Egm6341.s10.mtg21.djvu|p.21-2]]

Problem Statement
In the proof of the Higher order error of the Trapezoidal rule, show that,
 * {| style="width:100%" border="0" align="left"

\int_{-1}^{1}(-t)g^{(1)}(t)\,dt = \int_{-1}^{1}(-t)g(t)\,dt - [g(-1) + g(1)] $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$\displaystyle (Eq. 3)$$
 * }

Solution
Using integration by parts, we know that,
 * {| style="width:100%" border="0" align="left"

\int_{a}^{b}u\,dv = [u v]_{a}^{b} -\int_{a}^{b}v\,du $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$\displaystyle (Eq. 4)$$
 * }

Comparing the LHS of Eqn(3) and Eqn (4)
 * {| style="width:100%" border="0" align="left"

u = -t \,\,\,\, and\,\,\, \, dv = g^{(1)}(t)\,dt $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$\displaystyle (Eq. 5a)$$


 * }


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

du = \,-dt \,\,\, and\,\,\,\, v = g(t) \, $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$\displaystyle (Eq. 5b)$$
 * }

Applying Eqn(4),Eqn(5a) and Eqn(5b) in Eqn(3),
 * {| style="width:100%" border="0" align="left"

\int_{-1}^{1}(-t)g^{(1)}(t)\,dt = [(-t)(g(t))]_{-1}^{1} + \int_{-1}^{1}-g(t)\,dt $$
 * $$ \displaystyle
 * $$ \displaystyle
 * }


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

= - [(g(1)) - (g(-1))] + \int_{-1}^{1}g(t)\,dt $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$ \displaystyle
 * }

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

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

Subramanian Annamalai 18:21, 1 March 2010 (UTC)

= Proof: Higher order error Trapezoidal rule : Chain Rule = Ref: Lecture Notes [[media:Egm6341.s10.mtg21.djvu|p.21-2]]

Problem Statement
In the proof of the Higher order error of the Trapezoidal rule, show that,
 * {| style="width:100%" border="0" align="left"

g_k^{(i)} = \bigg(\frac{h}{2}\bigg) ^{(i)} \cdot f^{(i)}(x(t)) \,\,\,\,\,\,\,\, \in [x_k, x_{k+1}] $$
 * $$ \displaystyle
 * $$ \displaystyle
 * }

Solution
Recall that,
 * {| style="width:100%" border="0" align="left"

g_k(t) = f(x(t)) $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$\displaystyle (Eq. 6)$$
 * }

and


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

x(t):= t\cdot\bigg(\frac{h}{2}\bigg) + \frac{x_k + x_{k+1}}{2} \,\,\,\,\,\,\,\, t \in [-1 ,1] $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$\displaystyle (Eq. 7)$$
 * }

Applying Chain rule,
 * {| style="width:100%" border="0" align="left"

g_k^{(1)} = \frac{dg}{dt} = \frac{dg}{dx}\cdot\frac{dx}{dt} $$
 * $$ \displaystyle
 * $$ \displaystyle


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

= f^{(1)}(x(t))\cdot\frac{d}{dt}\big(t\cdot\bigg(\frac{h}{2}\bigg) + \frac{x_k + x_{k+1}}{2}\big) $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$ \displaystyle
 * }


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

= f^{(1)}(x(t))\cdot \frac{h}{2} $$
 * $$ \displaystyle
 * $$ \displaystyle
 * $$ \displaystyle
 * }

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

g_k^{(2)} = \frac{d}{dt}\bigg(\frac{dg}{dt}\bigg) $$
 * $$ \displaystyle
 * $$ \displaystyle


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

= \frac{d}{dt}\Bigg[(f^{(1)}(x(t))\cdot \bigg(\frac{h}{2}\bigg)\Bigg] $$
 * $$ \displaystyle
 * $$ \displaystyle
 * }


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

=\bigg(\frac{h}{2}\bigg)\cdot\frac{d}{dx}\frac{dx}{dt} $$
 * $$ \displaystyle
 * $$ \displaystyle
 * }


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

= f^{(2)}(x(t))\cdot \bigg(\frac{h}{2}\bigg)^2 $$
 * $$ \displaystyle
 * $$ \displaystyle
 * }


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

.

.

.

.


 * }

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

$$ \displaystyle g_k^{(i)} = \bigg(\frac{h}{2}\bigg) ^{(i)} \cdot f^{(i)}(x(t)) $$
 * style="width:25%; padding:10px; border:2px solid #8888aa" |
 * style="width:25%; padding:10px; border:2px solid #8888aa" |
 * style = |
 * }

Subramanian Annamalai 18:21, 1 March 2010 (UTC)

= Proof: Higher order Trapezoidal rule : Integration by parts =

Ref. Lecture notes [[media:Egm6341.s10.mtg20.djvu|p.21-2]]

Problem Statement
Prove that $$E:= \int\limits_{1}^{-1}(-t)g^{(1)}(t) dt = [P_2(t)g^{(1)}(t)]_{-1}^{1} - \int\limits_{1}^{-1} P_2(t) g^{(2)}(t) dt   $$

where

$$ P_1(t)= -t $$

and

$$ P_2(t)= \int\limits_{1}^{-1} P_1(t) dt $$

Solution
$$ E:= \int\limits_{1}^{-1} \underbrace{(-t)}_{P_1(t)} g^{(1)}(t) dt $$

$$  = \int\limits_{1}^{-1} P_1(t) g^{(1)}(t) dt  $$

$$  = g^{(1)}(t) \int P_1(t) dt - \int [\frac{d}{dt}g^{(1)}(t) \int P_1(t) dt]dt $$

$$  = [g^{(1)}(t) P_2(t))]_{-1}^{1} - \int [g^{(2)}(t) P_2(t)]dt $$

$$  where, P_2(t) = \int P_1(t) dt                        = \int (-t) dt                        = \frac{-t^2}{2} + {\alpha}\rightarrow_{Integration Constt.} $$

Abhishekksingh 07:55, 3 March 2010 (UTC)

= Comparison of various numerical integration codes for the function $$ f(x) = (e^x-1)/x$$ =

Install Chebfun and practice sum commond
Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-1]]

Problem Statement
Install chebfun and practice sum command to do integration$$\displaystyle I = \int_{0}^{1}\frac{e^{x}-1}{x}dx$$, and record the computational time

Solution
i) Chebfun installation The chebfun can be downloaded from http://www2.maths.ox.ac.uk/chebfun/files.html

We installed the "chebfun" successful based on the instruction which is also availabe in the above webpage.

In chebfun, 'sum' command will returns the definite integral of a chebfun over its range of defination.

The integral is normally calculated by an FFT-based variant of Clenshaw-Curtis quadrature.

This formula is applied on each fun and then the results are added up.

The chebfun will be used in the following problem to show comparision of speed among different numerical methods for integration.

For detailed information, one can refer the paper: Battles and Trefethen, An extension of Matlab to continuous functions and operators (SIAM J. Sci. Comp., 2004)

ii) Chebfun, Sum command

We use the command 'sum' in chebfun to compute I with error to $$ 10^{-10}$$

Result:

Integral value =1.317902151454404

Error =4.403588604873221e-012

Elapsed time is 0.139211 seconds.

-- By Min Zhong 12:00, 02 March 2010 (UTC)

Problem Statement
$$ \displaystyle Compare\, the\, second\, column\, of\, Romberg\, Table\, Romberg\, to\, Simpson's\, rule\, and \, derive any relationship,\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-1]]

Solution
$$ \displaystyle a)\, Formula\, for\, Composite\, Trap.\, rule\, is\, given\, by\, $$

$$ \int_{a}^{b}f(x)dx=\frac{b-a}{2n}*[f(x_{0})+2f(x_{1})+2f(x_{2})+.........+2f(x_{n-1})+f(x_{n})] $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg7.pdf|p.7-1]]

$$ \displaystyle b)\, Formula\, for.\, Romberg\, table\, is\, given\, by\, $$

$$ T_{k}(n)=\frac{2^{2k}T_{k-1}(n)-T_{k-1}(\frac{n}{2})}{2^{2k}-1} $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg19.djvu|p.19-2]]

$$ \displaystyle Since\, the\, 2nd\, column\, of\, Romberg\, table\, corresponds\, to\, k=1\, $$

$$ T_{1}(n)=\frac{4T_{0}(n)-T_{0}(\frac{n}{2})}{3} $$

$$ \displaystyle T_{0}(n)\, and\, T_{0}(\frac{n}{2})\, are\, composite\, trap.\, rule\, $$

$$ \displaystyle i)\, T_{0}(\frac{n}{2})= \frac{b-a}{\frac{n}{2}}[\frac{1}{2}f_{0}+f_{2}+...+f_{n-2}+\frac{1}{2}f_{n}]\, =\, \frac{b-a}{n}[f_{0}+2f_{2}+...+2f_{n-2}+2f_{n-1}+f_{n}] $$

$$ \displaystyle Where\, the\, number\, of\, pannels\, =\, \frac{n}{2}\,, \,\, h=\frac{b-a}{\frac{n}{2}}\, $$

$$ \displaystyle ii)\, T_{0}(n)= \frac{b-a}{n}[\frac{1}{2}f_{0}+f_{1}+f_{2}+...+f_{n-2}+f_{n-1}+\frac{1}{2}f_{n}] $$

$$ \displaystyle Where\, the\, number\, of\, pannels\, =\, n\,, \,\, h=\frac{b-a}{n}\, $$

$$ \displaystyle iii)\, T_{1}(n)= \frac{1}{3}\underbrace{[4T_{0}(n)-T_{0}(\frac{n}{2})]}_{call\, p} $$

$$ \displaystyle Where\, p\, =\, 4T_{0}(n)-T_{0}(\frac{n}{2})\, =\, \frac{b-a}{n}[f_{0}+4f_{1}+2f_{2}+...+2f_{n-2}+4f_{n-1}+f_{n}] $$

$$ \begin{align} \displaystyle Therefore,\, T_{1}(n)&=\frac{1}{3}p \\ &=\frac{1}{3}(\frac{b-a}{n}[f_{0}+4f_{1}+2f_{2}+...+2f_{n-2}+4f_{n-1}+f_{n}]) \\ &=\underbrace{\frac{h}{3}[f_{0}+4f_{1}+2f_{2}+...+2f_{n-2}+4f_{n-1}+f_{n}])}_{Exactly same with the formula of Comp. Simp. rule} \\\end{align}\ $$

--Heejun Chung 18:31, 5 March 2010 (UTC)

Integration:Using composite Trap rule, composite Simpson's rule and Romberg Table
Ref: Lecture notes [[media:Egm6341.s10.mtg25.djvu|p.25-1]]

Problem Statement
For $$\displaystyle I = \int_{0}^{1}\frac{e^{x}-1}{x}dx$$, compare $$\displaystyle I_{n}$$ to $$\displaystyle O(10^{-10})$$ and computational time for - Composite Trap rule, - Composite Simpson's rule and - Romberg table.

Solution
In what follows, we are using the follow matlab function as a given function $$\displaystyle f(x) = \frac{e^{x}-1}{x}$$.

i) Composit Trap. rule

Recall the formula of comp. Trap. rule.


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

T_{0}(n) = h\left[\frac{1}{2}f(x_{0}) + f(x_{1}) + f(x_{2}) + \cdot\cdot\cdot + f(x_{n-1}) + \frac{1}{2}f(x_{n}) \right] $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, realize the fact that if we want to double the number of sub-intervals, we don't have to calculate the functions at even point, i.e. $$\displaystyle f(x_{0}), f(x_{2}), \cdot\cdot\cdot, f(x_{n-2}), f(x_{n}) $$ because we already calculated the functions at these even point at the previous step. Therefore, all we have to do is evaluating the functions only at the odd points and combine them with the calculated values at the previous step. As a result, we can use the following inductive formula.
 * {| style="width:100%" border="0" align="left"

T_{0}(2n) = T_{0}(n) + h\sum_{i=1}^{n}f(x_{2i-1}). $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,trap} = T_{0}(n) = 1.317902151493210, \quad n = 2^{15} = 32768 $$ $$\displaystyle E_{n,trap} = I - I_{n,trap} = 3.880606946893295e-11 $$

ii) Composite Simpson's rule

Recall the formula of comp. Simpson's rule.


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

T_{1}(n) = \frac{h}{3}\left[f(x_{0}) + 4f(x_{1}) + 2f(x_{2}) + 4f(x_{3}) + \cdot\cdot\cdot + 4f(x_{n-3}) + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_{n}) \right] $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, realize the fact that if we want to double the number of sub-intervals, we don't have to calculate the functions at even point, i.e. $$\displaystyle f(x_{0}), f(x_{2}), \cdot\cdot\cdot, f(x_{n-2}), f(x_{n}) $$ because we already calculated the functions at these even point at the previous step. Therefore, all we have to do is evaluating the functions only at the odd points and combine them with the calculated values at the previous step. As a result, we can use the following inductive formula.
 * {| style="width:100%" border="0" align="left"

T_{1}(2n) = T_{1}(n) + h\sum_{i=1}^{n}4f(x_{2i-1}) - \underbrace{h\sum_{i=1}^{\frac{n}{4}}2f(x_{4i-2})}_{Already\,\, obtained\,\, when\,\, T_{1}(n) \,\, was\,\, calculated.}. $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,Simpson} = T_{1}(n) = 1.317902151460891, \quad n = 2^{7} = 128 $$ $$\displaystyle E_{n,Simpson} = I - I_{n,Simpson} = 6.486811088279865e-12 $$

iii) Romberg Table

Since the first column of Romberg table corresponds to Trap. rule, we can use the below formulas in order to obtain the first column of Romberg talbe.


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

T_{0}(n) = h\left[\frac{1}{2}f(x_{0}) + f(x_{1}) + f(x_{2}) + \cdot\cdot\cdot + f(x_{n-1}) + \frac{1}{2}f(x_{n}) \right] $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }


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

T_{0}(2n) = T_{0}(n) + h\sum_{i=1}^{n}f(x_{2i-1}). $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, we gradually set up the Romberg table in parallel, as soon as the result for $$\displaystyle T_0(2^j) $$ becomes available using the below formula.


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

T_{k}(2n) = \frac{2^{2k}T_{k-1}(2n) - T_{k-1}(n)}{2^{2k}-1}. $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,Romberg} = T_{3}(n) = 1.317902151489847, \quad n = 2^{3} = 8 $$ $$\displaystyle E_{n,Romberg} = I - I_{n,Romberg} = 3.544298188273842e-11 $$

iv) Compare computational time

We can compare computational times of above three methods using tic/toc commands in Matlab.

Matlab code

As above table shown, Romberg Table has the shortest computational time. Note that Simpson's rule also has the same order of magnitude of computational time with Romberg Table even though Romberg Table has slightly shorter computational time than Simpson's rule. Therefore,  both Romberg table and Simpson's rule are efficient ways to get a numerical integration value for the given function.

Yong Nam Ahn 22:49, 1 March 2010 (UTC)

Compare trapz, quad, clencurt MATLAB functions
Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Problem Statement
For $$\displaystyle I = \int_{0}^{1}\frac{e^{x}-1}{x}dx$$, compare $$\displaystyle I_{n}$$ to $$\displaystyle O(10^{-10})$$ and computational time by - Matlab : trapz, quad and clencurt

Solution
i) Trapz

We use the command 'trapz' in matlab to compute I with error to $$ 10^{-10}$$

Result

Integral value = 1.317902151493208

n = 32768

Error = 3.880429311209355e-011

Elapsed time is 0.753827 seconds

ii) Quad

We use the command 'quad' in matlab to compute I with error to $$ 10^{-10}$$

Result:

Integral value =1.317902151454412

Elapsed time is 0.043016 seconds.

iii) Clencurt

We use the command 'clencurt' in matlab to compute I with error to $$ 10^{-10}$$

Result:

Integral value = 1.317902151447634

Elapsed time is 0.001750 seconds.

iv) Compare computational time

As above table shown, clencurt has the shortest computational time and trapz takes longest time.

Therefore,clencurt is the more efficient way to get a numerical integration value for the given function.

-- By Min Zhong 12:00, 02 March 2010 (UTC)

Problem Statement
$$ \displaystyle Compare\, the\, time\, for\, above\, each\, method\, with\, comment\, on\, the\, results\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-1 & 25-2]]

Solution
''' The Romberg Table, Comp. Simpson's rule and Clencurt are efficient ways to get a numerical integration value for ''' $$\displaystyle I = \int_{0}^{1}\frac{e^{x}-1}{x}dx$$

--Heejun Chung 18:36, 5 March 2010 (UTC)

= Comparison of various numerical integration codes for the function $$f(x) = 1/(1+x^{2})$$ =

Apply chebfun to do integration
Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Problem Statement
For $$\displaystyle I = \int_{-5}^{5}\frac{1}{1+x^2}dx$$, compare $$\displaystyle I_{n}$$ to $$\displaystyle O(10^{-10})$$ and computational time by - chebfun: sum command,

Solution
 Chebfun

We use the command 'sum' in chebfun to compute I with error to $$ 10^{-10}$$

Result:

Integral value = 2.746801533890032 Elapsed time is 0.141973 seconds. Error=4.440892098500626e-016

Integration: Using composite Trap rule, composite Simpson's rule and Romberg Table
Ref: Lecture notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Problem Statement
For $$\displaystyle I = \int_{-5}^{5}\frac{1}{1+x^2}dx$$, compare $$\displaystyle I_{n}$$ to $$\displaystyle O(10^{-10})$$ and computational time for - Composite Trap rule, - Composite Simpson's rule and - Romberg table.

Solution
In what follows, we are using the follow matlab function as a given function $$\displaystyle f(x) = \frac{1}{1+x^{2}}$$.

i) Composit Trap. rule

Recall the formula of comp. Trap. rule.


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

T_{0}(n) = h\left[\frac{1}{2}f(x_{0}) + f(x_{1}) + f(x_{2}) + \cdot\cdot\cdot + f(x_{n-1}) + \frac{1}{2}f(x_{n}) \right] $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, realize the fact that if we want to double the number of sub-intervals, we don't have to calculate the functions at even point, i.e. $$\displaystyle f(x_{0}), f(x_{2}), \cdot\cdot\cdot, f(x_{n-2}), f(x_{n}) $$ because we already calculated the functions at these even point at the previous step. Therefore, all we have to do is evaluating the functions only at the odd points and combine them with the calculated values at the previous step. As a result, we can use the following inductive formula.
 * {| style="width:100%" border="0" align="left"

T_{0}(2n) = T_{0}(n) + h\sum_{i=1}^{n}f(x_{2i-1}). $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,trap} = T_{0}(n) = 2.746801533832637, \quad n = 2^{16} = 65536 $$ $$\displaystyle E_{n,trap} = I - I_{n,trap} = 5.739320130260239e-11 $$

ii) Composite Simpson's rule

Recall the formula of comp. Simpson's rule.


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

T_{1}(n) = \frac{h}{3}\left[f(x_{0}) + 4f(x_{1}) + 2f(x_{2}) + 4f(x_{3}) + \cdot\cdot\cdot + 4f(x_{n-3}) + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_{n}) \right] $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, realize the fact that if we want to double the number of sub-intervals, we don't have to calculate the functions at even point, i.e. $$\displaystyle f(x_{0}), f(x_{2}), \cdot\cdot\cdot, f(x_{n-2}), f(x_{n}) $$ because we already calculated the functions at these even point at the previous step. Therefore, all we have to do is evaluating the functions only at the odd points and combine them with the calculated values at the previous step. As a result, we can use the following inductive formula.
 * {| style="width:100%" border="0" align="left"

T_{1}(2n) = T_{1}(n) + h\sum_{i=1}^{n}4f(x_{2i-1}) - \underbrace{h\sum_{i=1}^{\frac{n}{4}}2f(x_{4i-2})}_{Already\,\, obtained\,\, when\,\, T_{1}(n) \,\, was\,\, calculated.}. $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,Simpson} = T_{1}(n) = 2.746801533879843, \quad n = 2^{9} = 512 $$ $$\displaystyle E_{n,Simpson} = I - I_{n,Simpson} = 1.018696238475059e-11 $$

iii) Romberg Table

Since the first column of Romberg table corresponds to Trap. rule, we can use the below formulas in order to obtain the first column of Romberg talbe.


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

T_{0}(n) = h\left[\frac{1}{2}f(x_{0}) + f(x_{1}) + f(x_{2}) + \cdot\cdot\cdot + f(x_{n-1}) + \frac{1}{2}f(x_{n}) \right] $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }


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

T_{0}(2n) = T_{0}(n) + h\sum_{i=1}^{n}f(x_{2i-1}). $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, we gradually set up the Romberg table in parallel, as soon as the result for $$\displaystyle T_0(2^j) $$ becomes available using the below formula.


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

T_{k}(2n) = \frac{2^{2k}T_{k-1}(2n) - T_{k-1}(n)}{2^{2k}-1}. $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,Romberg} = T_{2}(n) = 2.746801533889937, \quad n = 2^{8} = 256 $$ $$\displaystyle E_{n,Romberg} = I - I_{n,Romberg} = 9.281464485866309e-14 $$

iv) Compare computational time

We can compare computational times of above three methods using tic/toc commands in Matlab.

Matlab code

As above table shown, Simpson's rule has the shortest computational time. However, note that Romberg Table has almost same computational time with Simpson's rule even though Simpson's rule has slightly shorter computational time than Romberg Table. Therefore,  both Simpson's rule and Romberg Table are efficient ways to get a numerical integration value for the given function.

Yong Nam Ahn 22:49, 1 March 2010 (UTC)

Compare trapz, quad, clencurt functions of MATLAB
Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Problem Statement
For $$\displaystyle I = \int_{-5}^{5}\frac{1}{1+x^2}dx$$, compare $$\displaystyle I_{n}$$ to $$\displaystyle O(10^{-10})$$ and computational time for - chebfun: sum command, - Matlab : trapz, quad and clencurt

Solution
i) Trapz

We use the command 'trapz' in matlab to compute I with error to $$ 10^{-10}$$

Result:

Integral value = 2.746801533832626

n = 65536

Error = 5.740607988968804e-011

Elapsed time is 0.061356 seconds.

ii) Quad

We use the command 'quad' in matlab to compute I with error to $$ 10^{-10}$$

Result:

I_q= 2.746801533869947

Elapsed time is 0.029688 seconds.

iii) Clencurt

We use the command 'clencurt' in matlab to compute I with error to $$ 10^{-10}$$

Result:

Integral value = 2.746801533874087

n = 62

Elapsed time is 0.011186 seconds.

iv) Compare computational time

As above table shown, Quad has the shortest computational time among these commonds.

Therefore, Quad is a more efficient way to get a numerical integration value for the given function.

-- By Min Zhong 2:00, 03 March 2010 (UTC)

Problem Statement
$$ \displaystyle Compare\, the\, time\, for\, above\, each\, method\, with\, comment\, on\, the\, results\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Solution
 Both Simpson's rule and Romberg Table are efficient ways to get a numerical integration value for  $$\displaystyle I = \int_{-5}^{5}\frac{1}{1+x^2}dx$$

--Heejun Chung 18:43, 5 March 2010 (UTC)

= Comparison of various numerical integration codes for the given function $$\ f(\theta) = \frac{1-sin^{2}(\frac{\pi}{12})}{1-sin(\frac{\pi}{12})cos\theta}$$ =

Apply chebfun for integration
Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Problem Statement
Compute $$\displaystyle I = \int_{0}^{2\pi}\frac{1-sin^{2}(\frac{\pi}{12})}{1-sin(\frac{\pi}{12})cos\theta}d\theta$$, and computational time - chebfun: sum command,

Solution
 Chebfun

We use the command 'sum' in chebfun to compute I with error to $$ 10^{-10}$$

Result:

Integral value= 6.069090959564774

Error =3.552713678800501e-015

Elapsed time is 0.140592 seconds.

Integration: Using composite Trapezoidal rule, composite Simpson's rule and Romberg Table
Ref: Lecture notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Problem Statement
For $$\displaystyle I = \int_{0}^{2\pi}\frac{1-sin^{2}(\frac{\pi}{12})}{1-sin(\frac{\pi}{12})cos\theta}d\theta$$, compare $$\displaystyle I_{n}$$ to $$\displaystyle O(10^{-10})$$ and computational time for - Composite Trap rule, - Composite Simpson's rule and - Romberg table.

Solution
In what follows, we are using the follow matlab function as a given function $$\displaystyle f(\theta) = \frac{1-sin^{2}(\frac{\pi}{12})}{1-sin(\frac{\pi}{12})cos\theta}$$.

i) Composit Trap. rule

Recall the formula of comp. Trap. rule.


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

T_{0}(n) = h\left[\frac{1}{2}f(x_{0}) + f(x_{1}) + f(x_{2}) + \cdot\cdot\cdot + f(x_{n-1}) + \frac{1}{2}f(x_{n}) \right] $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, realize the fact that if we want to double the number of sub-intervals, we don't have to calculate the functions at even point, i.e. $$\displaystyle f(x_{0}), f(x_{2}), \cdot\cdot\cdot, f(x_{n-2}), f(x_{n}) $$ because we already calculated the functions at these even point at the previous step. Therefore, all we have to do is evaluating the functions only at the odd points and combine them with the calculated values at the previous step. As a result, we can use the following inductive formula.
 * {| style="width:100%" border="0" align="left"

T_{0}(2n) = T_{0}(n) + h\sum_{i=1}^{n}f(x_{2i-1}). $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,trap} = T_{0}(n) = 6.069090959564874, \quad n = 2^{4} = 16 $$ $$\displaystyle E_{n,trap} = I - I_{n,trap} = 1.039168751049147e-13 $$

ii) Composite Simpson's rule

Recall the formula of comp. Simpson's rule.


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

T_{1}(n) = \frac{h}{3}\left[f(x_{0}) + 4f(x_{1}) + 2f(x_{2}) + 4f(x_{3}) + \cdot\cdot\cdot + 4f(x_{n-3}) + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_{n}) \right] $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, realize the fact that if we want to double the number of sub-intervals, we don't have to calculate the functions at even point, i.e. $$\displaystyle f(x_{0}), f(x_{2}), \cdot\cdot\cdot, f(x_{n-2}), f(x_{n}) $$ because we already calculated the functions at these even point at the previous step. Therefore, all we have to do is evaluating the functions only at the odd points and combine them with the calculated values at the previous step. As a result, we can use the following inductive formula.
 * {| style="width:100%" border="0" align="left"

T_{1}(2n) = T_{1}(n) + h\sum_{i=1}^{n}4f(x_{2i-1}) - \underbrace{h\sum_{i=1}^{\frac{n}{4}}2f(x_{4i-2})}_{Already\,\, obtained\,\, when\,\, T_{1}(n) \,\, was\,\, calculated.}. $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,Simpson} = T_{1}(n) = 6.069090959564741, \quad n = 2^{5} = 32 $$ $$\displaystyle E_{n,Simpson} = I - I_{n,Simpson} = 2.930988785010413e-14 $$

iii) Romberg Table

Since the first column of Romberg table corresponds to Trap. rule, we can use the below formulas in order to obtain the first column of Romberg talbe.


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

T_{0}(n) = h\left[\frac{1}{2}f(x_{0}) + f(x_{1}) + f(x_{2}) + \cdot\cdot\cdot + f(x_{n-1}) + \frac{1}{2}f(x_{n}) \right] $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }


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

T_{0}(2n) = T_{0}(n) + h\sum_{i=1}^{n}f(x_{2i-1}). $$ where $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle x_{0} = a,\quad x_{n} = b $$.
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Also, we gradually set up the Romberg table in parallel, as soon as the result for $$\displaystyle T_0(2^j) $$ becomes available using the below formula.


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

T_{k}(2n) = \frac{2^{2k}T_{k-1}(2n) - T_{k-1}(n)}{2^{2k}-1}. $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Matlab code

Result

$$\displaystyle I_{n,Romberg} = T_{0}(n) = 6.069090959564874, \quad n = 2^{4} = 16 $$ $$\displaystyle E_{n,Romberg} = I - I_{n,Romberg} = 1.039168751049147e-13 $$

iv) Compare computational time

We can compare computational times of above three methods using tic/toc commands in Matlab.

Matlab code

As above table shown, Simpson's rule has the shortest computational time. However, we should note that the computational time does not really different from each other for three method. Especially, trap rule has a shorter computational time than Romberg Table (in general, computational time of Romberg Table is much shorter than that of Trap. rule.). Since the given function $$\displaystyle f(\theta) = \frac{1-sin^{2}(\frac{\pi}{12})}{1-sin(\frac{\pi}{12})cos\theta}$$ is a periodic function, the numerical integration of this function is converged very fast by Trap. rule. Therefore, '''Trap. rule is an efficient way to numerically integrate the given function $$\displaystyle f(\theta) = \frac{1-sin^{2}(\frac{\pi}{12})}{1-sin(\frac{\pi}{12})cos\theta}$$ because this function is a periodic function.'''

Yong Nam Ahn 22:49, 1 March 2010 (UTC)

Compare trapz, quad, clencurt MATLAB functions
Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Problem Statement
For $$\displaystyle I = \int_{0}^{2\pi}\frac{1-sin^{2}(\frac{\pi}{12})}{1-sin(\frac{\pi}{12})cos\theta}d\theta$$,

- Matlab : trapz, quad and clencurt

Solution
i) Trapz

We use the command 'trapz' in matlab to compute I with error to $$ 10^{-10}$$

Result:

Integral value =6.069090959564874

n =16

Error =1.039168751049147e-013

Elapsed time is 0.0536801 seconds.

ii) Quad

We use the command 'quad' in matlab to compute I with error to $$ 10^{-10}$$

Result:

Integral value = 6.069090959562212

Elapsed time is 0.023456 seconds.

iii) Clencurt

We use the command 'clencurt' in matlab to compute I with error to $$ 10^{-10}$$

Result:

Integral value = 6.069090959664361

n = 20

Elapsed time is 0.006226 seconds.

iv) Compare computational time

As above table shown, clencurt has the shortest computational time among these commonds .

Therefore, cleancurt is a more efficient way to get a numerical integration value for the given function.

-- By Min Zhong 12:00, 02 March 2010 (UTC)

Problem Statement
$$ \displaystyle Compare\, the\, time\, for\, above\, each\, method\, and\, comment\, on\, the\, results\, $$

Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|p.25-2]]

Solution
''' Comp. Trap. rule is an efficient way to numerically integrate the given function $$\displaystyle f(\theta) = \frac{1-sin^{2}(\frac{\pi}{12})}{1-sin(\frac{\pi}{12})cos\theta}$$ because this function is a periodic function.'''

--Heejun Chung 18:43, 5 March 2010 (UTC)

= Contributing Team Members =
 * 1) Subramanian Annamalai 18:21, 1 March 2010 (UTC) Authored: 1,3 and 5   Proof-read:6(c) 7(c) and 8(c)
 * 2) Yong Nam Ahn 22:49, 1 March 2010 (UTC) Authored: 6(c), 7(c) and 8(c),   Proof-read: 6(b), 6(e), 7(e) and 8(e)
 * 3) Abhishekksingh 08:02, 3 March 2010 (UTC) Authored: 2 and 5 , Proof-read : 6(a) , 6(d) ,7(a),7(d),8(a) and 8(d)
 * 4) Min Zhong 11:21, 3 March 2010 (UTC) Authored: 6(a), 6(d, 7(d) and 8(d),  Proof-read: 2 and 5
 * 5) Heejun Chung 18:43, 5 March 2010 (UTC) Authored: 6(b), 6(e), 7(e) and 8(e),  Proof-read: 1, 3 and 4

=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