What's your favourite root finder? « Next Oldest | Next Newest »

 ▼ Les Wright Posting Freak Posts: 1,368 Threads: 212 Joined: Dec 2006 01-22-2007, 06:05 PM There has been a lot of discussion in the Forum over the months about quadrature. Recently, I have turned my attention to the admittedly less glamourous topic of root-finding. Even though I am aware of the sophisticated algorithms underlying the SOLVE functions in the HP41 Advantage Pac, HP34C, HP15C, HP42S, HP33S, etc, I must confess to really preferring elegant user code in a case like this, especially if I can manipulate it to show interim results and show the path of convergence. Given this, the SV function from the PPC ROM is about my favourite right now. Only 32 steps, a straightforward implementation of the secant method (which has pretty good convergence for a wise initial guess), and the ability to show or print out interim results. The authors of the SV routine were apologetic about its simplicity and limitations, but I have no idea why. For its small size, it really does the job, and is an exemplar of the tight and economical coding that the PPC ROM was famous for. I have played around with some others, like the Illinois modification of the regula falsi method in the HP Standard Applications Pac and the sort-of Newton-Raphson method in the HP11C owner's manual (I say "sort of" since the derivative is actually only numerically approximated rather than directly provided and computed), but SV is my favourite. For anyone interested, here is the routine for HP41C/CV/CX. I have replaced the single synthetic step with corresponding user code: ``` 01 LBL "SV" 02 STO 07 03 1 04 % 05 RCL Z 06 X=0? 07 X<>Y 08 STO 09 09 CLST 10 LBL 04 11 RCL Z 12 STO 08 13 RCL 07 14 FS? 10 15 VIEW X 16 XEQ IND 06 17 ST* 09 18 ST- 08 19 RCL 09 20 RCL 08 21 X#0? 22 / 23 STO 09 24 X<> 07 25 ST+ 07 26 RND 27 RCL 07 28 RND 29 X#Y? 30 GTO 04 31 RCL 07 32 END ``` You store the name of the function to be solved in register 6 (ASTO 06). Then, place on the stack the desired initial step size (or just 0 if you want to use the default of 1% of the initial guess) and the initial guess, then XEQ "SV". If you want to see interim or print results, set flag 10. Convergence is determined when the last two rounded results are equal so the setting of FIX, SCI, or ENG determines when the routine stops. What's your favourite root finder? Les ▼ Karl Schneider Posting Freak Posts: 1,792 Threads: 62 Joined: Jan 2005 01-23-2007, 04:06 AM Hi, Les -- Quote: What's your favourite root finder? Well, I still prefer built-in machine code instead of added-on keystroke programs, so I'm going with the HP-32SII implementation of SOLVE: It runs fast. The variables can be identified and chosen directly (unlike the original implementation on the "pre-Pioneer" models) It's easy to set up (unlike on the HP-42S). It can be used with both keystroke programs and equations. There's no direct-solution method that sometimes gives unexpected answers (unlike the HP-33S) -- KS David Smith Posting Freak Posts: 1,788 Threads: 36 Joined: Aug 2007 01-23-2007, 11:17 AM A pig named "Sweetiepie". ▼ Les Junior Member Posts: 4 Threads: 1 Joined: Jan 1970 01-24-2007, 12:22 PM I knew pigs were smart, but that is really Some Pig. Les Wright Posting Freak Posts: 1,368 Threads: 212 Joined: Dec 2006 01-24-2007, 06:35 PM Does Sweetiepie prefer Newton-Raphson, Brent, secant, Pegasus, or Illinois? Les ▼ David Smith Posting Freak Posts: 1,788 Threads: 36 Joined: Aug 2007 01-25-2007, 04:24 PM Sweetiepie prefers truffels, but potatoes, carrots, and radishes will do. Hold the onions and garlic. Valentin Albillo Posting Freak Posts: 1,755 Threads: 112 Joined: Jan 2005 01-25-2007, 09:46 AM Hi, Les: For the particular and extremely frequent case of finding the roots of polynomial equations, nothing beats PROOT, as found in the HP-71B's Math ROM. It will find and output all roots, real and/or complex, of any real polynomial of any degree whatsoever you may feed to it, without any input from the user at all (save the polynomial's coefficients, of course), at extremely fast speeds (Saturn's assembly language), and without ever failing. It can deal with roots of high multiplicity, and with such exotic beasts as zero leading coefficients, infinities, denormalized values, and NaNs. No polynomial has ever been produced which makes PROOT fail, thus no need to "show interim results and show the path of convergence" in this case. For example: ```Find all roots of x5-3*x4+8.1x2-1.37 = 0 Assuming you want to input the coefficients in array P and want the resulting roots returned in array R: MAT INPUT P -> P(1)=? 1,-3,0,8.1,0,-1.37 MAT R=PROOT(P) @ MAT DISP R will *very* quickly return the five roots: x1 = 0.423520562722 x2 = -0.428212993 x3 = -1.3019202156 x4 = 2.15330632294 + 1.07962690848i x5 = 2.15330632294 - 1.07962690848i ``` As for the general case, FNROOT, also in the HP-71B's Math ROM, can solve just about anything quickly (assembly language again), and allows itself to appear in your function's definition, thus being capable of solving non-linear systems of up to 5 unknowns. As for watching interim results or the path of convergence, just include the appropriate DISP or PRINT statements within your function's definition and there you are. You can even include some counter to let you know the number of iterations being performed. For example: ```Find a root of Cos(x) = x in [0,1], watching convergence: 1 DEF FNF(X) @ Y=COS(X)-X @ DISP X,Y @ FNF=Y @ END DEF 2 FNROOT(0,1,FNF(FVAR)) >FIX 5 @ RADIANS >RUN 0.00000 1.00000 1.00000 -0.45970 0.50000 0.37758 0.72548 0.02270 0.73840 0.00115 ... 0.74722 -0.01363 0.73909 0.00000 0.73909 (accurate to 12 digits, just 5 shown) Same without watching convergence: STD FNROOT(0,1,COS(FVAR)-FVAR) 0.739085133215 ``` If insisting in user code, them I'm afraid I'm partial to ye olde couple of mine productions, namely the polynomial solver which can find all roots (real and/or complex) of a given polynomial (which can also have complex coefficients, which PROOT doesn't allow), written in HP-41C user code and available in several user-produced ROMs (modified to take advantage of 41C Mcode complex words in the same ROM), and, for the general case, a very short but efficient and convenient routine I submitted for consideration to be included in the PPC ROM, also in HP-41C user code. It's small, fast, quite general, can be called as a stand-alone, prompting program or as a non-prompting subroutine, and even let's the user specify both required precision and a limit for the number of iterations in case of no roots, aborting the procedure and passing the fact to the calling program. Best regards from V. Edited: 25 Jan 2007, 9:54 a.m.

 Possibly Related Threads... Thread Author Replies Views Last Post [HP Prime] Using n-root symbol and exponent problem uklo 7 1,545 11-11-2013, 01:39 AM Last Post: Alberto Candel Cubic root (-8) = 2 ? Gilles Carpentier 37 5,104 08-12-2013, 10:26 PM Last Post: jep2276 Square Root Simplifier for HP39gII Mic 4 1,117 03-11-2013, 08:25 AM Last Post: Eddie W. Shore Cube root on standard calculator Thomas Klemm 22 3,396 11-09-2012, 06:11 AM Last Post: Pierre ROOT bug? HP 48S/48G Eddie W. Shore 8 2,589 07-13-2012, 07:05 PM Last Post: Eddie W. Shore x root y on hp42s David Griffith 14 2,134 04-08-2012, 12:43 PM Last Post: Walter B 35s prompt for multi-character variables in program like "low footprint" root finder Chris C 8 1,599 02-14-2012, 06:52 PM Last Post: Chris C [WP34S] SLV returns "Failed" when it finds a root just fine. Les Wright 6 1,192 11-27-2011, 04:31 PM Last Post: Paul Dale WP34S: Valentin Albillo's Polynomial Root Finder Miguel Toro 28 3,574 11-23-2011, 07:39 PM Last Post: Miguel Toro Results of new root-seeking methods Namir 3 772 08-22-2011, 03:28 PM Last Post: C.Ret

Forum Jump: