Page 43 - Applied Numerical Methods Using MATLAB

32 MATLAB USAGE AND COMPUTATIONAL ERRORS
For instance, consider the following two formulas:

√
√ √ √ x
f 1 (x) = x( x + 1 − x), f 2 (x) = √ √ (1.2.6)
x + 1 + x
These are theoretically equivalent, hence we expect them to give exactly the
same value. However, running the MATLAB program “nm122.m” to compute
the values of the two formulas, we see a surprising result that, as x increases,
the step of f 1 (x) incoherently moves hither and thither, while f 2 (x) approaches
1/2 at a steady pace. We might feel betrayed by the computer and have a doubt
about its reliability. Why does such a flustering thing happen with f 1 (x)?It is
because the number of significant bits abruptly decreases when the subtraction
√ √
( x + 1 − x) is performed for large values of x, which is called ‘loss of
15
significance’. In order to take a close look at this phenomenon, let x = 10 .
Then we have
√ 7
x + 1 = 3.162277660168381 × 10 = 31622776.60168381
√ 7
x = 3.162277660168379 × 10 = 31622776.60168379

These two numbers have 52 significant bits, or equivalently 16 significant digits
15
−8
8
(2 52 ≈ 10 52×3/10 ≈ 10 ) so that their significant digits range from 10 to 10 .
Accordingly, the least significant digit of their sum and difference is also the
−8
eighth digit after the decimal point (10 ).
√ √
x + 1 + x = 63245553.20336761
√ √
x + 1 − x = 0.00000001862645149230957 ≈ 0.00000002
Note that the number of significant digits of the difference decreased to 1 from
16. Could you imagine that a single subtraction may kill most of the significant
digits? This is the very ‘loss of significance’, which is often called ‘catastrophic
cancellation’.

%nm122
clear
f1 = inline(’sqrt(x)*(sqrt(x + 1) - sqrt(x))’,’x’);
f2 = inline(’sqrt(x)./(sqrt(x + 1) + sqrt(x))’,’x’);
x=1;
format long e
for k = 1:15
fprintf(’At x=%15.0f, f1(x)=%20.18f, f2(x) = %20.18f’, x,f1(x),f2(x));
x = 10*x;
end
sx1 = sqrt(x+1); sx = sqrt(x); d = sx1 - sx; s = sx1 + sx;
fprintf(’sqrt(x+1) = %25.13f, sqrt(x) = %25.13f ’,sx1,sx);
fprintf(’ diff = %25.23f, sum = %25.23f ’,d,s);