|
|
This notebook is designed to accompany Chapter 6 of Introduction to Scientific Programming: Computational Problem Solving Using Mathematica and C by Joseph L. Zachary. In it, we will use Mathematica to learn how to solve the quadratic equation in floating-point arithmetic without loss of significance errors. (AW, 31Jul97)
You are no doubt familiar with the formula for solving quadratic equations. If you've forgotten what it is, Mathematica can help.
Solve[a * x^2 + b * x + c == 0, 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.
StyleBox[{xplus = rootplus[1, 5, 6], xminus = rootminus[1, 5, 6]}, ShowStringCharacters -> True]
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 x^2 + 200000x - 3 = 0 using floating-point numbers.
StyleBox[{xplus = rootplus[1., 200000., -3.], xminus = rootminus[1., 200000., -3.]}, ShowStringCharacters -> True]
One way of checking these answers is to plug them back into the equation. We start with xminus:
xminus^2 + 200000 * xminus - 3
Quite good. Now for xplus:
xplus^2 + 200000 * xplus - 3
Whoops. This is not so good. It appears that xplus has not been computed accurately at all.
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.
Simplify[rootplus[1, 200000, -3]]
InputForm[N[Sqrt[10000000003]]]
InputForm[-100000 + N[Sqrt[10000000003]]]
It should now be clear where the loss of significant digits occurred. The number Sqrt[10000000003] rounds to 100000.000015. When this gets added to the value -100000, about six digits are cancelled. The places of these digits are taken by meaningless numbers at the end. 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.
Simplify[rootminus[1, 200000, -3]]
InputForm[N[Sqrt[10000000003]]]
InputForm[-100000 - N[Sqrt[10000000003]]]
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 Mathematica, 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 = 200000, c = -3).
xplusnew = (-3/1)/rootminus[1., 200000., -3.]
Let's test this new root by plugging it back into the equation
xplusnew^2 + 200000 * xplusnew - 3
This time, the root has been computed to full precision.
Here is another example, in which we solve x^2 - 16340x + 2. Before doing any Mathematica 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 Mathematica computations.
xplus = rootplus[1., -16340., 2.]
xplus^2 - 16340 * xplus + 2
xminus = rootminus[1., -16340., 2.]
xminus^2 - 16340 * xminus + 2
xminusnew = (2/1)/xplus
xminusnew^2 - 16340 * 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 notebook.
(a) 2x^2 - 400000x - 1 = 0
(b) 5x^2 + 800000x + 2 = 0