# HP Forums

Full Version: Halley-Ostriwski Root Method
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
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```
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :