User:Egm6341.s10.team3.ynahn81/HW6-1

= (1) Computation of arc length an ellipse $$\displaystyle \overset{\frown}{PQ} $$ using successive numerical integration = Ref: Lecture notes [[media:Egm6341.s10.mtg32.djvu|p.32-1]]

Problem Statement
Compute arc length of ellipse $$\displaystyle \overset{\frown}{PQ} $$ to $$\displaystyle O(10^{-10}) $$ from $$\displaystyle \theta_{P} = 0 $$ to $$\displaystyle \theta_{Q} = \pi/3 $$. Successive numerical integration results should be used as stopping criterion, i.e. $$\displaystyle |I_{2n} - I_{n}| < O(10^{-10}) $$. 

i) Use error estimate for composite trapezoidal rule
 * {| style="width:100%" border="0" align="left"

$$ $$ to find $$\displaystyle n $$.  
 * $$\displaystyle
 * E_{n}| \leq \frac{(b-a)h^2}{12}M_{2}
 * E_{n}| \leq \frac{(b-a)h^2}{12}M_{2}
 * $$\displaystyle (Eq. 1)
 * }
 * }

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

\overset{\frown}{PQ} = \int_{\theta_{P}}^{\theta_{Q}} dl $$ $$ where $$\displaystyle dl = d\theta\left[r^2 + \left(\frac{dr}{d\theta}\right)^2\right]^{\frac{1}{2}} $$ and $$\displaystyle r(\theta) = \frac{(1-e^2)}{1-e\cos\theta} $$ to obtain the arc length from composite trapezoidal rule, Romberg table, clencurt in Matlab and sum in Chebfun.  
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 2)
 * }
 * }

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

\overset{\frown}{PQ} = \int_{\theta_{P}}^{\theta_{Q}} \left[1-e^2\sin^2\theta\right]^{\frac{1}{2}}d\theta. $$ $$ to obtain the arc length from composite trapezoidal rule, Romberg table, clencurt in Matlab and sum in Chebfun.  
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle (Eq. 3)
 * }
 * }

Solution
i) Find $$\displaystyle n $$ In (Eq. 1), recall that $$\displaystyle h = \frac{b-a}{n} $$ and $$\displaystyle M_{2} = \max_{\zeta \in [a,b]}|f^{(2)}(\zeta)| $$. Here, $$\displaystyle a = 0 $$ and $$\displaystyle b = \pi/3 $$.

Also, rearrange (Eq. 1) to get the following equation.
 * {| style="width:100%" border="0" align="left"

n^2 \leq \frac{(b-a)^3}{12}\frac{\max_{\zeta \in [a,b]}|f^{(2)}(\zeta)|}{|E_{n}|}. $$ $$   Now, put the below values into (Eq. 4).  $$\displaystyle a = 0 $$  <br\> $$\displaystyle b = \pi/3 $$ <br\> <br\> $$\displaystyle |E_{n}| = 10^{-11} (\because |I_{2n} - I_{n}| < O(10^{-11}) $$ <br\> <br\> $$\displaystyle M_{2} = \max_{\zeta \in [a,b]}|f^{(2)}(\zeta)| = \max_{\theta \in [0,\frac{\pi}{3}]}|dl^{(2)}(\theta)| = 0.2934 $$ for (Eq. 2) and <br\> <br\> $$\displaystyle M_{2} = \max_{\zeta \in [a,b]}|f^{(2)}(\zeta)| = \max_{\theta \in [0,\frac{\pi}{3}]}|\left(\sqrt{1-e^2\sin^2(\theta)}\right)^{(2)}| = 0.0670 $$ for (Eq. 3) <br\> <br\> <br\>
 * $$\displaystyle
 * $$\displaystyle
 * <p style="text-align:right;">$$\displaystyle (Eq. 4)
 * }
 * }

Then, $$\displaystyle n^2 \leq 2.8077e+09 \quad \Rightarrow \quad n \leq 52988.59 $$ for (Eq. 2) and<br\> <br\> $$\displaystyle n^2 \leq 6.4105e+08 \quad \Rightarrow \quad n \leq 25319.00 $$ for (Eq. 3).<br\> <br\> <br\>

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

$$\displaystyle n \leq 53988 $$ for (Eq. 2) and
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = |
 * }
 * }
 * {| style="width:100%" border="0" align="left"

$$\displaystyle n \leq 25319 $$ for (Eq. 3). <br\> <br\> <br\>
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style = |
 * }
 * }

