Has anyone here had occasion to ever work with the adaptation of Newton's Method found at p. 154 of the 11C Owners Manual?

I have been intrigued by the routine because it uses a numerical approximation to the derivative so the user doesn't have to provide it. It goes a little overboard in terms of convergence criteria--the routine stops when both the change in the root between iterations and evaluation of the function at the approximated root fall within specified tolerances.

I have ported the routine to the HP41 and have found certain circumstances where the numerical approximation to the derivative eventually leads to division by zero error just before convergence since the function value at x and x + delta are essentially the same. This can be rectified by setting the desired tolerance bigger by at least an order of magnitude.

That said, the routine is pretty pokey compared the secant method and modified regula falsi method I have mentioned in another thread, since there are two function evaluations for each iteration rather than just one. It is also unnecessarily long and convoluted--there is no need for two convergence criteria or loop counter. I have reduced the size somewhat by cutting back on this. I have inserted a step that shows or prints interim results as desired.

If anyone has had occasion to work with this particular routine from our literature, I would like to hear about your experience with it.

Les

Les,

I use a practical approximation to the derivative of f(x):

f'(x) = (f(x+h) - f(x)) / h

where h = 0.01 * (1 + |x|)

Newton's algorithm which is:

x = x - f(x) / f'(x)

becomes:

x(new) = x - h * f(x) / (f(x+h) - f(x))

I usually calculate f(x) first and store it since it appeas twice in the above equation.

The iterations continue until the absolute value of the refinement in x (that is the absolute value of h * f(x) / (f(x+h) - f(x))) is less than a tolerance value OR the iteration has exceeded the maximum limit.

Namir

What's nice about Newton's Method is that it will converge extremely rapidly -- but only if the function in question is well-behaved around the root, and you also need to give it a reasonably close starting value, or else it may bounce around or wander off to infinity.

When Isaac Newton devised his Method, he had a very specific purpose in mind, namely, to solve for certain parameters of planetary motion. Just Google for newton's method planetary motion and you'll get the idea. The point being, he understood the functions in question perfectly, and knew their behavior around the roots he was looking for, so he could be confident that his iteration scheme would converge -- and quickly, too, which was a pretty big deal for a person like him, who had to perform all his math by hand!

For a general-purpose numerical root finder, on the other hand, Newton's Method is simply inappropriate. If you don't have the time to do the kind of analysis that would prove that NM will converge, you need to use something more robust. Unfortunately, the literature about this kind of robust root-finding is a bit scarce -- I suppose it's because the people who need to integrate numerical root finding into some software product usually have the ability to analyze their problems to the point where they can identify one of the classic approaches to be adequate; which leaves only the "general purpose" kind of software, like Matlab or Mathematica or (of course!) HP calculators, whose developers tend to keep their numerical know-how to themselves...

You may find the Free42 implementation of SOLVE to be of some interest. It uses a mixed approach, using the secant method to find two values of X such that f(x1) and f(x2) have opposite signs (in case the user hasn't supplied two such values already), and then uses either exponential fitting, the secant method, or simple bisection, to narrow down the interval as far as it can, and switches between those modes whenever necessary. The Free42 code for this is not very easy to read, largely because of the way I had to structure it so that it could run in non-preemptive single-threaded environments, but I'll be happy to explain it in detail if you're interested.

- Thomas

Thanks, Namir, and the email that preceded it.

I turned out pretty quickly a rough little routine that works fine both for the 41 and, with a tiny bit of modification, the 11C. It is about 40 steps, and with a little more ingenuity I could make it shorter. I find this interesting since the 11C program in the manual that does the exact same thing is 68 steps. And your approximation to the derivative behaves a lot better and division by zero errors are largely avoided.

I have put in a flag 10 dependent way to print out interim results, analagous to IG and SV in the PPC ROM, since I like that sort of thing.

Thanks for sharing with me.

Les

Hello Les,

What do you think of a 35-step RPN program written by Newton himself? ;-)

If you don't like fiction like I do, proceed directly to page 9 for the HP-10C program listing. It runs on the HP-33C with no modification.

There's no need to credit the real author here as you can tell by the link:

http://membres.lycos.fr/albillo/calc/pdf/DatafileVA004.pdf

Regards,

Gerson.

Hi, Gerson:

