%% Richardson Extrapolation
%% (and the Shanks Transformation)
%
%% First Iteration
% Lewis Richardson (1881-1953), first to apply mathematics to weather prediction, sought methods to generate higher accuracy based on low-order formulae. Similar methods were used by Archimedes and Huygens. [See Burden and Faires, Numerical Analysis.]
% We want to compute M from formula N_1(h). Assuming that M-N_1(h) = K_1h+K_2h^2+O(h^3), we seek a higher order approximation,\hat N(h), such that M=\hat N(h) + \hat K_2h^2+O(h^3).
% Replace h by \frac h2 to obtain M=N_1\left(\frac h2\right) + \frac 12 K_1h+\frac 14K_2h^2+O(h^3). Note that M = 2M-M = 2N_1\left(\frac h2\right) - N_1(h) - \frac 12
%  K_2h^2+O(h^3). This gives a an approximation of second order, which we can write as  N_2(h) = 2N_1\left(\frac h2\right) - N_1(h)  = N_1\left(\frac h2\right)+\left[N_1\left(\frac h2\right)-N_1(h)\right].
% Example: We apply this to improve the centered finite difference approximation.

clear
f = @(x) cos(x.^2);
dfdx = @(x) -2*x.*sin(x.^2);

x = 1;
exact=dfdx(x)

%
% Compute the centered difference formula followed by extrapolation N_2(h).

% Apply centered difference
h = 0.5;
df(1,1) = (f(x + h) - f(x - h))/(2*h)    % N1(h)
err=abs(exact-df(1,1))

% Extrapolate Once
h=h/2;
df(2,1) =  (f(x + h) - f(x - h))/(2*h);  % N1(h/2)
df(2,2)=(2*df(2,1)-df(1,1))              % N2(h)
err=abs(exact-df(2,2))

%
% Do another extrapolation to find N_3(h).

h=h/2;
df(3,1) =  (f(x + h) - f(x - h))/(2*h);  % N1(h/4)
df(3,2)=2*df(3,1)-df(2,1);               % N2(h/2)
df(3,3)= df(3,2) +1/3*(df(3,2)-df(2,2))  % N3(h)
%df(3,3)=(4*df(3,2)-df(2,2))/3
err=abs(exact-df(3,3))

%
%% General Iteration
% We can repeat the process.  Starting with M = N_2(h) - \frac 12
%  K_2h^2-\frac 34K_3 h^3+O(h^4), We eliminate the h^2 term by halving h. This gives M = N_2\left(\frac h2\right) - \frac 18
%  K_2h^2-\frac 3{32}K_3 h^3+O(h^4). Eliminating h^2, we find 3M =4 N_2\left(\frac h2\right) -N_2(h) + \frac 3{8}K_3 h^3+O(h^4). So, we identify N_3(h)=N_2\left(\frac h2\right) +\frac 13\left[N_2\left(\frac h2\right) -N_2(h)\right].
% Carrying this on, we have the general process N_j(h)=N_{j-1}\left(\frac h2\right) +\frac 1{2 ^{j-1}-1}\left[N_{j-1}\left(\frac h2\right) -N_{j-1}(h)\right].
%
% Let Q_{i,j}=N_i(2^{j-i}h). We can place these in an array:
% \matrix{ Q_{1,1} &  & \cr Q_{2,1} & Q_{2,2} & \cr Q_{3,1} & Q_{3,2} & Q_{3,3} \cr Q_{4,1} & Q_{4,2}& Q_{4,3} & Q_{4,4} }
% We first compute the first column as Q_{i,1}=N_i(2^{1-i}h).. Successive columns are found from the previous column using extrapolation. The improved answers are found along the diagonal.

h=0.5;
f = @(x) cos(x.^2);
dfdx = @(x) -2*x.*sin(x.^2);

x = 1;
exact=dfdx(x)

% First column
for j=1:7
  H=h/2^(j-1);
  F(j)=(f(x+H)-f(x-H))/(2*H);
end

% Put extrapolations in Q
n=length(F);
Q=zeros(n,n);
Q(:,1)=F(:);

for col=2:n
  for row=col:n
    Q(row,col)=(2^(col-1)*Q(row,col-1)-Q(row-1,col-1))/(2^(col-1)-1);
  end
end

% Result
Q

% Plot Error
Error=diag(abs(Q-exact));

plot(log10(Error),'.','MarkerSize',20);
title('Log_{10} of the Error vs. Extrapolation Number')

%
% Example: We apply Richardson extrapolation to improve the Trapezoidal approximaton.

f=@(x) exp(sin(x));
a=0;
b=2;
n=4;
x = linspace(a,b,201)';
fill([x;b;a],[f(x);0;0],[1,0.9,0.9])
title('Area Under f(x)=exp(sin(x))')  
xlabel('x'), ylabel('f(x)')    

I = integral(f,a,b,'abstol',1e-14,'reltol',1e-14);
fprintf('Integral =    %.15f\n',I)
T=trapezoid(f,a,b,n);
fprintf('Trapezoidal = %.15f\n',T)
fprintf('Error       = %.15f\n',I-T)

