Halley-Ostriwski Root Method

[Namir]

The algorithm for the Halley-Ostrowski method that I developed is:

Code:

F0 = f(x)
Do 
  h = 0.01 * (1 + |x|)   
  Fp = f(x + h) 
  Fm = f(x - h) 
  Deriv1 = (Fp - Fm) / 2 / h 
  Deriv2 = (Fp - 2 * F0 + Fm) / h / h 
  Diff = F0 / Deriv1 / (1 - F0 * Deriv2 / Deriv1 / 2 / Deriv1) 
  z = x - Diff 
  Fz = f(z) 
  If |x – z| < h Then h = x - z 
  Deriv1b = (F0 - 2 * Fz) / (x - z) 
  Deriv2b = (Fp - 2 * Fz + Fm) / h / h 
  Diff2 = Fz / Deriv1b / (1 - Fz * Deriv2b /  
              Deriv1b / 2 / Deriv1b) 
  x = z – Diff2 
  F0 = f(x)
Loop Until |F0| < Toler  
Return X as the refined guess for the root.

The HP-41C implementation uses the following memory map:

Code:

R00 = X
R01 = FxToler
R02 = h
R03 = F0
R04 = Fp
R05 = Fm
R06 = Deriv1, Deriv1b
R07 = Deriv2, Deriv2b
R08 = z
R09 = Fz
R10 = Diff, Diff2

The HP-41 listing is:

Code:

001    LBL "HALOST"        
002    LBL A        
003    "X?"        
004    PROMPT        
005    STO 00        
006    CF 22        
007    "FXTOLER?"        
008    PROMPT        
009    FC? 22        
010    1.00E-08        
011    STO 01        
012    RCL 00        
013    XEQ E        
014    STO 03        # F0 = f(X)
015    LBL 00         # start of main loop
016    RCL 00        
017    PSE        
018    ABS        
019    1        
020    +        
021    0.001        
022    *        
023    STO 02        # h = 0.01 * (1 + |x|)
024    RCL 00        
025    RCL 02        
026    -        
027    XEQ E        
028    STO 05        # Fm = f(X-h)
029    RCL 00        
030    RCL 02        
031    +        
032    XEQ E        
033    STO 04        # Fp = f(X+h)
034    RCL 05        
035    -        
036    2        
037    /        
038    RCL 02        
039    /        
040    STO 06        # Deriv1 = (Fp - Fm) / 2 / h
041    RCL 04        
042    RCL 05        
043    +        
044    RCL 03        
045    2        
046    *        
047    -        
048    RCL 02        
049    X^2        
050    /        
051    STO 07        # Deriv2 = (Fp - 2 * F0 + Fm) / h / h 
052    RCL 03        
053    RCL 06        
054    /        
055    1        
056    RCL 03        
057    RCL 07        
058    *        
059    RCL 06        
060    X^2        
061    /        
062    2        
063    /        
064    -        
065    /        
066    STO 10        # Diff = F0 / Deriv1 / (1 - F0 * Deriv2 / Deriv1 / 2 / Deriv1) 
067    CHS        
068    RCL 00        
069    +        
070    STO 08        # z = x - Diff
071    XEQ E        
072    STO 09        # Fz = f(z)
073    RCL 00        
074    RCL 08        
075    -        
076    ABS        
077    RCL 02        
078    X<>y        
079    X<Y?        
080    GTO 01        
081    RCL 00        
082    RCL 08        
083    -        
084    STO 02        
085    LBL 01        
086    RCL 03        
087    RCL 09        
088    2        
089    *        
090    -        
091    RCL 00        
092    RCL 08        
093    -        
094    /        
095    STO 06        # Deriv1b = (F0 - 2 * Fz) / (x - z)
096    RCL 09        
097    RCL 06        
098    /        
099    1        
100    RCL 09        
101    RCL 07        
102    *        
103    RCL 06        
104    X^2        
105    /        
106    2        
107    /        
108    -        
109    /        
110    STO 10        # Diff2 = Fz / Deriv1b / (1 - Fz * Deriv2b / Deriv1b / 2 / Deriv1b)
111    CHS        
112    RCL 08        
113    +        
114    STO 00        
115    XEQ E        
116    STO 03        # F0 = f(X)
117    ABS        
118    RCL 01        # ABS(f(x)) > Toler
119    X<Y?        
120    GTO 00        # end man loop
121    "RT="        
122    ARCL 00        
123    PROMPT        
124    GTO A        
125    LBL E        # calculate f(x)
126    EXP        
127    LASTX        
128    X^2        
129    3        
130    *        
131    -        
132    RTN