User:Egm6341.s10.team3.ynahn81/HW4-7(c)

(c) Compare Comp Trap, Comp Simp, Romberg for $$\displaystyle f(x) = \frac{1}{1+x^{2}}$$
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.