{VERSION 2 3 "IBM INTEL NT" "2.3" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 256 "" 0 1 0 0 15 0 1 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 6 6 0 0 0 0 0 0 -1 0 }{PSTYLE "T itle" 0 18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }} {SECT 0 {PARA 18 "" 0 "" {TEXT -1 42 "Case Study: Solving the Quadrati c Equation" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 53 "This worksheet is designed to accompany Chapter 6 of " }{TEXT 256 87 "Introduction to Scientific Programming: Computational Problem \+ Solving Using Maple and C" }{TEXT -1 175 " by Joseph L. Zachary. In i t, we will use Maple to learn how to solve the quadratic equation in f loating-point arithmetic without loss of significance errors. (AW, \+ Jan 97)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 34 "Implementing the Quadratic Formula" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 124 "You are no doubt familiar with the formula for solv ing quadratic equations. If you've forgotten what it is, Maple can he lp." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "solve(a*x^2 + b*x + c, x);" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {EXCHG {PARA 0 "" 0 "" {TEXT -1 184 "Let's write a pair of functions t hat compute these two roots. One will correspond to the first root ( containing the + sign), and the other to the second root (containing t he - sign)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "rootplus := (a,b,c) -> (-b + sqrt(b^2-4*a*c))/(2*a); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "rootminus := (a,b,c) -> (-b - sqrt(b^2-4*a*c))/(2*a); " }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 209 "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 functi ons." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "xplus := rootplus(1,5,6);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "xminus := rootminus(1,5,6);" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 101 "It is always good practice to t est functions on simple problems with known answers before using them. " }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 20 "Loss of Significance" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 220 "Th e problem with the quadratic formula is that, when using floating-poin t 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." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 54 "Suppose we wish to compute both roots of the equation " }{XPPEDIT 18 0 "x^2 + 2000x - 3 = 0" "/,(*$%\"xG\"\"#\"\"\"*&\"%+?F'F%F'F'\"\"$! \"\"\"\"!" }{TEXT -1 55 " using floating-point numbers with ten-digit \+ mantissas." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "xplus := rootplus(1., 2000., -3.);" }}}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "xminus := \+ rootminus(1., 2000., -3.);" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 96 "One way of checking these answers is to p lug them back into the equation. We start with xminus:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "xminus^2 + 200 0*xminus - 3;" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Quite good. Now for xplus:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "xplus^2 + 2000*xplus - 3; " }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 93 "Whoops. This is not so good. It appears that xplus has not been \+ computed accurately at all." }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {EXCHG {PARA 0 "" 0 "" {TEXT -1 156 "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)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "Digits := 20;" } }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 141 "Compare these higher precision values to your original answers. What do you discover? Be sure to change Digits back to 10 before continui ng" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "Digits := 10;" }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 11 "Explanation" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 304 "We have learned that whereas xminus has been computed very accu rately, 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 computatio ns using rational arithmetic." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "rootplus(1, 2000, -3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "evalf(sqrt(1000003));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "-1000+evalf(sqrt(1000003));" }}} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 85 "It should now be clear where the loss of significant digits occurred. \+ The number " }{XPPEDIT 18 0 "sqrt(1000003" "-%%sqrtG6#\"(.++\"" } {TEXT -1 275 " rounds to 1000.001500. When this gets added to the \+ value -1000, about six digits are cancelled. The places of these dig its are taken by meaningless zeros at the end (which Maple chooses not to display). This leaves us with inadequate precision in the final answer. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 120 "Note that the same cancellation of digits does not take place in \+ the computation of xminus, which we can easily verify." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "rootminus(1, 2 000, -3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "evalf(sqrt(100 0003));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "-1000-evalf(sqrt (1000003));" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 " " {TEXT -1 331 "The two numbers that we are combining in the expressio n 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." } }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 24 "Maintaining Significance" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {EXCHG {PARA 0 "" 0 "" {TEXT -1 195 "The question is how to avoid the \+ loss of significance error. First we do a bit of math. It is easy t o show, using a piece of paper or Maple, that the product of xminus an d xplus is always c/a:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 43 "expand(rootplus(a,b,c) * rootminus(a,b,c));" } }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 127 "If xminus suffers from cancellation errors, xplus should be free of \+ any such errors. In that case we could compute xminus as" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 47 " \+ xminus = (c/a) / xplus." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 98 "Likewise, if xplus suffers from cancellat ion errors, xminus should be ok, in which case we compute" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 48 " \+ xplus = (c/a) / xminus." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 143 "Let us test this strategy by re-computin g the root above which suffered from cancellation errors. (Keep in mi nd that a = 1, b = 2000, c = -3)." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "xplusnew := (-3/1)/rootminus(1.,200 0.,-3.);" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 62 "Let's test this new root by plugging it back into the equ ation" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "xplusnew^2 + 2000*xplusnew - 3;" }}}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "This time, the root has been c omputed to full precision." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 43 "Here is another example, in which we solve " } {XPPEDIT 18 0 "x^2 - 1634x + 2" ",(*$%\"xG\"\"#\"\"\"*&\"%M;F&F$F&!\" \"\"\"#F&" }{TEXT -1 170 ". Before doing any Maple computations, pred ict whether it is xplus or xminus that will suffer from cancellation e rrors. Do this by examining the signs of the constants." }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 37 "Now we can do the \+ Maple computations." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "xplus := rootplus(1., -1634., 2.);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "xplus^2 - 1634*xplus + 2;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "xminus := rootminus(1., -163 4., 2.);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "xminus^2 - 1634 *xminus + 2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "xminusnew : = (2/1)/xplus;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "xminusnew ^2 - 1634*xminusnew + 2;" }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 9 "Exercises" }}{PARA 0 "" 0 "" {TEXT -1 363 "In each case, compute both roots of the equation using the standard q uadratic 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 t his worksheet." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 6 "(a) " }{XPPEDIT 18 0 "2x^2 - 4000x - 1 = 0" "/,(*&\"\"# \"\"\"*$%\"xG\"\"#F&F&*&\"%+SF&F(F&!\"\"\"\"\"F,\"\"!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 6 "(b) " }{XPPEDIT 18 0 "5x^2 + 8000x + 2 = 0" "/,(*&\"\"&\"\"\"*$%\"xG\"\"#F&F&*&\"%+!)F&F(F& F&\"\"#F&\"\"!" }}}}{MARK "0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 }