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 HP67 Standard Pac. A similar method which "avoids destructive cancellation" appeared in the HP15C Advanced Functions Handbook. Those "better" methods have not been consistently used in later HP 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 HP67 Standard Pac?
Avoiding Destructive Cancellation in the Solution of a Quadratic Equation
Page 31 of Maurice D. Weir's Calculus by Calculator (PrenticeHall, 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 = B^{2}  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 HP15C 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 + az^{2} . 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 ax^{2} + bx + c and a = 1E13, 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 HP 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 HP65 Standard Pac printed in 1974 or in the quadratic solution in the HP27 Owner's Handbook printed in 1976.
Pages 198 through 203 of the HP67 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 HP67 Standard Pac."
The Polynomial Evaluation program on pages 901 though 904 of the HP67 Standard Pac (1976) uses a routine which addresses destructive cancellation. It solves the test case successfully.
Pages 61 through 65 of the HP33E/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 HP41C 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 HP41 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 HP11C 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 205211 of the HP15C Advanced Functions Handbook which was published in 1982. A program for the HP15C is included which will solve the test case satisfactorily.
The QUAD option of the SOLVE menu of the HP28S (1988) does not solve the test problem correctly.
Pages 172 through 174 of my copy of the HP42S RPN Scientific Calculator Owner's Manual (10/88) present a demonstration of the conversion of a quadratic solver program from the HP41CV Owner's Manual for use on the HP42S. 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 SR56 Application Library (1976) and page IV82 of Personal Programming (the Owner's Manual for the TI58/59, 1977) only present quadratic solver programs using the textbook formula which will not solve the test case satidfactorily..
A program for the TI59 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 812 of TI PPC Notes (1991) provides a solution for the TI59 which uses the technique in Calculus by Calculator. It will yield the correct answer to the test case.
The POLY routine in the TI68 (1989), the QAD routine of the TI95 (1986) and the POLY routines in the TI85 and TI86 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 TI89 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.
hp33s Results
The following table presents results obtained with two different solution methods on the hp33s. The Textbook solution is from the Polynomial Root Finder program on pages 1520 through 1532 of the hp33s User's Guide. The Revised solution is an implementation of the technique for avoidance of destructive cancellation defined in the HP15C 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 HP15C Advanced Functions Handbook but with the first degree coefficient multiplied by minus two to change the definition of a quadratic from ax^{2}  2bx + c to ax^{2} + bx + c.
Coefficients Correct Textbook Revised
a 1E13 2E13 2E13 2E13
b 2 0.5 0 0.5
c 1
a 654323 1 1.000000921 1.000000000
b 1308644 0.9999969434 0.999996022 0.999996943
c 654321
a 11713 62.77179203 62.771792026 62.771792026
b 1470492 8.5375E05 i 62.771792026 8.5375E05 i
c 46152709
a 80841 12.2171175 12.217117552 12.217117552
b 1975288 1.374514E03 i 1.374045E03 i 1.374514E03 i
c 12066163
HP41 Results
The following table presents results obtained with three different solution methods on the HP41. The MathPac 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 HP41 program from pages 172 through 174 of the HP42S RPN Scientific Calculator Owner's Manual modified to complete the solution for complex roots. The Revised solution uses a program translated from the Revised hp33s solution in this submission.
Correct MathPac Textbook Revised
2E13 2E13 2E13 2E13
0.5 0 0 0.5
1 1.000000 0.999998472 0.999998472
0.9999969434 2.0000E05 i 0.999998472 0.999998472
62.77179203 62.77352410 62.77179203 62.77179203
8.5375E05 i 62.77006000 62.77179203 62.77179203
12.2171175 12.21711755 12.21711755 12.21711755
1.374514E03 i 1.414214E03 i 1.369104E03 i 1.377460537E03
The curious result from the MathPac 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.
An hp33s 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
K0010 x<>y
K0011 /
K0012 +/
K0013 0
K0014 1/x
K0015 STOP
K0016 GTO K
L0001 LBL L
L0002 x
L0003 x<>y
L0004 +/
L0005 2
L0006 /
L0007 ENTER
L0008 RD
L0009 x^2 x squared
L0010 
L0011 x<=0?
L0012 GTO M
L0013 SF 1
L0014 square root
L0015 x<>y
L0016 /
L0017 LASTx
L0018 R^ roll up
L0019 x<>y
L0020 /
L0021 STOP
L0022 GTO K
M0001 LBL M
M0002 +/
M0003 square root
M0004 x<>y
M0005 RD
M0006 x<>y
M0007 SGN
M0008 x
M0009 +
M0010 ENTER
M0011 R^
M0012 /
M0013 x<>yt
M0014 x=0?
M0015 STOP
M0016 x=0?
M0017 GTO K
M0018 RCL C
M0019 x<>y
M0020 /
M0021 STOP
M0022 GTO K
