Posted: 08/23/21

Due: 08/30/21, 11:55PM

1. Construct a convergence plot in logarithmic coordinates for the Leibniz series approximation of $\pi$. Estimate the rate and order of convergence.

   Solution. The convergence of the Leibniz series partial sum sequence ${\left\{L_n\right\}}_{n\in \mathbb{N}}$

   ${\displaystyle L_n=\sum _{k=0}^n\frac{{(-1)}^k}{2k+1},{lim}_{n\to \infty }{L}_n=\frac{\pi }{4},}$

   is shown in Fig. 1, and can be visually estimated to converge with order $p=1$. The rate of convergence is

   ${\displaystyle r={lim}_{n\to \infty }\frac{L_{n+1}}{L_n}={lim}_{n\to \infty }(1+\frac{{(-1)}^{n+1}}{(2n+1)L_n})=1.}$

   $•$$\circ$

   Figure 1. Convergence of Leibniz series

   ∴	N=200; n=20:20:N;

   ∴	function Leibniz(n) k=0:n sum((-1).^k ./ (2*k.+1)) end;

   ∴	err=log.(abs.(Leibniz.(n) .- π/4));

   ∴	clf(); plot(log.(n),err,"-o"); xlabel("log(n)"); ylabel("log(err)"); title("Leibniz series convergence"); grid("on"); plot([3,4],[-5,-7],"-"); plot([3,4],[-5,-6],"-");

   ∴	savefig(homedir()*"/courses/MATH661/images/H01Fig01.eps")

   ∴

2. Apply Aitken acceleration to the Leibniz partial sums sequence. Construct the convergence plot. Estimate the rate and order of convergence.

   Solution. The Aitken-accelerated sequence ${\left\{A_n\right\}}_{n\in \mathbb{N}}$ has general term

   ${\displaystyle A_n=L_n-\frac{{(L_n-L_{n-1})}^2}{L_n-2L_{n-1}+L_{n-2}},}$

   with convergence plot shown in Fig. 2, indicating faster than quadratic growth.

   $•$$\circ$

   Figure 2. Convergence of Aitken acceleration of Leibniz series

   ∴	N=40; n=1:N; L=Leibniz.(n); A=copy(L);

   ∴	for n=3:N ΔL = L[n] - L[n-1] Δ2L = L[n] - 2*L[n-1] + L[n-2] A[n] = L[n] - ΔL^2/Δ2L end

   ∴	errL=log.(abs.(L .- π/4)); errA=log.(abs.(A .- π/4));

   ∴	clf(); plot(log.(n),errL,"-o"); plot(log.(n),errA,"-x"); xlabel("log(n)"); ylabel("log(err)"); title("Aitken-accelerated Leibniz series convergence"); grid("on"); plot([2,4],[-6,-8],"-"); plot([2,4],[-6,-10],"-");

   ∴	savefig(homedir()*"/courses/MATH661/images/H01Fig02.eps")

   ∴

