HP Forums

Full Version: Results of new root-seeking methods
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

Hi All,

Fellow forum members Hugh Steer and Lyuka have opened an interesting door. Hugh mentioned the Ostrowski algorithm for finding nonlinear function roots. I went on a little online crusade to find more on the Ostrowski algorithm. I found several improvements (some just a few years old) on Newton's method and Ostrowski's methods, and a few improvements on Halley's methods. I came out with about 30 algorithms/variants starting with Newton, Halley, Ostrowski, and so on. I included the version of Ostrowski algorithm that Luyka mentioned in his web site.

I tested all these methods with many nonlinear functions. The results???

1. The simple Ostrowski algorithm ranked the best in the battery fo tests. Its simplicity (compared to the sophisticated enhancements given by various mathematicians) kept the number of function calls and iterations low.

2. The Ostrowski 4th Order (by Grau and Diaz-Barrero) was ranked second best.

3. The Ostrowski algorithm as implemented on Lyuka's web site came third. A special mention to the way the Ostrowski algorithm is implemented since it had only one function call per iteration!!! Very clever!!

Here is the basic scheme used by the simple Ostrowski algorithm:

y = x - f(x) / f'(x)
x = y - f(y)(x - y) / (f(x) - 2 f(y))

Here is the basic scheme used by he Ostrowski 4th Order (by Grau and Diaz-Barrero):

y = x - f(x) / f'(x)
u = (x - y) / (f(x) - 2 f(y))
z = y - f(y)
x = z - u * f(z)

Check Lyuka's web site algorithm version and implementation here

Edited: 22 Aug 2011, 11:11 p.m. after one or more responses were posted

I would rank the methods based on different criteria:

1. How many function evaluations are necessary for a given accuracy? The reasoning behind this is that the evaluation of f is more expensive than the operations that the algorithm does with the values. If algorithm (a) uses twice as many function evaluations as algorithm (b) it's ranked better only if the number of iterations is less then half of that of algorithm (b).

2. How stable is the algorithm for some nasty functions or with bad initial guesses? Is it able to solve more problems then some other algorithm? E. g. Newton will fail if started at or near an extremum.

Marcus,

I have ranked the various algorithms based on the number of function calls. Some of the sophisticated enhancements on Newton or Ostrowski require 4 or 5 function calls!! Several of these methods will give a refined root value (within the specified accuracy) in two or three iterations.

As far as stability, again some of the sophisticated enhancements become unstable depending on the combination of test function and the initial values.

Thank you Namir for pointing me out this stable and fast root-finding method. I ‘m now using it instead of a more classic but slow search by dichotomy algorithm.

The code for HP-42S (ref. OST &dash; Rev 2.8: Mar. 22 2011 &dash; Copyright(c) 1998-2011 Takayuki HOSODA, All rights reserved.) is great as it returns root, value,status and give the oportuniti to select with register to solve for.
But, it is a bit too much sophisticate for my own application; Such I compose a personal version which is much less sophisticated. In counterpart, usage is simplest and number of lines is reduced to 70 only. But there no more ‘status flag’, and the user have to always start with two guess. Then the code is finding a root, run forever or die.

Before a search can start, the two registers R04 and R05 have to be initialized respectively with expected accuracy eps (i.e. 1E-6 is a good enough) and the name of the function (i.e. “FCT1” as following exemple).

The function have to be programmed in the calculator as taking its only argument from the stack x register and return the value of the function at this point into the same x stack-register. The function code may use other stack registers (convenient for complicated function lucubrations) as the code always consider that at each call of the function, all content of the stack as well as LastX are lost, except the x register returning f(x).

-- - PROGRAM  -		  ----- STACK REGISTER -----
-- ------------ t: z: y: x:
01 LBL"FCT1 ~ x
02 E^X ~ exp(x)
03 LastX ~ exp(x) x
04 x^2 exp(x) x^2
05 3 ~ exp(x) x^2 3
06 * ~ ~ exp(x) 3.x^2
07 - ~ ~ ~ exp(x)-3.x^2
08 RTN ~ ~ ~ f(x)
-- ---------------

Registers to be initialized :
R04: eps
R05: "F"

Used Registers :
R00: i evaluations/iteration counter
R01: xn R06: yn = f(xn)
R02: xn+1 R07: yn+1 = f(xn+1)
R03: xn+2 R08: yn+2 = f(xn+2)
R04: eps expected accuracy
R05: "FCT1" name of function program/code
R09: -t last value of t used for computing next x_n+3 position.

-- - PROGRAM  -
-- ------------ -- -------------- -- --------------
01 LBL"OSTRO 26 RCL 01 51 STO 09
02 STO 02 27 RCL 03 52 *
03 STO 03 28 RCL 02 53 RCL 01
04 x<->y 29 STO 01 54 +
05 STO 01 30 - 55 1
06 ST+ 03 31 RCL 03 56 RCL 09
07 XEQ IND 05 32 STO 02 57 +
08 STO 06 33 RCL Z 58 /
09 RCL 02 34 - 59 STO 03
10 XEQ IND 05 35 / 60 ARCL 03
11 STO 07 36 RCL 06 61 AVIEW
12 2 37 / 62 XEQ IND 05
13 STO 00 38 RCL 08 63 STO 08
14 ST/ 03 39 RCL 06 64 ABS
15 RCL 03 40 - 65 RCL 04
16 XEQ IND 05 41 * 66 x<=y ?
17 STO 08 42 RCL 07 67 GTO 00
18 LBL 00 43 STO 06 68 RDN
19 ISG 00 44 * 69 RCL 03
20 NOP 45 RCL 08 70 END
21 CLA 46 STO 07
22 FIX 0 47 RCL 06
23 ARCL 00 48 -
24 ~": 49 /
25 FIX 6 50 CHS
-- -------------- -- -------------- -- --------------

Structure of program:
Instruction 01 to 17 :
- initialize search; input and store x0 and x1,
- compute y0=f(x0), y1=f(x1), as well as first intermediate point x2=(x0+x1)/2 and y2=f(x2).
- Initialized iteration counter (R00).

Instructions 18 to 67 :
- Main loop composed by two tasks , both combine in one code:
compute next point xn+3 from [italic]t
and translate registers so that value of xn-…and [italic]yn-... are at the correct registers for next iteration.
- Program loops until |f(xn)| is small enough, i.e. less than expected accuracy eps
- During root-search, the calculator display iteration number i and best root estimate xn)| so that user can observe progress.

Instructions 68 to 70:
- recall best root estimate in stack register x and end program.

Edited: 22 Aug 2011, 3:38 p.m.