+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (15C) Halley's Method (/thread-18082.html)
(15C) Halley's Method - Thomas Klemm - 02-26-2022 11:09 AM
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 IB: 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 RTNC: 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 RTND: 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 RTNE: 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 RTNRE: (15C) Halley's Method - Albert Chan - 02-26-2022 03:48 PM
(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.
RE: (15C) Halley's Method - Albert Chan - 06-11-2022 04:12 PM
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
RE: (15C) Halley's Method - Gil - 06-14-2022 01:20 AM
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?
RE: (15C) Halley's Method - Albert Chan - 06-14-2022 02:59 AM
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.
RE: (15C) Halley's Method - Gil - 06-14-2022 07:33 AM
Of course.
Thanks for reminding me.
Regards,
Gil
RE: (15C) Halley's Method - Albert Chan - 06-14-2022 06:22 PM
(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.
RE: (15C) Halley's Method - Albert Chan - 06-18-2022 03:37 PM
(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.