User:Egm6341.s2010.Team1/HW4

=Problem 1-Develop higher-order corrected trap.rule =

Find
Find $$ I_{n}\,$$ using$$ CT_{1},CT_{2},CT_{3}\,$$ for $$  n=2,4,8,16,...\,$$ untill  $$ I-I_{n}=O(10^{-6}) \,$$

Given
[[media:Egm6341.s10.mtg20.djvu|p.20-1]]

$$ CT_{K}(n)=CT_{K-1}(n)+a_{k}h^{2*k} \,$$

see [[media:Egm6341.s10.mtg6.pdf|p.6-5]] and [[media:Egm6341.s10.mtg19.djvu|p.19-2]]

Solution
From [[media:Egm6341.s10.mtg19.djvu|p.19-3]] and [[media:Egm6341.s10.mtg20.djvu|p.20-1]]

Case 1,find $$ I_{n}\, $$ using $$ CT_{1}\,$$

$$ CT_{1}(n)=\underbrace{\left( CT_{0}(n)\right)}_{T_{0}(n)} +a_{1}h^{2*1}\,$$

$$ I= \underbrace{\left( CT_{1}(n)\right)}_{I_{n}} +O(h^{4})\,$$

$$ CT_{2}(n)=\underbrace{\left( CT_{1}(n)\right)}_{T_{1}(n)} +a_{2}h^{2*2}\,$$

$$ I= \underbrace{\left( CT_{2}(n)\right)}_{I_{n}} +O(h^{6})\,$$

$$ CT_{3}(n)=\underbrace{\left( CT_{2}(n)\right)}_{T_{2}(n)} +a_{3}h^{2*6}\,$$

$$ I= \underbrace{\left( CT_{3}(n)\right)}_{I_{n}} +O(h^{8})\,$$

where $$ a_{1},a_{2},a_{3} \,$$ can be derived from [[media:Egm6341.s10.mtg19.djvu|p.19-3]]

$$ a_{i}=d_{i}[f^{2i-1}(b)-f^{2i-1}(a)]\,$$

$$ d_{i}=- \frac{B_{2i}}{(2i)!}\, $$

$$ d_{1}=- \frac{1}{12}, d_{2}= \frac{1}{720}, d_{3}= - \frac{1}{30240}\,$$

The above integral can be caculated using MATLAB codes and corresponding errors are tabulated.

From the results above, it can be easily seen that for $$ CT_{1}\,$$ the error can reach $$ O(10^{-6})\,$$ when n=4; for $$ CT_{2}\,$$ the error can already reach $$ O(10^{-7}\,$$) when n=2; for $$ CT_{3}\,$$ the error can already reach $$ O(10^{-9})\,$$ when n=2. So corrected trap.rule is more accurate than previous methods.

Matlab Code:

For $$ CT_{1},$$ and corresponding error

Matlab Code:

For $$ CT_{2},$$ and corresponding error Matlab Code:

For $$ CT_{3},$$ and corresponding error

Egm6341.s10.team1.lei16:31, 3 March 2010 (UTC)

=Problem 2 - Find which n,st error equal 0 =

Find
For which $$ n,\,$$ so that $$ E_{n}^{1}=0 \, $$

Given
[[media:Egm6341.s10.mtg20.djvu|p.20-2]] equation(1)

$$ f \,$$ is  periodic  on$$ [a,b] \Longleftrightarrow f^{(k)}(a)=f^{(k)}(b),   k=0,1,2,...\,$$ Clearly (1) is valid for odd numb.

Solution
From [[media:Egm6341.s10.mtg19.djvu|p.19-3]]

$$ a_{i}=d_{i}[f^{2i-1}(b)-f^{2i-1}(a)]\,$$

if (1) is valid, then

$$ a_{i}=0 \,$$

Recall [[media:Egm6341.s10.mtg18.djvu|p.18-2]] and [[media:Egm6341.s10.mtg18.djvu|p.18-3]]

