Re: What's your favourite root finder? Message #7 Posted by Valentin Albillo on 25 Jan 2007, 9:46 a.m., in response to message #1 by Les Wright
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.
|