On Convergence Rates of Root-Seeking Methods

[Gerson W. Barbosa]

(03-04-2018 11:58 PM)Mike (Stgt) Wrote:  A german proverb says 'no reply is also a reply'.
/M.

So is a reply from a non-expert.

Probably not the best methods, but lots of information I was not aware of:

https://en.wikipedia.org/wiki/Methods_of...uare_roots

a) Bakhshali method: interesting as the source is an ancient Hindu manuscript;

b) Taylor Series: quoting from Wikipedia:

Quote:As an iterative method, the order of convergence is equal to the number of terms used. With two terms, it is identical to the Babylonian method. With three terms, each iteration takes almost as many operations as the Bakhshali approximation, but converges more slowly.

Despite that remark, that's the one I chose to play with, because it is easy to follow. Also, it can be expanded at will if one wants to try it with more terms.

The first three terms:

\(\sqrt{x}= r\left ( 1+\frac{x-r^{2}}{2r^{2}}-\frac{\left ( x-r^{2} \right )^{2}}{8r^{4}} \pm \cdots \right )\)

Equivalently, after a manual simplification:

\(\sqrt{x}\approx \frac{3\left ( r^{2}+2x \right )-\left ( \frac{x}{r} \right )^{2}}{8r}\)

Decimal Basic program and output:

OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET t = TIME
INPUT x
LET r = 10^INT(LOG10(x)/2) + 1 ! initial estimate
PRINT r
DO 
   LET a = r
   LET r = (3*(r^2 + 2*x) - (x/r)^2)/(8*r)
   PRINT r
LOOP UNTIL ABS(r - a) < 1E-300
PRINT SQR(x) - r
PRINT TIME - t;" seconds"
END


? 17
 2 
 2.609375
 3.831457666622230815360207907180368367929063526645300241783992701606003096535603​1813886658895626223976443827895... 
 4.122244620637600419719958337875107818477074015229252100479776053793938950195150​0012131448438220412295740265140...
 4.123105625598878536178011647766885826496727434593597792311874268691149149019651​3723285709444226805259084101669...
 4.123105625617660549821409855974076830276053739312118674761681625480610941537298​4751649808177898214772303507785...
 4.123105625617660549821409855974077025147199225373620434398633573094954346337621​5935878636508106842966843263884...
 4.123105625617660549821409855974077025147199225373620434398633573094954346337621​5935878636508106842966845440409...
 4.123105625617660549821409855974077025147199225373620434398633573094954346337621​5935878636508106842966845440409...
 8.202714670745631983631354915192547370538210924268951508828665917189926704332598​08358492948699826676293788E-931...
2.02 seconds


This is a test for small argument as the initial estimation method will cause error for large arguments ( starting around 920 ).

RPL program

%%HP: T(3)A(R)F(.);
\<< DUP XPON 2. / ALOG 1. +
DO DUP2 DUP2 SQ SWAP DUP + + 3. * UNROT / SQ - OVER 8. * / SWAP DUP2 / 1. - ABS NIP .00000001 <
UNTIL
END NIP
\>>

Again, this will cause error for x = 0. No problem, we know Sqrt(0) = 0. Also, these are just tests :-)