Contents

Demonstration of rounding error issues

Example from the book N. Higham, Accuracy and stability of Numerical Algorithms, SIAM, 2002

We are interested in evaluating the function

$$ f(x)=\frac{1 - \cos x}{x^2}. $$

x = 1.2e-8;
c = cos(x);

fprintf('x = %20.16e \n',x);
fprintf('cos(x) = %20.16e (.999999999999999928)\n',c);
fprintf('1-cos(x) = %20.16e (.7200000000000000e-16)\n',1-c);
fprintf('(1-cos(x))/x^2 = %20.16e \n',(1-c)/x^2);

x = 1.2000000000000000e-08 
cos(x) = 9.9999999999999989e-01 (.999999999999999928)
1-cos(x) = 1.1102230246251565e-16 (.7200000000000000e-16)
(1-cos(x))/x^2 = 7.7098821154524766e-01

the result is calculated inaccurately and then amplified by $x^{-2}$.

Numerically stable evaluation

Using

$$ \cos(x) = 1 - 2 \sin^2\frac{x}{2}, $$

$f(x)$ can be evaluated in a numerically stable manner

y = x/2;
f = 0.5*(sin(y)/y)^2;


fprintf('0.5*(sin(x/2)/(x/2))^2 = %20.16e \n',f);

0.5*(sin(x/2)/(x/2))^2 = 5.0000000000000000e-01


Published with MATLAB® R2019b