$$ E_{n}^{1}=a_{1} h^{2*1}+a_{2} h^{2*2}+... =\sum _{i=1}^{\infty}a_{i}h^{2*i}\,$$

where $$ a_{i}=0\,$$, therefore no matter $$ h\,$$,$$ E_{n}^{1}=0\,$$

That means $$ E_{n}^{1}\,$$ does not depend on $$ n\,$$,if the $$ f \,$$ is  periodic  on$$ [a,b]\,$$

Egm6341.s10.team1.lei16:32, 3 March 2010 (UTC)

=Problem 3 -Discuss Pros and cons =

Find
Discuss Pros and cons:

1) $$ Taylor Series\,$$

2) $$ Compos. Trap.\,$$

3) $$ Compos. Simpson\,$$

4) $$ Romberg Table (Richardson)\,$$

5) $$ CT_{K}(n)\,$$

Solution
1) $$ Taylor  Series\,$$

$$ \bullet\,$$ Accuracy and speediness. $$ \bullet\,$$ Number of operations can be large. $$ \bullet\,$$ Need to know about smoothness of solution.It is senseless to use high order method when the solution has unbounded high order derivatives. $$ \bullet\,$$ 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.

2) $$ Compos.  Trap.\,$$

$$ \bullet\,$$ Approximate derivative by constant. $$ \bullet\,$$ Average (not discriminate) endpoints. $$ \bullet\,$$ Often converges very quickly for periodic functions, and tends to become extremely accurate when periodic functions are integrated over their periods. $$ \bullet\,$$ For various classes of functions that are not twice-differentiable, the trapezoidal rule has sharper bounds than Simpson's rule.

3) $$ Compos.  Simpson\,$$

$$ \bullet\,$$ In general has faster convergence than the trapezoidal rule However for various classes of rougher functions (ones with weaker smoothness conditions), the trapezoidal rule has faster convergence in general than Simpson's rule. $$ \bullet\,$$ Provides exact results for any polynomial of degree three or less.

4)$$ Romberg  Table\,$$

$$ \bullet\,$$ By using Richardson extrapolation repeatedly on the trap. rule is more efficient. $$ \bullet\,$$ Converge to high accuracy.

5) $$ CT_{K}(n)\,$$

$$ \bullet\,$$ More accurate than composite trap.rule and Romberg table. $$ \bullet\,$$ Converge more quickly than trap.rule.

=Problem 4-Verify change of variables=

Given
Higher Order Error for Trap. rule equation:


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle \ 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}))]
 * }
 * }

Find
Show that using change of variables this equation becomes


 * {| 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


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

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

And


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

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

Solution
Using the function
 * {| style="width:100%" border="0" align="left"

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

This can be changed into


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

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

This allows us to evaluate the function $$x(t)$$ at the endpoints to obtain


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle x(-1)=x_k
 * }
 * }

And


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle x(1)=x_{k+1}
 * }.
 * }.

Next using


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

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

we can transform


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle f(x(t))dx
 * }
 * }

into


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle g_k(t)dt
 * }
 * }

by using the derivative


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle \dfrac{d(x(t))}{dt}=\dfrac{h}{2}
 * }
 * }

Thus giving


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

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle d(x(t))=\dfrac{h}{2}dt
 * }.
 * }.

Placing all terms into the original equation with the new limits gives
 * {| style="width:100%" border="0" align="left"

$$ $$
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle 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))
 * }.
 * }.

Thus producing


 * {| 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= |

--Egm6341.s10.team1.toddmock 19:01, 5 March 2010 (UTC)
 * }
 * }

=Problem 5-Prove the Equality (integration by parts)=

Statement
Given $$ \displaystyle g_{k}(t):=\,f(x(t)):\,\,\,x\,\,\in\, [\,x_k,x_{k+1}\,]\,$$

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


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


 * Ref: Lecture Notes [[media:Egm6341.s10.mtg21.djvu|P21-2]]
 * }
 * }

Solution
Use Integration by parts [[media:Egm6341.s10.mtg5.pdf|P5-2]] $$ \displaystyle \int_{a}^{b}u^{'}vdt =\, [\,uv\,]_{a}^{b}-\,\int_{a}^{b}uv^{'}dt\,$$

Let $$ \displaystyle g^{(1)}(t)=\,u^{'}\,\,\,(-t)=\,v\,$$


 * $$ \displaystyle

\begin{align} \\\int_{-1}^{+1}g^{(1)}(t)(-t)dt &=\, [\,g(t)(-t)\,]_{-1}^{+1}-\,\int_{-1}^{+1}g(t)(-1)dt\; \\ &=\, -\,g(+1)-\,(g(-1))+\,\int_{-1}^{+1}g(t)dt\; \\ &=\, \int_{-1}^{+1}g(t)dt-\,[\,g(-1)+\,g(+1)\,]\,. \end{align} $$ EGM6341.s10.Team1.Ya-Chiao Chang 14:10, 7 March 2010 (UTC)

=Problem 6=

Find
Verify that $$g_k^{(i)}(t) = (\frac{h}{2})^i f^{(i)}(x(t)) \,$$ noting that $$g_k^{(i)}(t) = \frac{d^i}{dt^i} g_k(t) $$

Given
1) $$g_k(t) := f(x(t)), x \varepsilon [x_k, x_{k+1}]\,$$

2) $$x(t) := t\frac{h}{2}+\frac{x_k,x_{k+1}}{2}, t \varepsilon [-1,+1] $$ thus $$g_k(t) = f(t\frac{h}{2}+\frac{x_k,x_{k+1}}{2}) $$

Solution
$$g_k^{(1)}(t) = \frac{d}{dt} (\frac{h}{2}+\frac{x_k,x_{k+1}}{2}) f^{(1)}(t\frac{h}{2}+\frac{x_k,x_{k+1}}{2})\,$$

$$g_k^{(2)}(t) = [\frac{d}{dt} (\frac{h}{2}+\frac{x_k,x_{k+1}}{2})]^2 f^{(2)}(t\frac{h}{2}+\frac{x_k,x_{k+1}}{2})\,$$

$$g_k^{(3)}(t) = [\frac{d}{dt} (\frac{h}{2}+\frac{x_k,x_{k+1}}{2})]^3 f^{(3)}(t\frac{h}{2}+\frac{x_k,x_{k+1}}{2})\,$$

$$g_k^{(i)}(t) = [\frac{d}{dt} (\frac{h}{2}+\frac{x_k,x_{k+1}}{2})]^i f^{(i)}(t\frac{h}{2}+\frac{x_k,x_{k+1}}{2})\,$$

Since $$ \frac{d}{dt} (\frac{h}{2}+\frac{x_k,x_{k+1}}{2}) = \frac{h}{2}, \,$$

Egm6341.s10.team1.andrewdugan 21:48, 7 March 2010 (UTC)

=Problem 7=

Find
Verify that $$E = \int_{-1}^{+1}(-t)g^{(1)}(t)dt = [P_2(t)g^{(1)}(t)]_{-1}^{+1}-\int_{-1}^{+1}P_2(t)g^{(2)}(t)dt \,$$

Given
1) $$P_1(t):=-t \,$$ 2) $$P_2(t) :=\int P_1(t)dt $$

Solution
Using integration by parts, the integral becomes $$ \int_a^b udv = [udv]_a^b - \int_a^b vdu \,$$, where $$ u = g^{(1)}(t) \,$$ $$ du = g^{(2)}(t)dt \,$$ $$ v = P_2(t) \,$$ $$ dv = P_1(t)dt \,$$ $$ a = -1\,$$ $$ b = +1 \,$$

$$ \int_a^b udv = \int_{-1}^{+1}g^{(1)}(t)P_1(t)dt \,$$ $$ [udv]_a^b = [g^{(1)}(t)P_2(t)]_{-1}^{+1} \,$$ $$ \int_a^b vdu = \int_{-1}^{+1}P_2(t)g^{(2)}(t)dt \,$$

