%% Chapter 1
% _This is a MATLAB Live Editor version of the *Fundamentals of Numerical Computation* 
% Jupyter Notebooks found at_ <https://tobydriscoll.net/project/fnc/ _https://tobydriscoll.net/project/fnc/_>_._ 
% Example 1.1.2
% Recall the grade-school approximation to the number 𝜋.

p = 22/7
%%
format long
p
%% 
% Note that not all the digits displayed for |p| are the same as for 𝜋. As 
% an approximation, its absolute and relative accuracy are

abs(p-pi)            % absolute accuracy
abs(p-pi)/pi         % relative accuracy
-log10(abs(p-pi)/pi)  % accurate digits
% Example 1.1.3
% There is no double precision number between 1 and $1+\epsilon_\mbox{mach}$. 
% Thus the following difference is zero despite its appearance.

e = eps/2
(1.0 + e) - 1.0
%% 
% However, $1+\frac{\epsilon_\mbox{mach}}2$  is a double precision number, so 
% it and its negative are represented exactly:

1.0 + (e - 1.0)
%% 
% This is now the "correct" result. But we have found a rather shocking breakdown 
% of the associative law of addition!
% Example 1.3.2
% Here we show how to use |horner| to evaluate a polynomial. We first define 
% a vector of the coefficients of $p(x)=(x-1)^3=x^3-3x^2+3x-1,$ in descending 
% degree order. [Need hornereval.m im the current directory, rewritten so as not 
% to rewrite MATLAB's version.]

c = [1,-3,3,-1]
%help horner        % Uncomment to see MATLAB function
hornereval(c,1.6)
%% 
% The above is the value of 𝑝(1.6)p(1.6), up to a rounding error.
% Example 1.3.3
% Our first step is to construct a polynomial with six known roots.

r = [-2.0,-1,1,1,3,6]
p = poly(r)
%% 
% Now we use a standard numerical method for finding those roots, pretending 
% that we don't know them already. We also sort them from lowest to highest values.

r_computed = sort(roots(p))
%% 
% Here are the relative errors in each of the computed roots. [Since r_computed 
% is a column vector, we use its transpose.]

abs(r - r_computed')./r 
%% 
% It seems that the forward error is acceptably close to machine epsilon for 
% double precision in all cases except the double root at $x=1$. This is not a 
% surprise, though, given the poor conditioning at such roots.
% 
% Let's consider the backward error. The data in the rootfinding problem are 
% the polynomial coefficients. We can apply poly to find the coefficients of the 
% polynomial (that is, the data) whose roots were actually computed by the numerical 
% algorithm.

p_computed = poly(r_computed)
%% 
% We find that in a relative sense, these coefficients are very close to those 
% of the original, exact polynomial:

abs(p-p_computed)/p
%% 
% In summary, even though there are some computed roots relatively far from 
% their correct values, they are nevertheless the roots of a polynomial that is 
% very close to the original.
% Example 1.3.4

a = 1.0;
b = -(1e6+1e-6);
c = 1.0;
x1 = (-b + sqrt(b*b-4*a*c)) / (2*a)
%% 
% So far, so good. But:

x2 = (-b - sqrt(b*b-4*a*c)) / (2*a)
%% 
% The first value is correct to all stored digits, but the second has fewer 
% than six accurate digits:

-log10(abs(1e-6-x2)/1e-6 ) 
% Example 1.3.5

a = 1.0
b = -(1e6+1e-6)
c = 1.0
%% 
% First, we find the "good" root using the quadratic forumla.

x1 = (-b + sqrt(b*b-4*a*c)) / (2*a)
%% 
% Then we use the alternative formula for computing the other root.

x2 = c/(a*x1)
x2 - 1e-6