HP 41C solver problem

[Trond]

Hi folks!
I am currently testing out the equation solver (root finder) that came in the Standard Applications book. It seems to work, but then it got stuck on this:

0.0074X^3-100=0

My other two HP calculators find the answer relatively quickly (X=23.82), but the HP41CV just keeps going forever (well, several minutes until I stop it). Any idea why?

[Thomas Okken]

(12-07-2017 03:17 AM)Trond Wrote:  Hi folks!
I am currently testing out the equation solver (root finder) that came in the Standard Applications book. It seems to work, but then it got stuck on this:

0.0074X^3-100=0

My other two HP calculators find the answer relatively quickly (X=23.82), but the HP41CV just keeps going forever (well, several minutes until I stop it). Any idea why?

The function is smooth and monotonically increasing, and so provided your starting guesses bracket the one simple root, the ROOT program should always converge (or complain if the starting guesses are bad). I keyed in the program (taken from greendyk.nl) and your function, took a look at the function's plot with Wolfram Alpha, gave ROOT starting guesses of 0 and 40, and got the result: 23.8189580802. (Free42, not 41C, but that shouldn't make a difference).

My guess: you made a mistake entering the program. Double-check your code.

[Thomas Okken]

On second thought, trying this in Free42 rather than a 10-digit calculator could make a difference. The termination condition in ROOT might be a bit sharp. If you can't find any mistakes, try changing the 1e-10 in line 44 to something larger.

[Trond]

(12-07-2017 03:55 AM)Thomas Okken Wrote:  

(12-07-2017 03:17 AM)Trond Wrote:  Hi folks!
I am currently testing out the equation solver (root finder) that came in the Standard Applications book. It seems to work, but then it got stuck on this:

0.0074X^3-100=0

My other two HP calculators find the answer relatively quickly (X=23.82), but the HP41CV just keeps going forever (well, several minutes until I stop it). Any idea why?

The function is smooth and monotonically increasing, and so provided your starting guesses bracket the one simple root, the ROOT program should always converge (or complain if the starting guesses are bad). I keyed in the program (taken from greendyk.nl) and your function, took a look at the function's plot with Wolfram Alpha, gave ROOT starting guesses of 0 and 40, and got the result: 23.8189580802. (Free42, not 41C, but that shouldn't make a difference).

My guess: you made a mistake entering the program. Double-check your code.

I checked the code. Seems to be OK

I also checked a few similar equations. It found the solution quickly.

The moment I added the 0.0074 factor it slows down to a crawl, or it can't find it at all. Very bizarre.

[Trond]

Yup, it was able to find 2X^3-100=0

But not 0.0074X^3-100=0

The latter either takes so long that I lose my patience, or I wonder if a factor such as 0.0074 somehow can throw the possible solutions off enough (by rounding or something?) that it keeps searching forever. I am no expert in this however.

[AlexFekken]

According to my Prime:

f(23.8189580802) = -8*1E-10
f(23.8189580803) = 10*1E-10

Better do what Thomas has already suggested...

[Trond]

(12-07-2017 04:23 AM)Thomas Okken Wrote:  On second thought, trying this in Free42 rather than a 10-digit calculator could make a difference. The termination condition in ROOT might be a bit sharp. If you can't find any mistakes, try changing the 1e-10 in line 44 to something larger.

Worked! I changed it to 1 E-6

[Trond]

Thanks guys!
Glad that it wasn't only myself messing everything up.

By the way, the equation is based on a little approximation I made on whale weight:

W=CL^3, (Weight in metric tons, Length in meters) where C seems to be about 0.0074 for rorqual whales, but different for whales of different shape (almost double in the much chubbier bowhead or right whale. If you nitpick the exponent should vary slightly, but it is never far off from 3).

Playing around with all of those variables in this equation is easy in HP 42s (if I change one variable input, all of the others will change accordingly), but I am not sure if there's a program, or program combo, that can do it similarly in HP 41CV without constantly reprogramming each time I change a parameter?

[Thomas Okken]

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.

In the case of 0.0074 * x^3 - 100, it works badly because the first and second terms of that formula are both near 100 at the root, so the actual root loses two or three digits of precision, so the "less than 10^-10" test is useless: unless the function evaluates to exactly zero, it is guaranteed to evaluate to something greater than 10^-10.

A better test would be to check what happens to the interval being searched.

Consider:

Given guesses x1 and x2, where f(x1) < 0 and f(x2) > 0, the secant method provides a new guess x3, which is the point at which the straight line through (x1, f(x1)) and (x2, f(x2)) intersects the x axis.

If f(x3) = 0, you're done.
If f(x3) < 0, set x1 = x3, and repeat.
If f(x3) > 0, set x2 = x3, and repeat.

But what if x3 is less than or equal to x1, or greater than or equal to x2? Mathematically, that's impossible, but in finite-precision math, it is possible.

In that case, you can fall back on the bisection method: just try x3 = (x1 + x2) / 2 instead, and then do the same check on f(x3) = 0, f(x3) < 0, and f(x3) > 0, as above. In this case, x3 can never be less than x1 or greater than x2, and the only way it can be equal to x1 or x2 is if x1 and x2 are as close as they can be, to each other, given the floating-point representation in use, and so, x3 = x1 or x3 = x2 provides a robust termination condition.

I'm guessing the ROOT program doesn't use this kind of algorithm because they wanted to keep it simple.

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.

