%% Class 3: unavoidable error, numerically stable and numerically unstable algorithms
%% Introduction
% We want to compute f(x) for a given input value x.
% Assume that there are no measurement errors. In the ideal algorithm
% the given value is rounded to the closest machine number: \hat x = fl(x)
% then use \hat x and compute \tilde y:= f(\hat x)  "exactly" (with some extra precision)
% then round this extra precision result \tilde y to the closest machine number: \hat y=fl(\tilde y)
% As a result we obtain the unavoidable error
% |\epsilon_{\hat y}| \le |c_f(x)| \epsilon_M +\epsilon_M
% where c_f(x) = x\cdot f'(x)/f(x) denotes the condition number of f and \epsilon_M denotes the machine epsilon (which is 10^{-16} for double precision numbers which are used in Matlab).
% In a practical computation we break up the computation of f(x) into a sequence of operations which are available on our machine ("algorithm"). Then we run the algorithm on our computer using machine numbers and machine arithmetic.
%% Example
% We want to compute y=1-\cos(x) for x=5\cdot10^{-4}.
%% Algorithm 1
% The obvious way to compute y is

format compact; format long g
x = 5e-4;
yalg1 = 1-cos(x)     

% We call this Algorithm 1: first compute y_1:=\cos(x) and then y:=1-y_1.
% Matlab uses machine numbers with 16 digits accuracy and displays 15 digits.
% How accurate is this computed result?
% Is Algorithm 1 numerically stable, i.e., "optimal"?
% Can we find an algorithm which gives a much smaller error?
%% Finding the relative error
% In order to find the relative error we need the "exact result". One way is to use higher precision: If the original computation was in single precision we can get a more accurate result using double precision.
% Matlab uses double precision by default. One way to obtain a higher precision result is to use the symbolic toolbox and vpa which uses 30 digits accuracy by default:

X = sym(5*10^-4);     % symbolic input
Y = vpa(1-cos(X))     % evaluate using 30 digit accuracy
relerr_alg1 = double(yalg1-Y)/yex   % double() converts to machine number

% We see that the realtive error is about 3.4\cdot 10^{-10}.
%% Unavoidable error
% Here f(x)=1-\cos(x), so the condition number is
% c_f(x)=\frac{x\cdot f'(x)}{f(x)}=\frac{x\cdot \sin x}{1-\cos x}

cf = x*sin(x)/(1-cos(x))

% giving the unavoidable error |c_f(x)|\epsilon_M+\epsilon_M

epsM = 1e-16;
unav_err = cf*epsM + epsM

% The error of Algorithm 1 is 1.4\cdot10^{-10} which is much larger than the unavoidable error 3\cdot 10^{-16}. Therefore Algorithm 1 is numerically unstable.
%% Why do we get such a large error for algorithm 1?
% Note that we first compute y_1 :=\cos(x). The computed result is rounded to the closest machine number, this results in an error which can be as large as \epsilon_M.
% Then we compute y:=1-y_1. Note that y_1 \approx1, so we have subtractive cancellation, and the relative error present in y_1  will be multiplied by \left|\frac{y_1}{1-y_1}\right|:

y1 = cos(x);
factor = abs(y1/(1-y1))
factor*epsM

% This means that the roundoff error in y_1 may cause an error as large as 8\cdot 10^{-10} in y. This is the order of magnitude of the error we observed for algorithm 1.
%% Algorithm 2
% How can we compute a more accurate result in machine arithmetic? The idea is to use mathematical identities to rewrite the expression 1-\cos x so that the new expression does not suffer from subtractive cancellation.
% A useful trick is to multiply and divide by 1+\cos x:
% y=1-\cos x=\frac{(1-\cos x)(1+\cos x)}{1+\cos x}=\frac{1-\cos^2 x}{1+\cos x}
% So far this does not help: Evaluating 1-\cos^2 x still causes subtractive cancellation. But using \sin^2 x+\cos^2 x=1 we obtain
% y=1-\cos x=\frac{\sin^2 x}{1+\cos x}
% The new expression is now free from subtractive cancellation. This gives Algorithm 2:

yalg2 = sin(x)^2/(1+cos(x))
relerr_alg2 = double(yalg2-Y)/yex

% Note that this is of the order of the unavoidable error or less. Hence Algorithm 2 is numerically stable.