$$ \Rightarrow \int_{-1}^{+1}g^{(1)}(t)P_1(t)dt = [g^{(1)}(t)P_2(t)]_{-1}^{+1} - \int_{-1}^{+1}P_2(t)g^{(2)}(t)dt \,$$

Rembering that $$P_1(t)=-t \,$$,

Egm6341.s10.team1.andrewdugan 21:47, 7 March 2010 (UTC)

=Problem 8-Download and Install chebfun package=

Statement
Install chebfun, a Matlab code package wrote by Trefethen et al.(2004)


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


 * Ref: Lecture Notes [[media:Egm6341.s10.mtg25.djvu|P25-1]]
 * }
 * }

Install
Use this website to download Matlab package |chebfun

Test with two quick examples: (1)What is the integral of $$ \displaystyle sin(sin(x)) \,$$ from 0 to 10? (2)What is the maximum of $$ \displaystyle sin(x)+\,sin(x_{2})\,$$ over the same interval?

-> installed successfully

(3)calculate the integration of $$ \int_{0}^{1}\,\frac{e^{x}\,-1}{x}\mathrm{d} x $$

Use chebfun : Matlab code :

Get Ich = 1.317902151454403  Ech = 4.403366560268296e-012 (4) How long for using chebfun to calculate?

Get Elapsed time is 0.462024 seconds.

EGM6341.s10.Team1.Ya-Chiao Chang 14:10, 7 March 2010 (UTC)

EGM6341.s10.Team1.Ya-Chiao Chang 14:26, 7 March 2010 (UTC) - =Problem 9-Comparing speed of different algorithms=

Statement
1. Correct previous MATLAB code for efficiency 2. Compare the 2nd column of the Romberg table to the Simpson's rule. 3. Compute $$I=\int_{0}^{1}\frac{e^x-1}{x}dx$$ using the following methods and compare the execution times:


 * Composite trapezoidal rule
 * Composite Simpson's rule
 * Romberg table
 * chebfun sum command
 * MATLAB: trapz and quad commands

Solution
1. The following is the new (efficient) Romberg code:

2. The second column of the Romberg table is displayed below. The Simpson's rule is also shown (code shown below). The two methods are equivalent. The second column of the Romberg table is equal to the Simpson's rule.

3. The function $$I=\int_{0}^{1}\frac{e^x-1}{x}dx$$ is evualted. The computational time is computed in MATLAB using the tic and toc commands. The result is (code shown below):

EGM6341.s10.Team1.Kumanchik 21:54, 7 March 2010 (UTC)

= Problem 10: integration of ellipse arc length =

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

the arc length integral is:


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

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

compare MATLAB trapezoidal rule, Simpson's rule and Chebfun computational times.

Find
Total value of elliptical arc length (circumference) using built in MATLAB codes

Solution
For the Trapezoidal rule

For the Simpson's rule

And for the Chebfun function

As can be seen from the table above the quad function (simpson's rule) produces a result much faster than any of the other built in MATLAB functions --Egm6341.s10.team1.toddmock 19:02, 5 March 2010 (UTC)

=Problem 11-Determine computational time using various methods (Part 2)=

Statement
Consider the integral $$I=\int_{-5}^{5}\frac{1}{1+x^2}dx$$. Determine the computational time required to solve the integral numerically using,
 * Composite trapezoidal rule
 * Composite Simpson's rule
 * Romberg table
 * chebfun sum command
 * MATLAB: trapz and quad commands

Solution
The results are tabulated below (code also shown):

EGM6341.s10.Team1.Kumanchik 21:59, 7 March 2010 (UTC)

=Contributing members=

Egm6341.s10.team1.lei 19:05, 7 March 2010 (UTC) Author of Problems:1,2 and 3 Egm6341.s10.team1.andrewdugan 21:49, 7 March 2010 (UTC) Author of Problems 6,7 EGM6341.s10.Team1.Kumanchik 22:04, 7 March 2010 (UTC) Author problems 9-10 --Egm6341.s10.team1.toddmock 18:26, 9 March 2010 (UTC)