Chapter 1

```{card}
The number of digits displayed is controlled by `format`, but the underlying values are not affected by it.

format long
p = 22/7

The absolute and relative accuracies of the approximation are as follows.

abs_accuracy = abs(p - pi)
rel_accuracy = abs(p - pi) / pi

format short
accurate_digits = -log10(rel_accuracy)

In MATLAB, values are double-precision floats unless declared otherwise.

fprintf('1 has type: %s', class(1))
fprintf('1.0 has type: %s', class(1.0))

eps

Because double precision allocates 52 bits to the significand, the default value of machine epsilon is 2^-52.

log2(eps)

The spacing between adjacent floating-point values is proportional to the magnitude of the value itself. This is how relative precision is kept roughly constant throughout the range of values. You can get the adjusted spacing by calling eps with a value.

eps(1.618)

eps(161.8)

x = 161.8 + 0.1*eps(161.8);
x - 161.8

A common mistake is to think that $\epsilon_\text{mach}$ is the smallest floating-point number. It’s only the smallest relative to 1. The correct perspective is that the scaling of values is limited by the exponent, not the mantissa. The actual range of positive values in double precision is

format short e
[realmin, realmax]

There is no double-precision number between 1 and $1+\epsilon_\text{mach}$. Thus the following difference is zero despite its appearance.

( 1 + eps / 2 ) - 1

However, the spacing between floats in $[1/2,1)$ is $\macheps/2$, so both $1-\macheps/2$ and its negative are represented exactly:

1 - (1 - eps / 2)

This is now the expected result. But we have found a rather shocking breakdown of the associative law of addition!

The polynomial $p(x) = \frac{1}{3}(x-1)(x-1-\epsilon)$ has roots 1 and $1+\epsilon$. For small values of ε, the roots are ill-conditioned.

ep = 1e-6;
a = 1/3;             % coefficients of p...
b = (-2 - ep) / 3;   % ...
c = (1 + ep) / 3;    % ...in ascending order

Here are the roots as computed by the quadratic formula.

d = sqrt(b^2 - 4*a*c);
format long   % show all digits
r1 = (-b - d) / (2*a)
r2 = (-b + d) / (2*a)

format short e
err = abs(r1 - 1) ./ abs(1)
err = abs(r2 - (1 + ep)) ./ abs(1 + ep)

The condition number of each root is

Thus, relative error in the data at the level of roundoff can grow in the result to be roughly

eps / ep

This matches the observation pretty well.

Here we show how to use Function 1.3.1 to evaluate a polynomial. Let us define a vector of the coefficients of $p(x)=(x-1)^3=x^3-3x^2+3x-1$, in ascending degree order.

c = [1, -3, 3, 1]

Now we evaluate $p(1.6)$ using the function horner.

horner(c, 1.6)

The result above is the value of $p(1.6)$.

a = 1;  b = -(1e6 + 1e-6);  c = 1;
x1 = (-b + sqrt(b^2 - 4*a*c)) / (2*a)
x2 = (-b - sqrt(b^2 - 4*a*c)) / (2*a)

The first value is correct to all stored digits, but the second has fewer than six accurate digits:

error = abs(x2 - 1e-6) / 1e-6 
accurate_digits = -log10(error)

The instability is easily explained. Since $a=c=1$, we treat them as exact numbers. First, we compute the condition numbers with respect to $b$ for each elementary step in finding the “good” root:

Using (1.2.9), the chain rule for condition numbers, the conditioning of the entire chain is the product of the individual steps, so there is essentially no growth of relative error here. However, if we use the quadratic formula for the “bad” root, the next-to-last step becomes $u_4=(-u_3) - b$, and now $\kappa=|u_3|/|u_4|\approx 5\times 10^{11}$. So we can expect to lose 11 digits of accuracy, which is what we observed. The key issue is the subtractive cancellation in this one step.

We repeat the rootfinding experiment of Demo 1.4.1 with an alternative algorithm.

a = 1;  b = -(1e6 + 1e-6);  c = 1;

First, we find the “good” root using the quadratic formula.

x1 = (-b + sqrt(b^2 - 4*a*c)) / (2*a);

Then we use the identity $x_1x_2=\frac{c}{a}$ to compute the smaller root:

x2 = c / (a*x1)

This matches the exact root to the displayed digits; to be sure we have an accurate result, we compute its relative error.

abs(x2 - 1e-6) / 1e-6

r = [-2 ,-1, 1, 1, 3, 6]';
p = poly(r)

Now we use a standard numerical method for finding those roots, pretending that we don’t know them already. This corresponds to $\tilde{y}$ in Definition 1.4.1.

rr = sort(roots(p))

disp("Root errors:") 
abs(r - rr) ./ r

It seems that the forward error is acceptably close to machine epsilon for double precision in all cases except the double root at $x=1$. This is not a surprise, though, given the poor conditioning at such roots.

Let’s consider the backward error. The data in the rootfinding problem is the polynomial coefficients. We can apply poly to find the coefficients of the polynomial (that is, the data) whose roots were actually computed by the numerical algorithm. This corresponds to $\tilde{x}$ in Definition 1.4.1.

pp = poly(rr)

We find that in a relative sense, these coefficients are very close to those of the original, exact polynomial:

disp("Coefficient errors:") 
abs(p - pp) ./ abs(p)

In summary, even though there are some computed roots relatively far from their correct values, they are nevertheless the roots of a polynomial that is very close to the original.