 |
|
This worksheet is designed to accompany Chapter 6 of Introduction to Scientific Programming: Computational Problem Solving Using Maple and C by Joseph L. Zachary. In it, we will use Maple to learn how to solve the quadratic equation in floating-point arithmetic without loss of significance errors. (AW, Jan 97)
You are no doubt familiar with the formula for solving quadratic equations. If you've forgotten what it is, Maple can help.
> solve(a*x^2 + b*x + c, x);
Let's write a pair of functions that compute these two roots. One will correspond to the first root (containing the + sign), and the other to the second root (containing the - sign).
> rootplus := (a,b,c) -> (-b + sqrt(b^2-4*a*c))/(2*a);
> rootminus := (a,b,c) -> (-b - sqrt(b^2-4*a*c))/(2*a);
As a check on these functions, let's try to solve x^2 + 5x + 6 = 0. By factoring the left-hand side, it is easy to see that the two roots are -3 and -2. Let's confirm this with our newly created functions.
> xplus := rootplus(1,5,6);
> xminus := rootminus(1,5,6);
It is always good practice to test functions on simple problems with known answers before using them.
The problem with the quadratic formula is that, when using floating-point arithmetic, it may be subject to loss of significance errors of the kind discussed in Section 5.3.1 of the book. We can see this with an example.
Suppose we wish to compute both roots of the equation
using floating-point numbers with ten-digit mantissas.
> xplus := rootplus(1., 2000., -3.);
> xminus := rootminus(1., 2000., -3.);
One way of checking these answers is to plug them back into the equation. We start with xminus:
> xminus^2 + 2000*xminus - 3;
Quite good. Now for xplus:
> xplus^2 + 2000*xplus - 3;
Whoops. This is not so good. It appears that xplus has not been computed accurately at all.
Here is another way to check the accuracy of the two roots. Repeat the above calculations of xplus and xminus, but user a larger value for Digits (say 20).
> Digits := 20;
Compare these higher precision values to your original answers. What do you discover? Be sure to change Digits back to 10 before continuing
> Digits := 10;
We have learned that whereas xminus has been computed very accurately, xplus was not computed accurately at all. We lost quite a number of digits in the computation of xplus. Where did the loss of significant digits occur? To find out, let us repeat the computations using rational arithmetic.
> rootplus(1, 2000, -3);
> evalf(sqrt(1000003));
> -1000+evalf(sqrt(1000003));
It should now be clear where the loss of significant digits occurred.
The number
rounds to 1000.001500. When this gets added to the value -1000,
about six digits are cancelled. The places of these digits are taken
by meaningless zeros at the end (which Maple chooses not to display).
This leaves us with inadequate precision in the final answer.
Note that the same cancellation of digits does not take place in the computation of xminus, which we can easily verify.
> rootminus(1, 2000, -3);
> evalf(sqrt(1000003));
> -1000-evalf(sqrt(1000003));
The two numbers that we are combining in the expression above are of the same sign, so the leading digits are not cancelled. If b is negative, then xplus may suffer from cancellation error. If b is positive, then xminus may suffer from cancellation error. This provides a clue to the remedy of the problem, which we now discuss.
The question is how to avoid the loss of significance error. First we do a bit of math. It is easy to show, using a piece of paper or Maple, that the product of xminus and xplus is always c/a:
> expand(rootplus(a,b,c) * rootminus(a,b,c));
If xminus suffers from cancellation errors, xplus should be free of any such errors. In that case we could compute xminus as
xminus = (c/a) / xplus.
Likewise, if xplus suffers from cancellation errors, xminus should be ok, in which case we compute
xplus = (c/a) / xminus.
Let us test this strategy by re-computing the root above which suffered from cancellation errors. (Keep in mind that a = 1, b = 2000, c = -3).
> xplusnew := (-3/1)/rootminus(1.,2000.,-3.);
Let's test this new root by plugging it back into the equation
> xplusnew^2 + 2000*xplusnew - 3;
This time, the root has been computed to full precision.
Here is another example, in which we solve
. Before doing any Maple computations, predict whether it is xplus
or xminus that will suffer from cancellation errors. Do this by
examining the signs of the constants.
Now we can do the Maple computations.
> xplus := rootplus(1., -1634., 2.);
> xplus^2 - 1634*xplus + 2;
> xminus := rootminus(1., -1634., 2.);
> xminus^2 - 1634*xminus + 2;
> xminusnew := (2/1)/xplus;
> xminusnew^2 - 1634*xminusnew + 2;
In each case, compute both roots of the equation using the standard quadratic formula. Decide which root suffers from cancellation errors by both (a) looking at the signs of the constants and (b) plugging the computed roots back into the equation. Then compute improved values of the poorly computed roots by following the technique suggested in this worksheet.
(a)
(b)