3. Construct a convergence plot for the Ramanujan series approximation of $\pi$. Estimate the rate and order of convergence.

   Solution. The general term

   ${\displaystyle r_k=\frac{\left(4k\right)!(1103+26390k)}{{(k!)}^4{396}^{4k}}}$

   in the Ramanujuan series

   ${\displaystyle R_n=\frac{2\sqrt{2}}{9801}\sum _{k=0}^n\frac{\left(4k\right)!(1103+26390k)}{{(k!)}^4{396}^{4k}},{lim}_{n\to \infty }{R}_n=\frac{1}{\pi },}$

   decreases rapidly. Applying Stirling's formula, $\ln n!\cong n\ln -n$ (valid for large $n$), gives

   ${\displaystyle \ln r_k\cong \left(4k\right)\ln 4k-4k+\ln (1103+26390k)-4(k\ln k-k)-4k\ln 396,}$ ${\displaystyle \ln r_k\cong -4k\ln 99+\ln (1103+26390k)\Rightarrow }$ ${\displaystyle {\check{r}}_k={99}^{-4k}(1103+26390k)}$

   $•$$\circ$Since $r_0=1103\cong {10}^3$, terms smaller than ${10}^3ϵ$ ($ϵ$ is machine epsilon) would not appreaciably contribute to the sum

   ∴	ε=eps(Float64); kmax = log(10^3*ε)/(-4*log(99))

   $1.5851544170976106$

   ∴	ř(k)=99.0^(-4*k)*(1103+26390*k)

   ř

   ∴

   ř.(0:4)

   ${\displaystyle \left[\begin{array}{c}1103.0\\ 0.00028620772638853663\\ 5.839426693578139e-12\\ 9.056224058132163e-20\\ 1.2527103941143644e-27\end{array}\right]}$

   (1)

   ∴

   The above computation indicates that 2 terms of the sum should already give an approximation within machine precision. Floating point evaluation of $r_k$ should avoid intermediate computations close to overflow or undeflow. For example, the following evaluation is accurate in floating point arithmetic

   ∴	r(k) = (factorial(4*k)/(factorial(k))^4)*(1103.0+26390*k)/396.0^(4*k);

   ∴

   r.(0:4)

   ${\displaystyle \left[\begin{array}{c}1103.0\\ 2.6831974348925308e-5\\ 2.2453850201136646e-13\\ 1.9950749944958964e-21\\ 1.8393545314677564e-29\end{array}\right]}$

   (2)

   ∴

   In contrast, a different order of operations leads to incorrect results due to close-to-overflow conditions.

   ∴	rbad(k)=factorial(4*k)*(1103+26390*k)/(factorial(k))^4/396^(4*k);

   ∴	rbad.(0:4)

   ${\displaystyle \left[\begin{array}{c}1103.0\\ 2.6831974348925308e-5\\ -3.383976671091512e-11\\ 3.2409901671873466e-9\\ 8.916491639979272e-7\end{array}\right]}$

   (3)

   ∴

   Since there are only three partial sums that can be formed in floating point arithmetic, it is improper in this case to define an order of convergence based on the concept of a limit.

   ∴	R(n) = 2*sqrt(2.)/9801*sum(r.(0:n));

   ∴

   R.(0:2)

   ${\displaystyle \left[\begin{array}{c}0.31830987844047015\\ 0.31830988618379064\\ 0.3183098861837907\end{array}\right]}$

   (4)

   ∴	(R.(0:2) .- 1/π)*π

   ${\displaystyle \left[\begin{array}{c}-2.4326358923863425e-8\\ -1.743934249004316e-16\\ 0.0\end{array}\right]}$

   (5)

   ∴	eps(Float64)

   $2.220446049250313e-16$

   ∴

   From above computations, machine precision is achieved at $n=1$ (recall that machine epsilon is indistinguishable from zero, $1+ϵ=1$).

4. Construct a convergence plot for the continued fraction approximation of $\pi$. Estimate the rate and order of convergence.

   Solution. The continued fraction

   ${\displaystyle \pi +3=6+\frac{1^2}{6+},F_n=b_0+\stackrel{n}{\underset{k=1}{\text{K}}}\frac{a_k}{b_k},a_k={(2k-1)}^2,b_k=6,}$

   can be evaluated through multiple techniques:

   $•$$\circ$Iteration. The sequence of fractions is explicitly unrolled into a loop.

   ∴	function πIter(n) a(k) = (2*k-1)^2; b(k) = 6.0 F = b(n) for k = n:-1:1 F = b(k) + a(k)/F end return F end;

   ∴	πIter.(0:5)

   ${\displaystyle \left[\begin{array}{c}6.0\\ 6.166666666666667\\ 6.133333333333334\\ 6.145238095238096\\ 6.13968253968254\\ 6.142712842712843\end{array}\right]}$

   (6)

   ∴

   $•$$\circ$Recursion. The sequence of fractions is evaluated by a self-referential function. Note in the implementation the use of multi-method functions:

   • πRecur(k,n) is the recursive function with two arguments

   • πRecur(n) is a single-argument function that initializes the recursion

   ∴	function πRecur(k,n) a(k) = (2*k-1)^2; b(k) = 6.0 if (k==n) return b(n) else return b(k+1) + a(k+1)/πRecur(k+1,n) end end;

   ∴	πRecur(n) = πRecur(0,n);

   ∴	πRecur.(0:5)

   ${\displaystyle \left[\begin{array}{c}6.0\\ 6.166666666666667\\ 6.133333333333334\\ 6.145238095238096\\ 6.13968253968254\\ 6.142712842712843\end{array}\right]}$

   (7)

   ∴

   $•$$\circ$Composition. The continued fraction is evaluated through composition of functions

   ${\displaystyle T_n=t_0\circ t_1\circ \dots \circ t_n=\underset{k=0}{\overset{n}{\odot }}{t}_k=T_{n-1}\circ t_n}$

   ${\displaystyle t_0\left(z\right)=6+z,t_1\left(z\right)=\frac{1}{6+z},\dots ,t_k\left(z\right)=\frac{{(2k-1)}^2}{6+z},..}$

   ${\displaystyle F_n=T_n\left(0\right)}$

   ∴	function t(k,z) a(k) = (2*k-1)^2; b(k) = 6.0 if (k==0) return b(0) + z else return a(k)/(b(k)+z) end end;

   ∴	function T(n,z) if (n==0) return t(0,z) else return T(n-1,t(n,z)) end end;

   ∴	πComp(n) = T(n,0);

   ∴	n=5; [πIter(n) πRecur(n) πComp(n)]

   ${\displaystyle \left[\begin{array}{ccc}6.142712842712843& 6.142712842712843& 6.142712842712843\end{array}\right]}$

   (8)

   ∴

   From convergence study in Fig. 3, the continued fraction order of approximation is $p\cong 3$. This is often encountered: continued fraction approximations that lead to higher than second-order convergence.

   $•$$\circ$

   Figure 3. Convergence of continued fraction approximation of $\pi$

   ∴	N=200; n=20:20:N;

   ∴	err=log.(abs.(πIter.(n) .- (π+3)));

   ∴	clf(); plot(log.(n),err,"-o"); xlabel("log(n)"); ylabel("log(err)"); title("Continued fraction convergence"); grid("on"); plot([3,4],[-12,-13],"-"); plot([3,4],[-12,-14],"-"); plot([3,4],[-12,-15],"-");

   ∴	savefig(homedir()*"/courses/MATH661/images/H01Fig03.eps")

   ∴

