Newton and Halley's methods with enhanced derivatives estimation

[Dieter]

(10-06-2017 05:39 PM)Gerson W. Barbosa Wrote:  I've found an example which causes your first program to loop forever on the HP-75C.

This may always happen since, due to limited numeric accuracy, the approximation may oscillate around a certain value with the last digit(s) varying by a few units.

(10-06-2017 05:39 PM)Gerson W. Barbosa Wrote:  But then my tolerance factor might not be adequate enough:

Setting a suitable exit condition indeed is crucial. Since the Newton method has quadratic and the Hally method even cubic convergence, a desired relative error of 1E–12 does not mean that the last approximation has to agree in 12 significant digits with the previous one. If the last two values differ by a relative error of 1E–8, this leads to an error estimate of 1E–16 for Newton and 1E–24 for Halley's method. In theory, that is. ;-)

That's why I prefer an implementation like it is done in some routines in the WP34s firmware that require an iterative approach (e.g. the Lambert W function): If an accuracy near 1E–30 is required, the iteration might exit at, say, 1E–20 (Newton) or even 1E–15 (Halley). If you like you may add one final iteration, just to be sure.

An application similar to yours (finding the interest rate in a TVM problem for the HP65, c.f. HP65/67 Software Library) exits at 1E–7 on a 10-digit calculator. Due to the limited accuracy (the HP65 does not offer ln(1+x) or e^x–1) further digits will be rounded off anyway as soon as 1+i is calculated. So knowing the properties and limitations of the considered equation can be essential.

In your example for a 12-digit machine I would try T=1E–8. Could you give it a try and see how well the returned solution agrees with the true result? Maybe after one additional final iteration?

(10-06-2017 05:39 PM)Gerson W. Barbosa Wrote:  The program loops around 0.06596460. The exact 16-digit result is 6.596460097781873E-2.

With Excel's 15-digit BCD precision it's the enhanced Halley method which, if the set tolerance is too small, oscillates around

0,0659646009778178
0,0659646009778200
0,0659646009778178
0,0659646009778200
...

As you can see even if x varies in its last three digits the value of f(x) is essentially the same, so more than 14 decimals or 13 significant digits do not make much sense. That's why for n-digit precision I would try a relative tolerance near 10^–(2/3*n). Add one final iteration if it makes you feel better. But that's just the way I would do it – you should do your own tests. ;-)

Dieter