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.