[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

[Thomas Okken]

(12-07-2017 10:39 AM)Dieter Wrote:  

(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.

Yes. My source is http://www.greendyk.nl/hp41c-applications/index.php (page 41), which looks identical except for that one constant. Maybe an earlier edition of the book?

[Trond]

(12-07-2017 10:39 AM)Dieter Wrote:  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?
.. ...

I like it! Thanks for the code. Maybe I'll try Anderson-Björk later.

[Dieter]

(12-07-2017 01:56 PM)Thomas Okken Wrote:  My source is http://www.greendyk.nl/hp41c-applications/index.php (page 41), which looks identical except for that one constant. Maybe an earlier edition of the book?

Before I walk down to the, err... "archive room" and try to find the Standard Applications book that came with my early (tall key) 41C: maybe someone has it at hand and can look up this program.

But even the website you linked to says this:

Quote:The answer is labeled and displayed when the value of the function is less than 10^-10 closer tolerance can be obtained simply by keying in a different value when the program is entered.

But again: checking only the function result is a bad idea.

Dieter

[Guenter Schink]

(12-07-2017 06:10 PM)Dieter Wrote:  

(12-07-2017 01:56 PM)Thomas Okken Wrote:  My source is http://www.greendyk.nl/hp41c-applications/index.php (page 41), which looks identical except for that one constant. Maybe an earlier edition of the book?

Before I walk down to the, err... "archive room" and try to find the Standard Applications book that came with my early (tall key) 41C: maybe someone has it at hand and can look up this program.

But even the website you linked to says this:

Quote:The answer is labeled and displayed when the value of the function is less than 10^-10 closer tolerance can be obtained simply by keying in a different value when the program is entered.

But again: checking only the function result is a bad idea.

Dieter

My Standard application book (German) printed 1982 also has the value 1E-10. And yes it states that this value is subject to adjustment.

Günter

[Massimo Gnerucci]

(12-07-2017 06:10 PM)Dieter Wrote:  

(12-07-2017 01:56 PM)Thomas Okken Wrote:  My source is http://www.greendyk.nl/hp41c-applications/index.php (page 41), which looks identical except for that one constant. Maybe an earlier edition of the book?

Before I walk down to the, err... "archive room" and try to find the Standard Applications book that came with my early (tall key) 41C: maybe someone has it at hand and can look up this program.

But even the website you linked to says this:

Quote:The answer is labeled and displayed when the value of the function is less than 10^-10 closer tolerance can be obtained simply by keying in a different value when the program is entered.

But again: checking only the function result is a bad idea.

Dieter

My Standard application book (italian) printed 7.79 also has the value 1E-10 @line #44. And yes it states that this value is subject to adjustment.

[Massimo Gnerucci]

Having brought out of the drawer my almost 40 y.o. HP-41C box I was amazed, once more, by the quality of the materials (the box, manuals, module holder, calculator, case, even the Eveready N cells box!, are still in like pristine condition).
And I was reminded of fond memories, while devouring that complete and easily readable manual, coast to coast.

Thank you Hewlett·Packard from the Seventies.

[rprosperi]

(12-07-2017 08:15 PM)Massimo Gnerucci Wrote:  Having brought out of the drawer my almost 40 y.o. HP-41C box I was amazed, once more, by the quality of the materials (the box, manuals, module holder, calculator, case, even the Eveready N cells box!, are still in like pristine condition).
And I was reminded of fond memories, while devouring that complete and easily readable manual, coast to coast.

Thank you Hewlett·Packard from the Seventies.

Thank God Massimo... I thought I was the only geek that kept the original battery box!!

[Thomas Okken]

Checking TOS, the Spanish, Italian, and German versions of the Standard Applications book there all use the 1e-10 limit, while the English version, and the STRD 1C ROM image both have 1e-8.

The Spanish and Italian manuals are dated July 1979, while the English and German ones are dated August 1980.

The English manual includes this addendum:

Quote:A final paragraph for the description of "Root Finder" should read as follows:

This program calculates close approximations to a root, but may continue to iterate when the magnitude of the function evaluated at these approximations exceeds the tolerance.
You can check the progress of the solution by inspecting the current guesses in registers 1 and 2 using the VIEW function. You may find it convenient to assign VIEW to some key.

To monitor the progress of the root finder at any time press [R/S]
then [VIEW] 01 to see X1
and [VIEW] 02 to see X2
resume execution by [R/S]

[Dieter]

(12-08-2017 12:03 AM)Thomas Okken Wrote:  Checking TOS, the Spanish, Italian, and German versions of the Standard Applications book there all use the 1e-10 limit, while the English version, and the STRD 1C ROM image both have 1e-8.

Sounds like a bugfix.

(12-08-2017 12:03 AM)Thomas Okken Wrote:  

Quote:You can check the progress of the solution by inspecting the current guesses in registers 1 and 2 using the VIEW function. You may find it convenient to assign VIEW to some key.

Assign ?!? – VIEW is on the HP41 standard keyboard (Shift [R/S]).

Dieter

[Thomas Okken]

(12-08-2017 07:02 AM)Dieter Wrote:  

Quote:You can check the progress of the solution by inspecting the current guesses in registers 1 and 2 using the VIEW function. You may find it convenient to assign VIEW to some key.

Assign ?!? – VIEW is on the HP41 standard keyboard (Shift [R/S]).

The point seems to be to steer people toward using VIEW to check those two registers, rather than RCL, which would require one less keystroke but also alter the stack and mess up the iteration.

Fixing the program properly would have been better!