 The Museum of HP Calculators

HP Forum Archive 20

 Results of new root-seeking methodsMessage #1 Posted by Namir on 22 Aug 2011, 12:08 p.m. 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

 Re: Results of new root-seeking methodsMessage #2 Posted by Marcus von Cube, Germany on 22 Aug 2011, 12:42 p.m.,in response to message #1 by Namir 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.

 Re: Results of new root-seeking methodsMessage #3 Posted by Namir on 22 Aug 2011, 2:30 p.m.,in response to message #2 by Marcus von Cube, Germany 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.

 Re: Results of new root-seeking methodsMessage #4 Posted by C.Ret on 22 Aug 2011, 3:28 p.m.,in response to message #1 by Namir 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 ‐ Rev 2.8: Mar. 22 2011 ‐ 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. Go back to the main exhibit hall