(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
 Thomas Klemm Senior Member Posts: 2,059 Joined: Dec 2013
(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)
halley-15c.zip (Size: 2.65 KB / Downloads: 7)
02-26-2022, 03:48 PM (This post was last modified: 02-26-2022 04:42 PM by Albert Chan.)
Post: #2
 Albert Chan Senior Member Posts: 2,551 Joined: Jul 2018
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.
06-11-2022, 04:12 PM
Post: #3
 Albert Chan Senior Member Posts: 2,551 Joined: Jul 2018
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
06-14-2022, 01:20 AM
Post: #4
 Gil Senior Member Posts: 586 Joined: Oct 2019
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?
06-14-2022, 02:59 AM
Post: #5
 Albert Chan Senior Member Posts: 2,551 Joined: Jul 2018
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.
06-14-2022, 07:33 AM
Post: #6
 Gil Senior Member Posts: 586 Joined: Oct 2019
RE: (15C) Halley's Method
Of course.
Thanks for reminding me.
Regards,
Gil
06-14-2022, 06:22 PM
Post: #7
 Albert Chan Senior Member Posts: 2,551 Joined: Jul 2018
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.
06-18-2022, 03:37 PM
Post: #8
 Albert Chan Senior Member Posts: 2,551 Joined: Jul 2018
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.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)