# 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.