Thanks for directing Les (and forum's visitors in general) to my Newton's Anniversary commemorative short tale, I hope anyone bold enough to actually read it will be suitably pleased, as Datafile readers were at the time.

Perhaps Les will also find interesting this Root Finder routine for the HP-41C, which I submitted more than 25 years ago to PPC to be considered for inclusion in the then future PPC ROM. Per my stated objectives, it was intended to be:

• interactive, with full prompting (8 steps), and also callable as a non-prompting subroutine (35 steps, including the final END)

• very short (8 + 35 steps) and fast, with user-selectable precision by setting the desired display mode.

• suitably "ruggerized", so that a zero derivative wouldn't cause and error or premature termination but would simply redirect the search elsewhere

• controlled to avoid infinite looping in case of no roots, returning back to the calling program with a condition set so that the calling program can test for it and proceed as needed (perhaps try another guess or try)

• fully optimized, including a worthwhile trick or two (such as using two D-R in succession to multiply the current guess by a very small increment (thus actually being a relative increment) in order to compute a fairly good approximation to the derivative. This is just 2 bytes, instead of "{some constant} *", which would be several bytes for the constant, plus one for the multiplication (or else you'd need to previously store the constant in some register, which would be even more expensive)
  

I'm including here a reasonably legible scan with the listing, the description, the instructions, and the examples.

Edited: 29 Jan 2007, 8:51 a.m.

Thank you so much!

It took me just a few minutes to transcribe the routine into a textfile and convert it to both RAW ans barcode format. Here is the listing in text form for conveniece:

 01 LBL "RF"

 02 "NAME?"

 03 AON

 04 PROMPT

 05 AOFF

 06 ASTO 00

 07 "X0?"

 08 PROMPT

 09 LBL "RFP"

 10 CF 00

 11 STO N

 12 50

 13 STO O

 14 LBL 00

 15 RCL N

 16 1

 17 D-R

 18 D-R

 19 +

 20 XEQ IND 00

 21 STO M

 22 RCL N

 23 XEQ IND 00

 24 X=0?

 25 GTO 01

 26 ST- M

 27 RCL M

 28 X=0?

 29 SIGN

 30 /

 31 D-R

 32 D-R

 33 ST- N

 34 RND

 35 X=0?

 36 GTO 01

 37 DSE O

 38 GTO 00

 39 SF 00

 40 LBL 01

 41 RCL N

 42 CLA

 43 END

I really like being able to encourage faster convergence by using a less precise setting of FIX, which will affect the decisions made by the steps "RND X=0?" that crop up a couple of times. This is sometimes necessary with some functions or some initial guess to get the routine to stop, since due to rounding the last (usually 10th) SD of the root or function value at the root wiggles around a bit between iterations. Choosing FIX 7 or FIX 8 or even something less will fix this, and typically seems to result in convergence to the root that gets you as close to a zero function value as you can given the accumulated rounding error.

I wrote my own little routine with what Namir provided, but it is still 40 steps or so without the prompting stuff at the beginning. I have made your routine my favourite Newton routine and have enshrined it on its very own magnetic card. In the tradition of the PPC ROM I actually prefer the nonprompting subroutine. You should be proud that RFP is only three steps longer than the SV roujtine that is in the PPC ROM. Like SV for my own use I have added two steps to RFP (not in the listing above) that will show or print interim results if flag 10 is set.

Many thanks for directing us to this!

Les

I have also noticed that RFP is short enough to be adapted to the 33C while leaving enough space for functions to be solved. For that adaptation I would shorten it further by removing the 50 iteration criterion--actually hard to implement on the 33C without the DSE function. With all of the steps saved there you could insert a PSE if desire to watch for convergence oneself.

Edited: 30 Jan 2007, 4:33 a.m.

Quote:

---

I have also noticed that RFP is short enough to be adapted to the 33C while leaving enough space for functions to be solved

---

Actually, maybe not. The 33C doesn't have an RND function, which is crucial to the conciseness of the routine.

Hi, Les:

Thanks a lot for your kind words and appreciation of my humble (and old!) root-finding routine for the HP-41C. Despite its simplicity, it certainly delivers and I've used it myself as a subroutine called from a number of programs with perfectly
satisfactory results.

As I see you're keen on numerical analysis topics such as definite integration and root finding, I fully recommend that
you don't miss my future articles on advanced techniques for
both, belonging to a new series due next in Datafile. Very
powerful stuff, as you shall see.

P.S.: BTW, this next January-February issue of Datafile does
feature no less than three articles by me, introducing two
new regular series, among many other worthwhile contents.
I'll post a brief abstract as soon as the issue is released.

Quote:

---

It runs on the HP-33C with no modification.

---

In the absence of the trusty RND function on the 33C, does the conditional test x=0? compare the ACTUAL value of the X-register to zero, or does it compare the value rounded according to the setting of FIX?

I ask because with that excellent routine there will be times when the difference computed at step 21 will not be exactly zero, and the routine won't converge according to that criterion.

Maybe I should just try it out myself!

Les

Hi again, Les:

The routine featured in my short story is perfectly good and optimal for the HP-10C only, where it must cope with the lack of RAM, the lack of subroutine capability (as it must call your f(x) twice per iteration), and the lack of sorely needed conditionals, which must be mimicked by using "X=0?" and such.

For any other HP calc model not hindered by such appalling restrictions, this routine is sub-optimal and you would do much better with some adaptation of my HP-41C routine instead,
for instance.

I tried it out, and found that, alas, on the 33C, zero means zero internally, not the displayed value.

This equation turns up in Acton's Numerical Methods that Usually Work. On the 33C:

ENTER, x^2, 1, +, 2, /, 29.672, /, -

This fits just barely after the 34 step routine on the 33C, though I guess one could save steps by prestoring the constant 29.672.

This is basically a quadratic with two real roots, approx 0.017 and 59.3. If you start with a guess of 31, the routine won't converge but will wiggle back and forth between 59.32714430 and 59.32714429. Start with a better guess of 59.3, and the routine stops at the correct 10-digit root of 59.32714431. If the 33C had a RND function, the routine could be made to stop by setting FIX 7.

The rounding error can be fussy. If I start with the ostensibly even better guess of 59.327, I get the same wiggling as I got with a starting guess of 31. But if I choose an even worse guess of 30, it converges and stops at 59.32714431!

Weird. (Actually, not weird, as I know it has to do with how the rounding error accumulates according to how the root is approached.)

Thanks for sharing these amazing little routines. I note that they are half the length of the routine in the 11C manual, and I think are more efficient.

Les