Post Reply 
(15C) Halley's Method
02-26-2022, 11:09 AM (This post was last modified: 06-18-2022 10:01 AM by Thomas Klemm.)
Post: #1
(15C) Halley's Method
Halley's Method

References

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


Attached File(s)
.zip  halley-15c.zip (Size: 2.65 KB / Downloads: 7)
Find all posts by this user
Quote this message in a reply
02-26-2022, 03:48 PM (This post was last modified: 02-26-2022 04:42 PM by Albert Chan.)
Post: #2
RE: (15C) Halley's Method
(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 Cf expressions are not the same.
First form denominator + is really ± (2 roots of quadratic)
Second form pick the smaller correction, exactly what we wanted.
Find all posts by this user
Quote this message in a reply
06-11-2022, 04:12 PM
Post: #3
RE: (15C) Halley's Method
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
Find all posts by this user
Quote this message in a reply
06-14-2022, 01:20 AM
Post: #4
RE: (15C) Halley's Method
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?
Find all posts by this user
Quote this message in a reply
06-14-2022, 02:59 AM
Post: #5
RE: (15C) Halley's Method
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.
Find all posts by this user
Quote this message in a reply
06-14-2022, 07:33 AM
Post: #6
RE: (15C) Halley's Method
Of course.
Thanks for reminding me.
Regards,
Gil
Find all posts by this user
Quote this message in a reply
06-14-2022, 06:22 PM
Post: #7
RE: (15C) Halley's Method
(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.
Find all posts by this user
Quote this message in a reply
06-18-2022, 03:37 PM
Post: #8
RE: (15C) Halley's Method
(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: xm ≈ 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

xm ≈ -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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: