% Math 5610/6860 - Numerical Analysis I
% Computer lab #2, problem 1

% number of terms to compute in the sequence
n = 8;
% a vector that will contain all computed terms in the sequence
x = zeros(n,1);

% compute the terms in the sequence
x(1)=2; x(2)=3;
for k=2:n-1,
    x(k+1) = x(k) - (x(k)^2-2)/(x(k)+x(k-1));
end;

% compute absolute error
abserr = abs(x-sqrt(2));


% plot the absolute error
%
% * you can try to do this with n=20 and you will see that for n>8 the iterates
%   are within machine precision of the limit so the convergence is very fast
%
% * I am using semilogy which means the scale in the y direction is logarithmic,
%   this shows how small the errors are more clearly.
figure(1);
semilogy(abserr);
title('error at step n');
xlabel('n');

% plot of the log errors at n (x axis) vs log errors at n+1 (y axis)
%
% first notice that if x_n converges to xstar with a order alpha convergence
% rate then:
%
%  | x_{n+1} - xstar | <= C | x_n - xstar |^alpha
%
% taking logs:
%
%  log | x_{n+1} - xstar | <= log(C)  + alpha * log | x_n - xstar |
%
% which explains why the log log plot of the errors at n vs the errors at n+1
% is a line with slope the convergence rate.
%
% Using a loglog plot (the matlab command for plotting with a logarithmic scale
% on both the x and y directions) is not recommended here if you want to fit
% the line with matlab's built-in linear fitting routines (through the
% Tools->Basic Fitting menu in the figure). If you this Matlab will try to fit
% the data as it is, and will fail miserably. It is better to compute the logs
% manually  and then do the data fit as is done below. The convergence rate you
% get should be around 1.6.
figure(2);
plot(log(abserr(1:n-1)),log(abserr(2:n)));


% A crude way of estimating the convergence rate is to estimate the slope of
% the line by looking a tthe slope between any two consecutive points on the
% line. Instead of writign a for loop we can use Matlab's notation to do it:
diff(log(abserr(2:n)))./diff(log(abserr(1:n-1)))
% the diff command is very useful when you want to compute the consecutive
% differences of a vector, see help diff for more info on what it does.


% Here is how you would do the linearfit by hand. This was not covered during
% the computer lab, but I can explain it to anyone interested (you should have
% seen something like this in your linear algebra class, I am doing the fit by
% solving a linear least squares problem)

% form linear system
A = [ ones(n-1,1), log(abserr(1:n-1)) ];
b = log(abserr(2:n));

y = (A'*A)\A'*b; % solve normal equations
% y(1) == intercept, y(2) == slope
fprintf('the convergence rate is about %f\n',y(2));

% superimpose the line on the plot
hold on; % allows to add more curves to a plot without erasing previous plots
plot(log(abserr(1:n-1)),y(2)*log(abserr(1:n-1)) + y(1),'r');

% add a legend
legend('error','convergence rate model');
xlabel('log absolute error at n'); ylabel('log absolute error at n+1')

% Finally if we did not know that the sequence converges to sqrt(2) we can
% simply replace the true value by an iterate with larger n (if the computation
% is very expensive, for example if we are solving a PDE on a domain, sometimes
% we just take n+1). The same convergence trends should be observable with this
% ansatz.

% One final note: We could have done this with any base log. In matlab:
%  log(x)   corresponds to ln(x)    or the natural logarithm (base e)
%  log10(x) corresponds to log10(x) or base 10 logarithm
