User:EGM6341.s10.Team1.Kumanchik/HW6

=Problem 1: Compute the arc length of an ellipse on $$[0^{\circ},60^{\circ}]$$=

Statement
An arc is denoted $$\widehat{PQ}$$ and lies on an ellipse defined by,
 * $$r(\theta)=\frac{a(1-e^2)}{1-e \cos{(\theta)}}$$

where $$a=1$$ and $$e=\sin{(\frac{\pi}{12})}$$. The equation for the arc length is,
 * $$PQ=\int_a^bdl=\int_{\theta_1}^{\theta_2}d\theta \sqrt{r^2+\left (\frac{dr}{d\theta} \right )^2}$$

1.) Estimate the number of nodes required for the composite trapezoidal estimate of the integral to have an error $$O(10^{-10})$$ 2.)  Using successive numerical integrations, find the actual number of nodes required for the trapezoidal technique by considering the stopping criterion $$\left |I_{2n}-I_n \right |<10^{-10}$$ 3.) How many nodes are required for the Romberg table (verify vs trapezoidal), clencurt.m (verify with a similar stopping criterion), and sum (chebfun - compare to trapezoidal).

Solution
The integrating function is,
 * $$F=\sqrt{r^2+\left (\frac{dr}{d\theta} \right )^2}=\frac{1-e^2}{(1-e\cos{\theta})^2}\sqrt{1-2e\cos{\theta}+e^2}$$

1.) An estimate for the error in the composite trapezoidal technique is,
 * $$\left | E_n \right |\leq \frac{(b-a)^3}{12n^2}M_2$$

Rearranging for the number of nodes, $$n$$, gives,
 * $$n=\sqrt{\frac{(b-a)^3}{12(10^{-10})}M_2}$$

where $$M_2=max\left | F^{(2)}(\zeta) \right |, \zeta\in [a,b]$$ and $$[a,b]=[0,60\frac{\pi}{180}]$$. The following MATLAB code solves for $$M_2$$ and determines $$n<16,755$$

2.) Numerically integrating $$F$$ using the trapezoidal rule detemines that at $$n=16,384$$ we get $$\left | I_{16,384}-I_{8,192} \right | < 10^{-10}$$. The following MATLAB code is used.  Notice that in step sizes of 2n, the previous summation can be utilized as shown in the code.

3.) Numerically integrate using the Romberg, clencurt, and sum (chebfun).

Romberg gives $$O(10^{-10})$$ compared to itself at $$n=2^7=128$$, see code.

For the sum (chebfun) the code is shown next. Comparing this to the composite trapezoidal result shows a difference $$O(10^{-11})$$ so is considered accurate enough.

And finally for clencurt.m is shown last. The result is compared to itself and at n=10 has an error less than $$O(10^{-10})$$

=Problem 10: Use Kessler's code to generate his table=

Statement
Use Kessler's code to generate his table and then follow the code line-by-line to get $$(P_2,P_3)(P_4,P_5)(P_6,P_6)$$. Comment on the code.

Solution
Kessler's table was generated with the code below for n=8:

The code computes the Trapezoidal error polynomials evaluated at t=1. The polynomial equation is,
 * $$p_{2k}(t)=c_1\frac{t^{2k}}{(2k)!}+c_3\frac{t^{2k-2}}{(2k-2)!}+c_5\frac{t^{2k-4}}{(2k-4)!}+...+c_{2k-3}\frac{t^4}{4!}+c_{2k-1}\frac{t^2}{2!}+c_{2k+1}\frac{t^0}{0!}$$

When evaluated at t=1, the equation is,
 * $$p_{2k}(1)=\frac{c_1}{(2k)!}+\frac{c_3}{(2k-2)!}+\frac{c_5}{(2k-4)!}+...+\frac{c_{2k-3}}{4!}+\frac{c_{2k-1}}{2!}+\frac{c_{2k+1}}{0!}$$

