% This concerns cancellation errors. And is related to the HW problem
% Chapter 1, problem 27.
%
% Natasha asked, "How can you tell if you might have cancellation errors
% when you don't know exactly what the correct answer is? I.e., if you
% were trying to calculate cos(2pi) and you didn't already know that
% cos(2pi) = 1..."
%
% The answer is, "Well, you have a program that takes x and gives out y.
% You want to know what happens if you perturb x slightly and then give
% it to the program. If small perturbations are magnified hugely then
% you may have cancellation errors. In any case, if small perturbations
% are magnified hugely then you definitely have an ill-conditioned
% problem (this is how you define "ill-conditioned") and you need to
% find some other way of computing y."
%
% In the following, we take x to be 2pi, 4pi, 6pi, 8pi. We first
% estimate how far in the power series we should go before the tail
% becomes of the size of round-off. (eps is round-off). We then
% compute the partial sum, adding up in reverse order to minimize
% round-off error.
x = 2*pi;
%x = 4*pi;
%x = 6*pi;
%x = 8*pi;
tol = 10*eps;
test = 10*tol;
n = 1;
while test > tol,
% want to see how far out we have to go until the tail of the
% series becomes small:
test = abs(x^(2*n)/gamma(2*n+1));
n = n+1;
end
my_sum = 0;
for i=0:n,
% do the sum from smallest to largest:
ii = n-i;
my_sum = my_sum + (-1)^ii*x^(2*ii)/gamma(2*ii+1);
end
% okay, so that gives us an approximation of cos(x). Now we want to see
% what happens with small changes in x.
for j=1:99
xx(j) = x + (j-50)*100*eps;
my_sum(j) = 0;
for i=0:n,
% do the sum from smallest to largest:
ii = n-i;
my_sum(j) = my_sum(j) + (-1)^ii*xx(j)^(2*ii)/gamma(2*ii+1);
end
end
n
max(my_sum)-min(my_sum)
(xx(99)-xx(1))
(max(my_sum)-min(my_sum))/(xx(99)-xx(1))
% for x=2pi, n = 21, (max(my_sum)-min(my_sum))/(xx(99)-xx(1)) = 2.0e-2
% for x=4pi, n = 31, (max(my_sum)-min(my_sum))/(xx(99)-xx(1)) = 1.0e+1
% for x=6pi, n = 40, (max(my_sum)-min(my_sum))/(xx(99)-xx(1)) = 4.3e+3
% for x=8pi, n = 49, (max(my_sum)-min(my_sum))/(xx(99)-xx(1)) = 2.1e+6
%
% so for x=8pi, we have xx ranging over [8pi-1e-12,8pi+1e-12]
% and my_sum ranging over [1-2.3e-6,1+2.3e-6]. So a small perturbation
% in the input has a huge effect on the output. This is a signal that
% there might be cancellation errors. Which there are.