[update 2] Urgghh - the time-consumption really stems only from the construction of the h-order polynomial. The time for finding the roots (only 10 to 20 times Newton-iteration because of my nice start-intervals) is then completely neglectable.
So the effectiveness-aspect of my question has been solved, as it depends on the false premise of a costly Newton-iteration.
Only the more qualitative question of the existence of some procedure for finding only the few real roots of a polynomial remains open (but may be have the answer to the negative).

I have a family of polynomials of increasing order (see my discussion at the earlier MSE-thread), and have the strong hypothese, that each of that polynomials has either 1 or 2 real roots, and always exactly one positive - is then there an algorithm to find just that one real root?

I looked at "polsturm" in Pari/GP which allows to determine the number of real roots, but not also their value - I still have to use "polroots" which computes all complex roots.

Well, I have also a hypothese of the rough value of that root, so I tried also the inbuilt solver with a reasonable small interval.

But with polynomials of order 200 both methods need fairly time, and with order of 1000 the solver needed about 10 minutes...
Perhaps there is some Horner-like scheme? So my question is:

Q: Is there an efficient algorithm for finding the real roots of a high.order polynomial if I know there are at most 2?

[update] Using the solver-method I get that time-consumption $$ \small \begin{array} {r|r} \text{order} & \text{msec} \\ \hline \\ 62 & 31 \\ 125 & 296 \\ 250 & 3292 \\ 500 & 39967 \\ 1000 & \gt 400000 \end{array}$$

This is the relevant part of the code in Pari/GP

glpc =[0]   \\ introduce a global accessible vector holding the coefficients
            \\ of the current polynomial

{evalpoly(x)= local(f,su);  \\ the Horner-scheme evaluator of the polynomial
    su=0;f=1;             \\ works on the global vector glpc[]

{mypolroot(h)= local(xl,xu,res);  \\ creates the polynomial of order h
                                 \\ finds the real root in the interval
                                 \\ depending on h by my hypothese
  glpc = polcoeffs (makepoly(h));  \\ put the coeffs of the h'th polynomial 
                                  \\ into glpc
   xl=floor(h*gm/2); xu = ceil( h*gm/2+5); \\ "magic" lower and upper bound
                                           \\ for the positive real root
   res = solve(x= xl,xu , evalpoly(x) );

You say you have a good guess for the location of the root -- you've tried Newton's method?

If Newton's method doesn't converge, i.e. your guess is too inaccurate, you can take advantage of the knowledge that there is only a single positive root to bracket the root inside the interval $[a,b]$, with $f(a)f(b) < 0$. Worst case you can use $a=0$ and find a $b$ by successive doubling. Then you can run the Ridders' method. This method is guaranteed to converge, and at worst linearly, though it converges much faster in practice once you're near the root.

EDIT: Also, I feel like I should mention that I've used the Jenkins-Traub algorithm (C code here) in my research projects in the past, as the first (and almost always last) stop for finding polynomial roots. But finding all roots seems overkill when you're only looking for one and can bracket it effectively.

| cite | improve this answer | |
  • $\begingroup$ Yes, thanks for the hint! Well, I'm aware of the newton-method, and I think the inbuilt-solver (Pari/GP) in fact uses this method. The time consuming aspect is surely the repeated evaluation of a 500- or 1000-coefficients polynomial (I've optimized this already with a Horner-schema). I thought, if there is a standard method for finding only the real roots, then this might be more efficient... $\endgroup$ – Gottfried Helms Sep 25 '13 at 13:20
  • $\begingroup$ Evaluating a 1000-coefficient polynomial, especially in Horner form, should be child's play for a modern computer. Can you get more detailed information from Pari, i.e. the residual at every iterations and the time needed to evaluate $f$ (and its derivatives) at every iteration? $\endgroup$ – user7530 Sep 25 '13 at 13:25
  • $\begingroup$ I've added the Pari/GP-code in my question. Also I'll try whether a self-made Newton-iteration makes things faster... $\endgroup$ – Gottfried Helms Sep 25 '13 at 13:31
  • 2
    $\begingroup$ I expect a basic Newton-Raphson solver written in a compiled language will be thousands of times faster than Pari/GP. I think Pari/GP is the bottleneck, here. $\endgroup$ – bubba Sep 25 '13 at 13:58
  • 1
    $\begingroup$ Well, I've seen that I had a complete false impression about the consumtion of time of the two parts of my algorithm. The Newton-iteration actually is in the milliseconds, and in fact the timing problem occurs in the construction of the polynomials. Sorry for the noise... $\endgroup$ – Gottfried Helms Sep 25 '13 at 15:23

The Chebfun system (see here) takes only a fraction of a second to calculate all the real roots of a polynomial whose degree in the thousands. The code is written in Matlab, so presumably it would be even faster if it were written in a compiled language.

The root-finding technique is based on computing the eigenvalues of companion matrices. However, a naive implementation gives an algorithm that is $O(n^3)$ in the degree of the polynomial. To improve matters, Chebfun recursively approximates high-degree polynomials by low degree ones. Here "low" means less than 100. See Trefethen's book for details.

I don't know if there are algorithms that take advantage of the fact that the polynomial has only 1 or 2 real roots. But the Chebfun one seems like it might be fast enough for your needs, anyway. A Newton-Raphson solver with good starting points would be even faster.

Also, I agree with user7530's assessment -- evaluating the polynomial is not likely to be the bottleneck on a modern computer, or even on a modern cell-phone or wrist-watch :-)

| cite | improve this answer | |
  • 1
    $\begingroup$ Wow. I had the idea of doing something using eigenvalues - but didn't follow that idea because of the $O(n^3)$ -complexity which should be my barrier (being naive, as I am...). But that examples at Chebfun look real great. Thanks, I'll look at it later today! $\endgroup$ – Gottfried Helms Sep 25 '13 at 13:57

If you are absolutely sure that there is only one positive real root, then a combination of Newtons method and interval nesting a la bisection method or regula falsi should work. Pure Newton may fail by entering an infinite stable loop.

First determine the outer root radius R by one of the various formulas and take [0,R] as initial interval. Compute the values. Apply Newtons method as long as it stays in the interval, shrink the interval at each step by replacing the boundary with the same sign of the function value. If Newton fails to stay in the interval, compute the root of the secant as next iteration point.

One could enrich this by a counter of how long the same side of the interval gets changed to insert a bisection step every 5 or 10 iterations to really reduce the interval. Feel inspired by the Illinois variation of regula falsi.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.