(15C) Halley's Method

[Thomas Klemm]

Halley's Method

References

• Halley's Irrational Formula
• Halley-Verfahren
• Halley's method
• HP Article VA031 - Boldly Going - Going Back to the Roots
• 35s - find roots of 3rd and higher order equations
• Solver for the HP-15C
  

Formula

We start with Halley's Irrational Formula:

\(C_f(x)=x-\frac{2f(x)}{f{'}(x)+\sqrt{[f{'}(x)]^2-2f(x)f{''}(x)}}\)

However, we reduce the fraction with \(f{'}(x)\) to get:

\(C_f(x)=x-\frac{\frac{2f(x)}{f{'}(x)}}{1 + \sqrt{1 - \frac{2f(x)f{''}(x)}{[f{'}(x)]^2}}}\)

This allows us to reuse the following expression:

\(\frac{2f(x)}{f{'}(x)}\)

Also the expression avoids cancellation if \(f(x) \to 0\).

This is not the case for the alternative expression:

\(C_f(x)=x-\frac{1 - \sqrt{1 - \frac{2f(x)f{''}(x)}{[f{'}(x)]^2}}}{\frac{f{''}(x)}{f{'}(x)}}\)

Definitions

These definitions are used:

• \(z = x + i y\)
• \(f^{+} = f(z + h) = f^{+}_x + i f^{+}_y\)
• \(f^{-} = f(z - h) = f^{-}_x + i f^{-}_y\)
• \(f = f(z) = f_x + i f_y\)
• \(f{'} = f{'}(z) h \approx \frac{f^{+} - f^{-}}{2}\)
• \(f{''} = f{''}(z) h^2 \approx f^{+} - 2f + f^{-}\)
  

Registers

Intermediate results are kept in these registers:

\(\begin{matrix}
R_0 & h \\
R_1 & x \\
R_2 & y \\
R_3 & f^{-}_x \\
R_4 & f^{-}_y \\
R_5 & f_x \\
R_6 & f_y \\
I & \text{program} \\
\end{matrix}\)

Description

These are the steps of program A:

002 - 004: store \(z\)
005 - 062: calculate the next approximation
006 - 011: \(f(z - h)\)
012 - 016: \(f(z)\)
017 - 018: \(f(z + h)\)
019 - 029: \(f{'}\)
030 - 036: \(f{''}\)
037 - 040: \(2f \div f{'}\)
041 - 043: \(f{''} \div f{'}\)
045 - 053: \(dz\)
054 - 056: \(z = z - dz\)
057 - 062: loop until it is good enough
063 - 066: return \(z\)
067 - 075: calculate \(f(z + w) | w \in \{-h, 0, h\}\)

The programs B - E are examples from Valentin's article.

Example

The step-size h is stored in register 0, while the name/number of the program is specified in the variable I.
So let's assume we want to solve program B with step-size h = 10^-4 and starting guess 2:

RCL MATRIX B
STO I

EEX 4 CHS
STO 0

2 A

Intermediate values of |dz| are displayed:

(( running ))
0.148589460
(( running ))
0.002695411
(( running ))
0.000000017
(( running ))
0.000000000
(( running ))
1.854105968

Should that annoy you just remove the PSE-command in line 058.

Programs

A: Halley's Method

Code:

   001 {    42 21 11 } f LBL A
   002 {       44  1 } STO 1
   003 {       42 30 } f Re↔Im
   004 {       44  2 } STO 2
   005 {    42 21  0 } f LBL 0
   006 {       45  0 } RCL 0
   007 {          16 } CHS
   008 {       32  1 } GSB 1
   009 {       44  3 } STO 3
   010 {       42 30 } f Re↔Im
   011 {       44  4 } STO 4
   012 {           0 } 0
   013 {       32  1 } GSB 1
   014 {       44  5 } STO 5
   015 {       42 30 } f Re↔Im
   016 {       44  6 } STO 6
   017 {       45  0 } RCL 0
   018 {       32  1 } GSB 1
   019 {          36 } ENTER
   020 {          36 } ENTER
   021 {       45  3 } RCL 3
   022 {       45  4 } RCL 4
   023 {       42 25 } f I
   024 {          40 } +
   025 {          34 } x↔y
   026 {       43 36 } g LSTx
   027 {          30 } −
   028 {           2 } 2
   029 {          10 } ÷
   030 {          34 } x↔y
   031 {       45  5 } RCL 5
   032 {       45  6 } RCL 6
   033 {       42 25 } f I
   034 {          36 } ENTER
   035 {          40 } +
   036 {          30 } −
   037 {           1 } 1
   038 {       43 36 } g LSTx
   039 {       43 33 } g R⬆
   040 {          10 } ÷
   041 {       43 33 } g R⬆
   042 {       43 36 } g LSTx
   043 {          10 } ÷
   044 {          34 } x↔y
   045 {          20 } ×
   046 {       43 36 } g LSTx
   047 {          33 } R⬇
   048 {          30 } −
   049 {          11 } √x̅
   050 {          40 } +
   051 {          10 } ÷
   052 {       45  0 } RCL 0
   053 {          20 } ×
   054 {    44 30  1 } STO − 1
   055 {       42 30 } f Re↔Im
   056 {    44 30  2 } STO − 2
   057 {       43 16 } g ABS
   058 {       42 31 } f PSE
   059 {       45  0 } RCL 0
   060 {       43 11 } g x²
   061 {    43 30  8 } g TEST x<y
   062 {       22  0 } GTO 0
   063 {       45  1 } RCL 1
   064 {       45  2 } RCL 2
   065 {       42 25 } f I
   066 {       43 32 } g RTN
   067 {    42 21  1 } f LBL 1
   068 {       45  1 } RCL 1
   069 {       45  2 } RCL 2
   070 {       42 25 } f I
   071 {          40 } +
   072 {          36 } ENTER
   073 {          36 } ENTER
   074 {          36 } ENTER
   075 {       22 25 } GTO I

B: Find a root of : \(x^x = \pi\)

Code:

   076 {    42 21 12 } f LBL B
   077 {          14 } yˣ
   078 {       43 26 } g π
   079 {          30 } −
   080 {       43 32 } g RTN

C: Find all roots of: \(( 2 + 3i ) x^3 - (1 + 2i ) x^2 - ( 3 + 4i ) x - ( 6 + 8i ) = 0\)

Code:

   081 {    42 21 13 } f LBL C
   082 {           2 } 2
   083 {          36 } ENTER
   084 {           3 } 3
   085 {       42 25 } f I
   086 {          20 } ×
   087 {           1 } 1
   088 {          36 } ENTER
   089 {           2 } 2
   090 {       42 25 } f I
   091 {          30 } −
   092 {          20 } ×
   093 {           3 } 3
   094 {          36 } ENTER
   095 {           4 } 4
   096 {       42 25 } f I
   097 {          30 } −
   098 {          20 } ×
   099 {           6 } 6
   100 {          36 } ENTER
   101 {           8 } 8
   102 {       42 25 } f I
   103 {          30 } −
   104 {       43 32 } g RTN

D: Attempt to find a complex root of: \(x^3 - 6x - 2 = 0\)

Code:

   105 {    42 21 14 } f LBL D
   106 {       43 11 } g x²
   107 {           6 } 6
   108 {          30 } −
   109 {          20 } ×
   110 {           2 } 2
   111 {          30 } −
   112 {       43 32 } g RTN

E: Solve Leonardo di Pisa's equation: \(x^3 + 2 x^2 + 10 x - 20 = 0\)

Code:

   113 {    42 21 15 } f LBL E
   114 {           2 } 2
   115 {          40 } +
   116 {          20 } ×
   117 {           1 } 1
   118 {           0 } 0
   119 {          40 } +
   120 {          20 } ×
   121 {           2 } 2
   122 {           0 } 0
   123 {          30 } −
   124 {       43 32 } g RTN