HP 41C solver problem

[Dieter]

(12-07-2017 05:46 AM)Thomas Okken Wrote:  OK, so the problem with the ROOT program from the HP-41C Standard Applications book is that its termination condition is not very good. It checks whether the result is zero, or whether its absolute value is less than 10^-10, and that doesn't always work well.

I assume we are talking about this program. Here the test is whether |f(x)|<1E–8. Which is less problematic than 1E–10 but still not a good idea at all.

Defining an exit condition for root finders can be tricky. Both the value of f(x) as well as the last change in x can be checked, i.e. the amount by which the current approximation differs from the previous one. In both cases a fixed value to compare against is a bad idea: For x²–2 the root may have been found if f(x)<1E-8, but what about 100000*(x²–2) or 0,000001*(x²–2)? The same is true for the current approximation: if x is around 5000 a change in the last digit is order of 1E–6 while for x~0,005 this is 1E–12.

So you better check the relative change in x, i.e. (xnew–xold)/xnew. If the absolute value of this drops below a given limit the program may exit and return xnew. For a 10-digit calculator this relative limit may be 1E–9, but since most root finding algorithms converge better than linear, I prefer an alternate method: if the relative change of x is less than 1E–6, do one last iteration an exit afterwards. The idea is: if the current relative error is 1E–6, the next iteration most probably will change x only in the, say, 9th or 12th digit, so the current approximation already has 8...11 digit accuracy.

The Standard Applications program uses a modified regula falsi method. It looks like the Illinois algorithm described here: calculate xnew from x1,y1 and x2,y2 and check if there is a sign change between x2 and xnew. If true, continue with the standard regula falsi method; if not, divide y2 by 2 first. This avoids certain problems mentioned in the linked Wikipedia article. But since the program only checks whether |f(x)|<1E–8 it may loop forever even if the root already has been found.

(12-07-2017 05:46 AM)Thomas Okken Wrote:  A better test would be to check what happens to the interval being searched.
(...alternate method...)
I'm guessing the ROOT program doesn't use this kind of algorithm because they wanted to keep it simple.

Here the Illinois method is an improvement over a standard secant method / regula falsi. The problem is that the program does not check the accuracy of x at all. It should terminate as soon as the last two approximations agree in, say, 8 significant digits.

(12-07-2017 05:46 AM)Thomas Okken Wrote:  I'd say the easiest way of making the ROOT program more useful is to insert a VIEW X just before the XEQ IND 03 on line 39. That way, you can see what it's doing, and it's easy to tell when the program gets stuck in an endless loop while actually having converged.

Definitely, yes.

For the time being the original program may be changed like this. It's still not perfect but it does converge. ;-)

Code:

.. ...
35 RCL 06
36 *
37 -
38 STO 04
39 LastX     // last absolute change in x    
40 RCL 04
41 X=0?
42 SIGN      // if x=0 check absolute change (sign(0)=1 on the HP-41)
43 /         // else check relative change
44 ABS
45 1E-8
46 X>Y?
47 GTO 04
48 VIEW 04   // show current approximation
49 RCL 04
50 XEQ IND 03
51 STO 07
52 RCL 06
53 *
54 X>0?
.. ...

Of course you can add an X=0? GTO 04 after line 51 to check if f(x) has become zero.
Edit: But this is not strictly required since in this case the next correction term automatically becomes zero and the program exits.

For the given function and initial guesses of 10 and 30 the program converges after seven iterations with X=23,81895808.
The Anderson-Björk method (cf. Wikipedia link) may converge slightly faster. It can be implemented easily by changing the code at LBL 01:

Code:

LBL 01
1
RCL 07
RCL 06
/
-
X<=0?
,5
ST* 05
GTO 02

Dieter