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

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!

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


% 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!
