HP Forums

Hi All,

I am searching for good pseudo-code or source code or very clear step-by-step explanation for the Lehmer-Schur algorithm used in calculating polynomial roots by applying a special version of root-bracketing method.

Any links or info (beyond Wikipedia and the general algorithm description in Acron's book "Numerical Methods that Work") is greatly appreciated.

Thanks!

Namir

Edited: 22 Aug 2010, 3:02 p.m.

Thank you Google! :-)

Follow the links...

Greetings,
Massimo

Thanks .. I have already seen those in my net search. I am hoping to find better code. I will dig into the material if nothing else comes along.

I have developed a complex-plain search scheme that works very well and can find all roots, including multiple ones. I am curious about the Lehmer-Schur algorithm as it basically seems to use a geometrical search *based on the description in Acton's book). The LS algorithm seems a bit elaborate.

Namir

Edited: 22 Aug 2010, 3:41 p.m.

I have developed a complex-plain search scheme

Is that a plain search scheme for the complex plane? Will you have a fancy search scheme for the complex plane later?

[Sorry, couldn't resist. GDR!]

Eric,

Actually I have developed a very simple scheme in Matlab that I will translate to other languages that do not support complex math. I also have a more efficient version for languages that DO support complex math. In latter version I perform polynomial deflation and detect complex conjugate roots (when the polynomial coefficients are all real). The more advanced version can detect multiple roots, while the simple one calculates distinct roots only. Both methods work with polynomial coefficients that are real and complex.

Namir

Edited: 23 Aug 2010, 11:26 a.m.

It doesn't seem like it would be too hard to implement on a TI89, which can do complex numeric integrals (I think the HP50 for some reason doesn't like it if integrals are complex).

I imagine you start with

a < real part < b

m < imaginary part < n

You believe there is one (or more) zeros in that rectangle.

Then you divide the rectangle into quarters, a to (a+b)/2 to b, and m to (m+n)/2 to n

Then you just numerically integrate f'/f around the rectangle. So, for the upper left one, it could be the sum of integrals

from a + m i to b + m i

then from b + m i to b + n i

then from b + n i to a + n i

and finally from a + n i back to a + m i

Then divide this number by 2 * pi * i. It should be an integer greater than or equal to zero. Because we're working with polynomials, there should be no poles, and the integer should equal the number of roots in the rectangle.

If the number of roots is greater than zero, quarter this new rectangle and try again.

If the number of roots equals zero, forget about this rectangle and move on to the next.

I'll see if I can whip up an example program on the TI89 later.

I can give one example that's simple enough to do by hand.

x-1=0.

Suppose the rectangle is bounded by

0<real<2
-1<imaginary<+1

In this case, f'/f = 1/(x-1)

The integral of this is log(x-1) = log|x-1| + theta * i

Along the bottom side of the rectangle, the integral will be

(log(sqr(2))- pi i / 4) - (log(sqr(2) - 3 pi i / 4)

= pi i / 2

Likewise, each other side contributes pi i / 2. So the total integral is 2 pi i. Divide by 2 pi i to get 1, the number of roots in the rectangle.

In this case, the next step, the quartering of the rectangle, will fail, because the root is exactly at the center. So all four rectangles will have zero for their winding number. But in the unlikely event that that happens, you should know that's where the root exactly is.

Actually, I'm not even sure if the TI89 likes to do numeric complex integrals. I can finesse it into doing some, but I'm not sure if it's really doing them numerically, or if it's doing them exactly, then converting the exact result to a number.

I guess you could always use your own numeric integration techniques.

For the previous example, if you used either the trapezoidal rule or the midpoint rectangle rule, with one rectangle / trapezoid on each side, then each side would give an approximate integral of i. Then the whole loop integral would be 4i, and the winding number would be 2/pi ~ 0.6. That's not exactly 1, but I guess your program could be prepared for approximations like that.

If it's of any help, in my copy of "Numerical Recipies 3rd Ed.", for a detailed description of Lehmer-Schur it refers me to "Newton Methods for Nonlinear Problems" (P. Deuflhard, pub. Springer 2004).

The Deuflhard books looks interesting (and so is it's price). I may bite the bullet and buy it.

I have been a fan of solving nonlinear equation since the mid 70s. I have always focused on non-polynomial functions and solved for real roots. Armed with Matlab, RPL, and other piece of software, solving for ALL roots of polynomial now seems possible. So I am cataloging iterative and non-iterative methods to solve for the roots of polynomials. The Matlab roots() function is a good tool to check my code, since it yields all of the roots of a polynomial. I think the roots() function is based on the Jenkins-Traub algorithm.

Namir

Crawl,

Can you rewrite your pseudo-code using indentations, loops, and decision making constructs?? It looks very interesting.

Of course a TI-89 example will be splendid since it will spell out things more clearer than pseudo-code.

Thanks!

Namir

Edited: 24 Aug 2010, 8:06 p.m.

Who needs matlab? :) Wave [link:http://www.wolframalpha.com/input/?i=x^4+-3i*x^2-84] the magic Wolfram Alpha wand [/link] and find the roots of

x^4+-3i*x^2-84

Link
http://www.wolframalpha.com/input/?i=x^4+-3i*x^2-84

Note: the forum link function seems not to be working.

Edited: 24 Aug 2010, 10:25 p.m.

I have used the WolframAlpha site to solve for the roots of polynomial. It works very well. My real intent is to program the various algorithms in several languages. I am starting with Matlab and will work my way into different languages, including RPL.

I was impressed to find that the HP-65 math pac had programs to solve cubic and quartic polynomials.

Namir

Edited: 25 Aug 2010, 10:32 a.m.

Okay. I think I'll try it for the HP50g after all, since it looks like with either calculator I'd have to numerically integrate.

My batteries are dead, though, and I'm charging rechargeables for it as I type (the first time I'll try that). So I won't be able to get around to it until tomorrow at the earliest.

Update: Looks like I might have to go back to the TI89. I can't figure out for the life of me how to do complex integration on the Hp50g (not counting "making my own routine" as an option -- I don't want to do that much work)