ii) Based on (Eq. 2) In what follows, we are using the follow matlab function as a given function $$\displaystyle dl(\theta) $$.

<br\> <br\>

 (a) Composite trapezoidal rule  <br\> Recall the formula of composite trapezoidal 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}). $$ <br\> <br\> <br\> <br\>
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

<br\> The obtained number of nodes to achieve $$\displaystyle O(10^{-10}) $$ is $$\displaystyle n = 32768 $$ which satisfies the result of part i) $$\displaystyle n \leq 53988 $$. <br\> The integration result, final error and computational time are summarized in Table at the last part of this problem. <br\> <br\>

 (b) Romberg Table  <br\>

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

Result

<br\>

The integration result, final error and computational time are summarized in Table at the last part of this problem.

<br\> <br\>

 (c) clencurt  <br\>

To achieve the desired tolerance of error, we use the following Matlab code. <br\>

The integration result, final error, final number of nodes and computational time are summarized in Table at the last part of this problem. <br\> <br\>

 (d) chebfun sum command  <br\> The 'sum' command in chebfun performs integration of a given function up to the machine accuracy. Therefore, we can achieve the tolerance of error of numerical integration for the given function $$\displaystyle dl $$ up to the machine accuracy by using the following Matlab code.

The integration result, final error, final number of nodes and computational time are summarized in Table at the last part of this problem. <br\> <br\>

 (e) Summary  <br\>

In the previous HW, recall that composite trapezoidal rule was a very efficient method to perform integration due to the periodicity of the given function. In this HW, however, note that composite trapezoidal rule requires a very large number of nodes and long computation time. The reason is that the given function does not show a periodicity within a given range $$\displaystyle \theta \in [0, \pi/3] $$  even though it shows a periodicity for the range of $$\displaystyle \theta \in [0, 2\pi] $$.

As shown above, both Romberg table and clencurt require only 32 nodes to achieve the desired tolerance of error and the computation times of the both methods are very short. Therefore, both Romberg Table and clencurt are efficient methods to integrate the given function. <br\> <br\> <br\>

iii) Based on (Eq. 3) In what follows, we are using the follow matlab function as a given function $$\displaystyle dl(\theta) $$.

<br\> <br\>

 (a) Composite trapezoidal rule  <br\> We use the same Matlab code with that of part ii).

The obtained number of nodes to achieve $$\displaystyle O(10^{-10}) $$ is $$\displaystyle n = 16384 $$ which satisfies the result of part i) $$\displaystyle n \leq 25319 $$. <br\> The integration result, final error and computational time are summarized in Table at the last part of this problem. <br\> <br\>

 (b) Romberg Table  <br\> We use the same Matlab code with that of part ii).

<br\>

The integration result, final error and computational time are summarized in Table at the last part of this problem.

<br\> <br\>

 (c) clencurt  <br\> We use the same Matlab code with that of part ii).

The integration result, final error, final number of nodes and computational time are summarized in Table at the last part of this problem. <br\> <br\>

 (d) chebfun sum command  <br\> The 'sum' command in chebfun performs integration of a given function up to the machine accuracy. Therefore, we can achieve the tolerance of error of numerical integration for the given function $$\displaystyle sqrt(1 - e^2\sin^2\theta) $$ up to the machine accuracy by using the following Matlab code.

The integration result, final error, final number of nodes and computational time are summarized in Table at the last part of this problem. <br\> <br\>

 (e) Summary  <br\>

As we already explained in summary of part i), composite trapezoidal rule is not a efficient method to integrate the given function due to the same reason (i.e. The given function does not show a periodicity within a given range $$\displaystyle \theta \in [0, \pi/3] $$). Romberg table and clencurt require only 16 nodes to achieve the desired tolerance of error and the computation times of the both methods are very short. Therefore, both Romberg Table and clencurt are efficient methods to integrate the given function.

In addition, note that the integration result of part iii) is different from the integration result of part ii) although the given ranges of the integrations are with each other (i.e. $$\displaystyle \theta \in [0, \pi/3] $$). As the figures show in problem statement, the definition of $$\displaystyle \theta $$ of part ii) is different from that of part iii). Therefore, it is very clear that the integration results of the both parts are not same with each other even though the range of $$\displaystyle \theta $$ is same.