|Which quadratic solution should we use?|
Message #1 Posted by Palmer O. Hanson, Jr. on 3 Mar 2007, 11:07 p.m.
A "more efficient and accurate" method for the solution of quadratic equations appeared in the HP-67 Standard Pac. A similar method which "avoids destructive cancellation" appeared in the HP-15C Advanced Functions Handbook. Those "better" methods have not been consistently used in later H-P applications.
Is there a reason for not using the "more efficient and accurate" methods?
Is there an earlier reference on the "more efficient and accurate" methods than that in the HP-67 Standard Pac?
Avoiding Destructive Cancellation in the Solution of a Quadratic Equation
Page 31 of Maurice D. Weir's Calculus by Calculator (Prentice-Hall, 1982) gives the typical algebra textbook solution for a quadratic equation but then states:
"... The program QUADS calculates the real roots (if they exist) for any specified quadratic equation. The formula used is a modification of Equation 1.11 [the textbook form] to reduce roundoff error."
Analysis of the QUADS program reveals that the solution sequence is:
1. Calculate the determinant D = B2 - 4AC
2. Calculate the value E = (sqr(D) + abs(B)/2A
3. Change the sign of E if B<0
4. Calculate the real roots as C/AE and E.
A more thorough treatment of the subject with a similar solution method appears in the appendix "Accuracy of Numerical Calculations" to the HP-15C Advanced Functions Handbookalso published in 1982. The discussion on page 191 presents the textbook formula and an algebraically equvalent calculation sequence which avoids "destructive cancellation". (Readers of that reference should note that for some reason the authors chose to define the quadratic as c - 2bz + az2 . The more common definition is used here.) The discussion states
"... Such a program will be listed later (page 205) and must be used in those instances, common in engineering, when the smaller root y is needed accurately despite the fact that the quadratic's other unwanted root is relatively large."
A demonstration problem is provided. For a quadratic equation of the form ax2 + bx + c and a = 1E-13, b = -2 and c = 1 the exact solution is r1 = 2E13 and r2 = 0.5. Implementation of the textbook formula on most HP and TI machines will yield r1 = 2E13 and r2 = 0. .
The rest of this submission presents all the gory details which led up to this submission.
Some H-P History on a Better Quadratic Solution Methoid
This history is limited to the hardware and documentation available from my own collection.
The method was not included in the quadratic solution in the HP-65 Standard Pac printed in 1974 or in the quadratic solution in the HP-27 Owner's Handbook printed in 1976.
Pages 198 through 203 of the HP-67 Owner's Handbook and Programming Guide (Revision C 12/76) presented a quadratic solver program which used the textbook formula. The end of the discussion includes the following statement "For a more efficient and accurate method of finding the roots of a quadratic equation, see the Polynomial Evaluation program in your HP-67 Standard Pac."
The Polynomial Evaluation program on pages 9-01 though 9-04 of the HP-67 Standard Pac (1976) uses a routine which addresses destructive cancellation. It solves the test case successfully.
Pages 61 through 65 of the HP-33E/33C Owner's Handbook and Programming Guide (April 1979) presents a quadratic solver program which uses the textbook formula. It can not be expected to and does not solve the test case successfully.
Pages 179 through 183 of my copy of the HP-41C Owner's Handbook and Programming Guide (Revision B 8/79) present a quadratic solver program which uses the textbook formula. The program cannot solve the test case satisfactorily. A note at the bottom of page 179 does offer the following caution "Some values of a, b, and c may result in misleading answers because their solutions require greater than twelve digits of accuracy."
The "Polynomial Solutions/Evaluation" program of the HP-41 Math Pac (11/79) includes a quadratic solver which uses the textbook formula. It cannot solve the test case satisfactorily.
Pages 206 through 210 of the HP-11C Owner's Handbook and Problem Solving Guide (Revision B 9/81) describes how to develop a quadratic solver using the textbook formula.
A thorough discussion of the method for avoiding destructive cancellation appears on pages 191 and 205-211 of the HP-15C Advanced Functions Handbook which was published in 1982. A program for the HP-15C is included which will solve the test case satisfactorily.
The QUAD option of the SOLVE menu of the HP-28S (1988) does not solve the test problem correctly.
Pages 172 through 174 of my copy of the HP-42S RPN Scientific Calculator Owner's Manual (10/88) present a demonstration of the conversion of a quadratic solver program from the HP-41CV Owner's Manual for use on the HP-42S. The solution uses the textbook formula and will not solve the test case satisfactorily..
Some TI History on A Better Quadratic Solution Method
The "Solution of Quadratic Equation" in the SR-56 Application Library (1976) and page IV-82 of Personal Programming (the Owner's Manual for the TI-58/59, 1977) only present quadratic solver programs using the textbook formula which will not solve the test case satidfactorily..
A program for the TI-59 by Bill Skillman which uses a method which will solve the test case satisfactorily was published in late 1978 in V3N12P2 of 52 Notes.
My article "Improved Solution for the Quadratic" on pages 8-12 of TI PPC Notes (1991) provides a solution for the TI-59 which uses the technique in Calculus by Calculator. It will yield the correct answer to the test case.
The POLY routine in the TI-68 (1989), the QAD routine of the TI-95 (1986) and the POLY routines in the TI-85 and TI-86 all yield the correct answer to the test case but there is no discussion of the use of an improved method in the manuals. The TI-89 yields inconsistent results. The Polynomial Root Finder of the Apps Desktop returns roots of 2.E13 and 0.51. However, if the test case is entered through the home screen as either a" factor" problem or as a "solve" problem a correct answer will be returned.
The following table presents results obtained with two different solution methods on the hp-33s. The Textbook solution is from the Polynomial Root Finder program on pages 15-20 through 15-32 of the hp-33s User's Guide. The Revised solution is an implementation of the technique for avoidance of destructive cancellation defined in the HP-15C Advanced Functions Handbook. A listing appears at the end of this submission. The coefficients are from test cases 3 through 6 on page 207 of the HP-15C Advanced Functions Handbook but with the first degree coefficient multiplied by minus two to change the definition of a quadratic from ax2 - 2bx + c to ax2 + bx + c.
Coefficients Correct Textbook Revised
a 1E-13 2E13 2E13 2E13
b -2 0.5 0 0.5
a 654323 1 1.000000921 1.000000000
b -1308644 0.9999969434 0.999996022 0.999996943
a 11713 62.77179203 62.771792026 62.771792026
b -1470492 8.5375E-05 i 62.771792026 8.5375E-05 i
a 80841 12.2171175 12.217117552 12.217117552
b -1975288 1.374514E-03 i 1.374045E-03 i 1.374514E-03 i
The following table presents results obtained with three different solution methods on the HP-41. The Math-Pac solution uses a version of the textbook formula but with the coefficient of the second degree term set to one. The Textbook solution uses the HP-41 program from pages 172 through 174 of the HP-42S RPN Scientific Calculator Owner's Manual modified to complete the solution for complex roots. The Revised solution uses a program translated from the Revised hp-33s solution in this submission.
Correct Math-Pac Textbook Revised
The curious result from the Math-Pac for the second problem results from the definition of the coefficient of the second degree term as unity. That means that the coefficient of the first degee term is entered as 654321/654323 and the constant term as -1308644/654323. I obtain a similar curious result if I enter those coefficients in either the textbook solution or in the revised solution.
2E13 2E13 2E13 2E13
0.5 0 0 0.5
1 1.000000 0.999998472 0.999998472
0.9999969434 2.0000E-05 i 0.999998472 0.999998472
62.77179203 62.77352410 62.77179203 62.77179203
8.5375E-05 i 62.77006000 62.77179203 62.77179203
12.2171175 12.21711755 12.21711755 12.21711755
1.374514E-03 i 1.414214E-03 i 1.369104E-03 i 1.377460537E-03
An hp-33s Program Which Avoids Destructive Cancellation
K0001 LBL K
K0002 CF 1
K0003 INPUT A
K0004 INPUT B
K0005 INPUT C
K0006 RCL A
K0007 x ne 0 ? not equal to
K0008 GTO L
K0009 RD roll down
K0016 GTO K
L0001 LBL L
L0009 x^2 x squared
L0012 GTO M
L0013 SF 1
L0014 square root
L0018 R^ roll up
L0022 GTO K
M0001 LBL M
M0003 square root
M0017 GTO K
M0018 RCL C
M0022 GTO K