HP 41C solver problem

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