I figured out how to do it on the TI89, though.

Say you want to integrate e^(-x^2) from -1+i to 1+i. (You can be sure the calculator isn't doing this symbolically, since that would be impossible)

Just change it to

e^(-(x+i)^2), from -1 to 1. Then it works.

It seems as though the problem the TI89 has is with complex arguments in the LIMITS, so if you can avoid that (it doesn't matter if the integrand is complex), you can make it work.

I'm a little annoyed by this since I prefer programming on the HP50g, but whatever.

Thanks Crawl!

I am going to try to code some of cubic routines on the HP50G next ewek. I think I have come across some bugs in Matlab since two functions that use basically the same equations (expressed somewhat differently) give different results IN SOME CASES!!!

I would like to see if the 50G can do the job AND does not show similar bugs.

Namir

If you're interested in cubic equations, I posted a solution to the quartic equation (and the cubic) for the HP50g here. It works for symbolic or numeric input. I believe it has no bugs, but if you could find one I'd be interested.

Some day I might enter the simplified (from Mathematica) solution to the quartic into the HP50g.

Okay, I have to go back to the HP for the other thing. I figured out how to do complex integrals on it. Not only do the limits need to be real, the integrand does, too. But if you put a complex integrand in the RE() or IM() functions, it'll work.

On the other hand, the roadblock I was running into on the TI was that f(x+i)->g(x) doesn't work. Let's say f(x) = 1/(x-1), like in the earlier example. Then obviously what you want is g(x)=1/(x+i-1). However, what you actually get is g(x)=f(x+i). In other works, it doesn't evaluate the function f before making the assignment to g. You obviously could do that manually in normal calculator mode, but since program mode doesn't work the same as calculator mode (that's the advantage to programming in RPL), I don't know a workaround.

And ... it's done! And I checked it with at least one example.

You can download the programs from here. WIND calculates the Winding Number of the curve on the given rectangle for the given function. LEHMER is the main program, and each time it's run it quarters the previous rectangle.

First, you need to define F(X) to be f'(x)/f(x), with f(x) the function you're trying to zero.

In our example, f(x)=x^4-2*x^2+10. Then F(X)=(4x^3-4x)/(x^4-2*x^2+10)

I so happen to know that the exact solution to this equation is x = sqrt(1+3i)~1.442+1.0397i. Let's say that we know the solution is between x = 1+i and 2 + 2i.

So, you put those numbers on the stack. I believe that I intended the lower left corner of the rectangle to be on the stack first, but I also think I tried to make the program work even if the numbers were reversed.

Then just press LEHMER to iterate.

This is what it looks like:

The program assumes that there is one, and only one, root in the given rectangle, and also no poles. The case for which multiple zeroes are in the rectangle can be left as an exercise for the reader.

One thing you'll note about this program (really the method's fault more than my programs): It is INCREDIBLY slow! And somewhat error prone. I would definitely not recommend using it on a real HP50g. I did, and the first iteration took about a half an hour. Even on an emulator, it seems as when you approach the root, the integrations become slower. For the GIF I posted, even on an emulator, the last iteration took quite a while. So be forewarned.

Also, the first time I tested it, I checked to make sure the calculated winding number was within 0.1 of 1. It ran for a while, and terminated with "ERROR". I single stepped through the program, and the calculated winding number was something like 1.15. So I changed the margin of error to 0.3. It still should have to be closer to 1 than 0.

I don't know if this method has much use. It seems to me to be much worse than bisection for real numbers -- and bisection is a pretty poor algorithm. At least bisection is robust -- if there's a root in the interval, you're almost guaranteed to find it. This method has all kinds of problems.

Edited: 27 Aug 2010, 11:32 p.m.