The coefficients, $$c$$, are computed by solving the below equation using $$c_1=-1$$ as a starting point,
 * $$\frac{c_1}{(2k-1)!}+\frac{c_3}{(2k-3)!}+\frac{c_5}{(2k-5)!}+...+\frac{c_{2k-5}}{5!}+\frac{c_{2k-3}}{3!}+\frac{c_{2k-1}}{1!}=0$$

For illustration, the first few coefficients and polynomial equations are shown below.

Notice that the factorials in the denominators of the polynomials differ by 1 from those in the coefficients. Also, to solve for the next coefficient, the preceding coefficients only need to be summed and then moved to the other side of the equal sign, i.e $$c_7=-(\frac{c_1}{7!}+\frac{c_3}{5!}+\frac{c_5}{3!})$$.

The following code will be commented line-by-line to explain the algorithm of Kessler:

Finally, the code will be followed to generate $$(P_2,P_3)(P_4,P_5)(P_6,P_6)$$: The polynomials are written as,
 * $$p_2(t)=c_1\frac{t^2}{2!}+c_3$$
 * $$p_3(t)=c_1\frac{t^3}{3!}+c_3t$$
 * $$p_4(t)=c_1\frac{t^4}{4!}+c_3\frac{t^2}{2!}+c_5$$
 * $$p_5(t)=c_1\frac{t^5}{5!}+c_3\frac{t^3}{3!}+c_5t$$
 * $$p_6(t)=c_1\frac{t^6}{6!}+c_3\frac{t^4}{4!}+c_5\frac{t^2}{2!}+c_7$$
 * $$p_7(t)=c_1\frac{t^7}{7!}+c_3\frac{t^5}{5!}+c_5\frac{t^3}{3!}+c_7t$$

Following the code to get $$c_{1-7}$$, n=3 because 3 loops will generate all 7 coefficients f=1, g=2, cn=-1, cd=1 entering loop which runs 3 times f=1*2*3 (3! = 6) enter fracsum with (1,6) and will sum the fractions w/o leaving a decimal div = gcd(1,6) = 1 n = 1 d = 6 dsum = 1 enter loop (runs once) dsum = 6 leave loop nsum = 6*(1/6)= 1 div = gcd (1,6) = 1 nsum = 1 dsum = 6 leave fracsum with (1,6) cn = [-1; 1] cd = [1; 6] f = [6;1] g = [4;2] enter fracsum with ([-3;1],[6;6]) div = [3;1] n=[-1;1] d=[2;6] enter loop dsum = 2 then 6 leave loop nsum=6*-1/2+6*1/6=-3+1=-2; div=2 nsum=-1; dsum=3; leave fracsum with (-1,3) pn(1,1)=-1; pd(1,1)=3; move on to next loop k=2

At this point it is apparent that fracsum will sum up fractions without reducing them to decimals Also, the second fracsum call does not affect cn so we will skip calculating it And, cn=[-1;1] cd=[1;6] so we have c1 and c3

now k=2 f=[6;1].*[4;2].*[5;3] = [24;2].*[5;3] = [120;6] enter fracsum with (-[-1;1],[1;6].*[120;6])=([1;-1],[120;36]) this is like 1/120-1/36 = -7/360 leave fracsum cn = [-1;1;-7]; cd = [1;6;360]; f = [120;6;1] g = [2+4;[4;2]] = [6;4;2] skip to k=3 f=[120;6;1].*[6;4;2].*[7;5;3] = [720;24;2].*[7;5;3] = [5040;120;6] enter fracsum with ([1;-1;7],[1;6;360].*[5040;120;6]) = ([1;-1;7],[5040;720;2160]) this is 1/5040-1/720+6/2160 = 31/15120 leave fracsum cn = [-1;1;-7;31]; cd = [1;6;360;15120]

The result is,
 * $$c_1=-1 \!$$
 * $$c_3=-\frac{1}{6}$$
 * $$c_5=-\frac{-7}{360}$$
 * $$c_7=-\frac{31}{15120}$$

Plugging into the expressions for $$p_1(t) -> p_7(t)$$ generates the correct result.