5. Completely state the mathematical problem of computing the absolute value

   ${\displaystyle \left|x\right|,x\in ℝ.}$

   Find the absolute and relative condition numbers.

   Solution. The mathematical problem is

   ${\displaystyle f:ℝ\to {ℝ}_{+},f\left(x\right)=\left|x\right|}$

   From the overview of mathematical problem conditioning [1], the absolute condition number introduced in [2] is given by

   ${\displaystyle \hat{\kappa }=\left|f^{\text{'}}\left(x\right)\right|=1}$

   where $f$ is differentiable. The relative condition number is also equal to one

   ${\displaystyle \kappa =\frac{\left|x\right|\left|f^{\text{'}}\left(x\right)\right|}{\left|f\left(x\right)\right|}=1.}$

   The only remaining point is behavior at $x=0$. By definition

   ${\displaystyle \hat{\kappa }={lim}_{\epsilon \to 0}{sup}_{\left|\delta x\right|⩽\epsilon }\frac{|f(0+\delta x)-f\left(0\right)|}{\left|\delta x\right|}={lim}_{\epsilon \to 0}{sup}_{\left|\delta x\right|⩽\epsilon }\frac{\left|\delta x\right|}{\left|\delta x\right|}=1.}$

   The relative condition number is

   ${\displaystyle \hat{\kappa }={lim}_{\epsilon \to 0}{sup}_{\left|\delta x\right|⩽\epsilon }\frac{|f(0+\delta x)-f\left(0\right)|\left|\epsilon \right|}{\left|\delta x\right|\left|f\left(\epsilon \right)\right|}=1.}$

   Taking the absolute value is a well-conditioned mathematical problem.

6. Completely state the mathematical problem of computing a difference.

   ${\displaystyle x_1-x_2,x_1,x_2\in ℝ.}$

   Find the absolute and relative condition numbers.

   Solution. The mathematical problem is

   ${\displaystyle f:ℝ\times ℝ\to ℝ,f\left(x\right)=x_1-x_2,}$

   with absolute condition number

   ${\displaystyle \hat{\kappa }\left(𝒙\right)={lim}_{\epsilon \to 0}{sup}_{||\delta 𝒙||⩽\epsilon }\frac{|\delta x_1-\delta x_2|}{||\delta 𝒙||}.}$

   Choosing the 1-norm,

   ${\displaystyle \hat{\kappa }\left(𝒙\right)={lim}_{\epsilon \to 0}{sup}_{||\delta 𝒙||⩽\epsilon }\frac{|\delta x_1-\delta x_2|}{\left|\delta x_1\right|+\left|\delta x_2\right|}⩽{lim}_{\epsilon \to 0}{sup}_{||\delta 𝒙||⩽\epsilon }\frac{|\delta x_1|+|\delta x_2|}{\left|\delta x_1\right|+\left|\delta x_2\right|}=1.}$

   The relative condition number is

   ${\displaystyle \hat{\kappa }\left(𝒙\right)={lim}_{\epsilon \to 0}{sup}_{||\delta 𝒙||⩽\epsilon }\frac{|\delta x_1-\delta x_2|}{\left|\delta x_1\right|+\left|\delta x_2\right|}\cdot \frac{\left|x_1\right|+\left|x_2\right|}{|x_1-x_2|}.}$

   Note that the relative condition number becomes arbitrarily large close to the $x_1=x_2$ bisector. Interpret the above to signify:

   • When $x_1\cong x_2\Rightarrow x_1-x_2\cong 0$, and an absolute condition number of $\hat{\kappa }=1$ indicates that errors in the input are not amplified;

   • However, relative to the expected result $x_1-x_2\cong 0$, relative errors are not bounded

     ${\displaystyle \hat{\kappa }\sim 𝒪\left(\frac{1}{|x_1-x_2|}\right).}$

Bibliography