%
% Compute the first column using Trapezoid by double the number of subintervals.

% First column
for j=1:7
  N=n*2^(j-1);
  T(j)=trapezoid(f,a,b,N);
end

%
% Fill in the rest of the extrapolation matrix and plot the error.

% Put extrapolations in Q
M=length(T);
Q=zeros(M,M);
Q(:,1)=T(:);

for col=2:M
  for row=col:M
    Q(row,col)=(2^(col-1)*Q(row,col-1)-Q(row-1,col-1))/(2^(col-1)-1);
  end
end

% Result
Q
% Plot Error
Error=diag(abs(Q-I));

plot(log10(Error),'.','MarkerSize',20);
title('Log_{10} of the Error vs. Extrapolation Number')

%
% It is interesting to see the errors in the Trapezoidal formula.

Error=diag(abs(T-I));

plot(n*2.^(0:6),log10(Error),'.','MarkerSize',20);
title('Log{10} of the Error vs. Number of Subintervals')

%
%% Shanks Transformation
% According to Wikipedia, "the Shanks transformation is a non-linear series acceleration method to increase the rate of convergence of a sequence. This method is named after Daniel Shanks, who rediscovered this sequence transformation in 1955." It is usually used to speed up convergence of partial sums of infinite series. So, if A_n=\sum_{k=0}^n a_k, then the new sequence given by
% S(A_n)=\frac{A_{n+1}A_{n-1}-A_n^2}{A_{n+1}-2A_n+A_{n-1}} = A_{n+1}-\frac{(A_{n+1}-A_n)^2}{(A_{n+1}-A_n)-(A_n-A_{n-1})} often converges more rapidly.
%
% Example: Wikipedia gives the Gregory series for \pi = 4\sum_{k=0}^\infty \frac{(-1)^k}{2k+1}.  Six-figure accuracy requires summing about 400,000 terms. We apply the Shanks transformaton several time to show we get this accuarcy with only a few terms.

N=20;
Ncol=4;
psum(1) =4;
for k=1:N
    psum(k+1)=psum(k) +4*(-1)^k/(2*k+1);
end

format long
psum

M=length(psum);
S=zeros(M,Ncol);
S(:,1)=psum(:);

for col=2:Ncol
  for row=col:M-(col-1)
   % S(row,col)=(S(row+1,col-1)*S(row-1,col-1)-S(row,col-1)^2)/...
   %     (S(row+1,col-1)-2*S(row,col-1)+S(row-1,col-1));
    S(row,col)=S(row+1,col-1)-(S(row+1,col-1)-S(row,col-1))^2/...
        ((S(row+1,col-1)-S(row,col-1))-(S(row,col-1)-S(row-1,col-1)));
  end
end

% Result
S

%
% We plot the logarithm of the absolute error for three applications of this transform.

AbsError =abs(pi-psum); 
plot(log10(AbsError),'.','MarkerSize',20);
title('Log_{10}|error| for Shanks Transformation Applied to the Gregory Series')
hold on
for col=2:Ncol
    AbsError =abs(pi-S(col:M-(col-1),col)); 
plot(col:M-(col-1),log10(AbsError),'.','MarkerSize',20);
end
legend('A','S(A)','S^2(A)','S^3(A)','Location','southwest')
hold off
axis([0 20 -8 1])

%
% Example: We now apply the Shanks transformaton to the sequence of Trapezoidal evaluatons used earlier to see how well it works in this case. We see that the second form of the Shanks transformation is stabler than the first!

M=length(T);
Ncol=4;
S=zeros(M,Ncol);
S(:,1)=T(:);

for col=2:Ncol
  for row=col:M-(col-1)
    S(row,col)=(S(row+1,col-1)*S(row-1,col-1)-S(row,col-1)^2)/...
        (S(row+1,col-1)-2*S(row,col-1)+S(row-1,col-1));
    % The second formula is more stable in this case
    S(row,col)=S(row+1,col-1)-(S(row+1,col-1)-S(row,col-1))^2/...
        ((S(row+1,col-1)-S(row,col-1))-(S(row,col-1)-S(row-1,col-1)));
  end
end

% Result
S
I
figure
hold on
for col=1:Ncol
    AbsError =abs(I-S(col:M-(col-1),col)); 
plot(col:M-(col-1),log10(AbsError),'.','MarkerSize',20);
end
title('Log_{10}|error| for Shanks Transformation Applied to Trapezoidal Rule')
legend('A','S(A)','S^2(A)','S^3(A)','Location','northeast')
hold off

% We see that one application of the Shanks transformation works fairly well similar to six Richardson extrapolations.
%% The End
%

function [T,t,y] = trapezoid(f,a,b,n)
%TRAPEZOID   Trapezoid formula for numerical integration.
% Input:
%   f     integrand (function)
%   a,b   interval of integration (scalars)
%   n     number of interval divisions
% Output:
%   T     approximation to the integral of f over (a,b)
%   t     vector of nodes used
%   y     vector of function values at nodes

h = (b-a)/n;
t = a + h*(0:n)';
y = f(t);
T = h * ( sum(y(2:n)) + 0.5*(y(1) + y(n+1)) ); 
end