% catastrophic cancellation demos. Lec 4. Barnett 4/4/14

format long g

% EG1:
y = (1+1e-10)-1

y_ex = 1e-10;   % exact ans (at least, to within emach relative error)

abs(y - y_ex)/abs(y)   % defn of relative error in ans (see lec 1)

% note it is close to the 2e-6 we gave as upper bound. It is not close to
% emach relative error!


% EG2:
f = @(x) (1 - cos(x))./x.^2;     % want to evaluate this f for small x
x = 1e-3:1e-3:1;figure; plot(x,f(x),'-');         % looks ok, but...
x = 1e-7:1e-7:1e-4; figure; plot(x,f(x),'-');      % terrible!
x = 1e-10:1e-10:1e-7; figure; plot(x,f(x),'-');      % crazy sawtooth wrong

g = @(x) 2*(sin(x/2)./x).^2;    % reformulate via math identity
figure; plot(x,g(x),'-');  % now relatively accurate to e_mach.





% EG3:

% quadratic eqn: one small, one large root
x1 = 1e-6; x2 = 1e6;  % known roots
a = 1; b = x1+x2; c = x1*x2;   % make the coeffs of quadratic eqn
xp = (b + sqrt(b^2-4*a*c))/2*a
xm = (b - sqrt(b^2-4*a*c))/2*a % only 5 digits rel acc!
xm = c/(a*xp)  % perfect!




% quadratic eqn: nearby roots
d = 1e-6; x1 = 1-d; x2 = 1+d;  % known roots
a = 1; b = x1+x2; c = x1*x2;   % make the coeffs of quadratic eqn
xp = (b + sqrt(b^2-4*a*c))/2*a
xm = (b - sqrt(b^2-4*a*c))/2*a % only 10 digits rel acc
% problem: given a,b,c as floating-point #s we cannot do better than this!
