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}^{b}dl=\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$

% Use error estimate of composite trapezoidal rule
% to find n
syms x;
e=sin(pi/12);
% Integrating function
F=sym((1-e^2)/(1-e*cos(x))^2*(1-2*e*cos(x)+e^2)^0.5);
% Second deriv of F
ddF=diff(F,'x',2);
a=0;
b=60*pi/180;
X=(a:.01:b)';
% Y=ddF is evaluated at X
Y=subs(ddF,x,X);
% Max determined
M2=max(abs(Y));
% Upper bound on n
n=((b-a)^3/12/1e-10*M2)^0.5;

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.

% Compositve Trapezoidal technique
% start n=2
n=2;
h=(b-a)/n;
xn=a:h:b;
% evaluate the function F at all nodes
fn=subs(F,x,xn);
% sum these 3 nodes (n=2)
pre=0.5*(fn(1)+fn(n+1))+sum(fn(2:n));
% divide by stepsize to get the integration
Tn=h*pre;
while 1
    % double the n
    n2n=2*n;
    h=(b-a)/n2n;
    % calculate all new nodes (not the ones we previously calculated)
    x2n=xn(1:n)+h;
    % evaluate the function F at these new nodes
    f2n=subs(F,x,x2n);
    % sum up these new nodes with the previously summed nodes
    pre2n=pre+sum(f2n);
    % calculate the new integration
    T2n=h*pre2n;
    % compare to last integration
    if abs(T2n-Tn) < 1e-10
        % job done
        break;
    end
    % reset and try again
    Tn=T2n;
    n=n2n;
    xn=a:h:b;
    pre=pre2n;
end

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.

% start at n=2^2=4
n=2;
% calculate the Romberg table
tbl=myRomb(f,a,b,2^n);
% store the most accurate result in the table
V=tbl(1,n);
while 1
    n=n+1;
    % calulate a new Romberg table
    tbl=myRomb(f,a,b,2^n);
    Vo=tbl(1,n);
    % compare most accurate result of this table to old table
    if abs(V-Vo) < 10^-10
        % job done
        break;
    end
    % reset and try again
    V=Vo;
end

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.

fc=chebfun(f,[a b]);
chbans=sum(fc);

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})$

% clencurt
% start at n=1
n=1;
% compute teh clencurt for n=1
[xc,w]=clencurt(n);
xc=0.5*((xc+1)*a-(xc-1)*b);
w=0.5*(b-a)*w;
Y=subs(F,x,xc);
% evaluate the integral for n=1
C=sum(w'.*Y);
while 1;
    % increment n
    n=n+1;
    [xc,w]=clencurt(n);
    xc=0.5*((xc+1)*a-(xc-1)*b);
    w=0.5*(b-a)*w;
    Y=subs(F,x,xc);
    % evaluate the new integral
    Co=sum(w'.*Y);
    % compare to previous result
    if abs(Co-C)<10^-10
        % job done
        break;
    end
    % reset and try again
    C=Co;
end

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:

Kessler table

	numerator	denominator		numerator	denominator
$c_{1}$	-1	1		-1	3	$p_{2}(1)$
$c_{3}$	1	6		1	45	$p_{4}(1)$
$c_{5}$	-7	360		-2	945	$p_{6}(1)$
$c_{7}$	31	15120		1	4725	$p_{8}(1)$
$c_{9}$	-127	604800		-2	93555	$p_{10}(1)$
$c_{11}$	73	3421440		1382	638512875	$p_{12}(1)$
$c_{13}$	-1414477	653837184000		-4	18243225	$p_{14}(1)$
$c_{15}$	8191	37362124800		3617	162820783125	$p_{16}(1)$
$c_{17}$	-16931177	762187345920000

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.

$k$	$c_{2k+1}$	$p_{2k}$
1	${\frac {c_{1}}{3!}}+{\frac {c_{3}}{1!}}=0$	$p_{2}(1)={\frac {c_{1}}{2!}}+{\frac {c_{3}}{0!}}$
2	${\frac {c_{1}}{5!}}+{\frac {c_{3}}{3!}}+{\frac {c_{5}}{1!}}=0$	$p_{4}(1)={\frac {c_{1}}{4!}}+{\frac {c_{3}}{2!}}+{\frac {c_{5}}{0!}}$
3	${\frac {c_{1}}{7!}}+{\frac {c_{3}}{5!}}+{\frac {c_{5}}{3!}}+{\frac {c_{7}}{1!}}=0$	$p_{6}(1)={\frac {c_{1}}{6!}}+{\frac {c_{3}}{4!}}+{\frac {c_{5}}{2!}}+{\frac {c_{7}}{0!}}$

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:

function [c,p]=traperror(n)
%compute the coefficients p_2, p_4, ... , p_{2n} associated with the
%trapezoidal rule error expansion. Also compute c_1, c_3, ...
f=1; g=2;
% f will be built up term-by-term to represent the odd-numbered factorials that appear in the denominators of the coefficient equation
% g will be used in conjunction with f to divide out the last multiplier in the f factorial thus making it an even factorial
cn=-1; cd=1;
% cn and cd are the numerators and denominators of all the coefficients computed.
% Since c1 = -1, this term is initialized here.
for k=1:n
% loop starts
f=f.*g.*(g+1);
% Though seemingly complex, this progressively builds f into odd-numbered factorials for computing the coefficients.
% Note that it will build the factorial vector in descending order.
[newcn,newcd]=fracsum(-1*cn,cd.*f);
% First the fraction is negated (moved to the right-hand-side of the equation).
% fracsum sums the inputted fraction.
% It is more complex so that it preserves the whole number ratio of the fraction rather than computing a decimal number.
cn=[cn;newcn]; cd=[cd;newcd];
% The result is inputted into the coefficient vectors (numerator and denominator respectively).
% Note this builds a vector that represents the coefficients in ascending order.
f=[f;1];
% Here the coefficient factorial vector is updated for use with g in this iteration
g=[2+g(1);g];
% Here the g vector is updated to be used with f next
[newpn,newpd]=fracsum((g-1).*cn,f.*cd);
% Here the polynomial is computed by adding all the terms (coefficients multiplied by the factorial).
% g is adjusted here to properly cancel the highest multiplier on f to create the proper factorial for use in the polynomial equation
pn(k,1)=newpn; pd(k,1)=newpd;
% The polynomial of the k-th term is stored here as numerator and denominator rather than as a decimal.
end
c=int64([cn cd]);
p=int64([pn pd]);
This operation is for display purposes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [nsum,dsum]=fracsum(n,d)
% n - numerator, d - denominator. n./d are all the fractions to be summed
div=gcd(round(n),round(d));
% Here the Greatest Common Denominator is found to simplify the fraction
n=round(n./div);
d=round(d./div);
% This does the simplification
dsum=1;
% dsum - represents the least common multiple of every single fraction to be summed
for k=1:length(d)
dsum=lcm(dsum,d(k));
end
% Mission complete, dsum is now the least common multiple of all fractions to be summed
nsum=dsum*sum(n./d);
% If a fraction is represented as Num/Den = Decimal, then the Num = Den*Decimal.
% Since dsum is the least common multiple of all fractions summed, nsum is guaranteed to be a whole number
div=gcd(round(nsum),round(dsum));
% The sum may have a common divisor so it is found here
nsum=nsum/div;
% The numerator is divided
dsum=dsum/div;
% The denominator is divided
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%Note the last line had to be added to generate the correct table
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

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_{3}t$

$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_{5}t$

$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_{7}t$

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.