(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

[Albert Chan]

(02-26-2022 11:09 AM)Thomas Klemm Wrote:  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 also avoid cancellation if Re(f'(x)) < 0, since Re(√z) ≥ 0

Comment: two C_f expressions are not the same.
First form denominator + is really ± (2 roots of quadratic)
Second form pick the smaller correction, exactly what we wanted.

[Albert Chan]

Halley's Irrational Formula may cause sqaure root of negative number problem.

(01-07-2020 06:13 PM)Albert Chan Wrote:  Example, for N=30, PV=6500, PMT=-1000, FV=50000, we have I = 8.96% or 11.10%

Does anyone knows what returns should be reported for above investments ?

Halley's rational formula for rate search work, without issue.

lua> n,pv,pmt,fv = 30,6500,-1000,50000
lua> i1, i2 = edge_i(n,pv,pmt,fv)
lua> i1, i2
-0.02     0.15384615384615385
lua> g1 = loan_rate2(n,pv,pmt,fv,i1)
lua> g2 = loan_rate2(n,pv,pmt,fv,i2)
lua> for i = 1,6 do print(i, g1(), (g2())) end
1   0.05032222218969118    0.11874137324380304
2   0.08034524496312957    0.11136968662548538
3   0.0889289402959179     0.11100337818114
4   0.0896051698287111     0.11100328640549177
5   0.08960585626123155    0.11100328640549177
6   0.08960585626123234    0.11100328640549177

Halley's irrational formula failed, even on first try, at i = -0.02
Halley's modified slope = (f' + sign(f')*sqrt((f')^2 - 2*f*f'') / 2

lua> f0 = 1356.1628502335457 -- i = -0.02
lua> f1 = -26468.743095455087
lua> f2 = 280416.32885371713
lua> f1*f1 - 2*f0*f2
-59986054.527367234

Perhaps it is because of negative guess rate ? No !
Even with guess i = 0.07, we still hit square root of negative number problem.

lua> f0 = 53.131798377782616 -- i = 0.07
lua> f1 = -4261.51788308051
lua> f2 = 171458.18285005045
lua> f1*f1 - 2*f0*f2
-59228.53500474244

[Gil]

About Albert Chan's problem and corresponding equation

6500*R^30+50000=(R^30-1)/(R-1)*1000,
with R=1+interest_rate/100,

is it possible to know for sure / to find that there exactly two "analytical" solutions, but without making a graph or trying guesses?

[Albert Chan]

Hi, Gil

We just count sign changes, see Descartes rules of signs
N+1 Cash flows: [PV, PMT, ... , PMT, PMT+FV] = [6500, -1000, ..., -1000, 49000]

2 sign changes, thus 0 or 2 positive roots (for R = 1+i)
Based on rate edges, (PMT/FV , PMT/-PV) = (-1/50, 2/13)

2 Rates (if exist at all): -1/50 < i < 2/13

If we can show we have 1 root, we must have 2.

[Gil]

Of course.
Thanks for reminding me.
Regards,
Gil

[Albert Chan]

(06-14-2022 01:20 AM)Gil Wrote:  6500*R^30+50000=(R^30-1)/(R-1)*1000,
with R=1+interest_rate/100,

Above setup were solving NFV = 0, for R

NFV = 6500*R^30 + 50000 - (R^30-1)/(R-1)*1000

For comparison, I was solving for f = NPMT/n = 0

f = NFV * (R-1)/(R^30-1) = (56500 / (R^30-1) + 6500) * (R-1) - 1000

Solving for f=0 for R is more stable, and, likely converge faster.

From the edges, solve f=0 with Newton's method, we have 1-sided convergence.
In other words, starting from 2 edges, we can get both roots (if existed).

Newton's method for NFV=0 starting from edges, for this example, we get only 1 root.
Interestingly, Halley's Irrational formula for NFV=0 work for small edge, R = 1+i = 0.98

0.98
1.0810815846
1.08936323285
1.08960584869
1.08960585626

For R < 1.12175, argument inside square root are positive.
Thus, it will not work on the other edge, with R = 1+i = 1+2/13 ≈ 1.15385

We have only very small guess window to get the other root, R ≈ 1.11100

R(where NFV'>0) to R(where (NFV')^2 - 2*NFV*(NFV'') ≥ 0) = (1.10110, 1.12175)

I would stay away using Halley's Irrational formula.
To avoid issues, it required very close guess, more trouble than it is worth.

[Albert Chan]

(06-14-2022 02:59 AM)Albert Chan Wrote:  2 sign changes, thus 0 or 2 positive roots (for R = 1+i)
Based on rate edges, (PMT/FV , PMT/-PV) = (-1/50, 2/13)

2 Rates (if exist at all): -1/50 < i < 2/13

If we can show we have 1 root, we must have 2.

I was thinking of a fast way to locate extremum, check f sign, to decide number of roots.

Let decay factor, g = n*x / ((1+x)^n - 1)

f = (pv+fv)/n*g + pv*x + pmt

f' = (pv+fv)/n*g' + pv

f' = 0      → g' = -n*pv/(pv+fv)

g' is too complex to be invertible back for x
Instead, we approximate g with simple exponential decay function.
Matching g value and slope at x=0, we have g ≈ exp((1-n)/2*x)
exp is invertible with ln, we have:

Extremum: x_m ≈ m * ln(-m*n*pv / (pv+fv))      , where m = 2/(1-n)

Reusing previous example, n=30, pv=6500, pmt=-1000, fv=50000

f = 56500*x/((1+x)^30-1) + 6500*x - 1000

x_m ≈ -2/29 * ln(2/29 * 30 * 6500/56500) ≈ 0.09899

f(x = 0.09899) ≈ -6.46                        // true extremum, f(x = 0.1002) ≈ -6.52
f(x = 0) = (pv+fv)/n + pmt = 883.33

We have locked in 1 roots ⇒ we have 2 roots.