Finding all complex polynomial roots in a specified range efficiently

Question

I need to find the roots of a rational polynomial that are near i. In the following code, I try that two different ways. First, I use a constraint to only find roots in the right region. Second, I find all the roots and discard the ones not in the right region. Surprisingly, the second, wasteful method is much faster. Why?

z[0, c_] := c
z[n_, c_] := z[n - 1, c]^2 + c;
r = 1/10;
prec = 10;
expr = Expand[z[10, c] - z[6, c]];
Print[Timing[aa = c /. N[Solve[{expr == 0, Abs[c - I] < r}, c], prec];]];
Print[Timing[bb = Select[c /. N[Solve[expr == 0, c], prec], Abs[# - I] <= r &];]];

(* {4365.47,Null} {176.39,Null} *)

Clearly Solve[] is doing something terribly inefficient in its handling of constraints. I know I could use FindRoots[] to find roots one at a time and divide them out of the polynomial, then I wouldn't know when I'd found all the roots in the region and could stop looking. But when I tried that, I ran into this problem.

So how can I efficiently find all the roots within radius r of i?

And what the heck is going with Solve[]?

Answer 1

Here is an entirely different numeric method, from this Usenet post. The idea is to use the Cauchy Integral formula for root counting, isolation, and refinement. I think the original code had some problems so below is modified (and thus probably has different problems..)

SetOptions[NIntegrate, WorkingPrecision -> 150, AccuracyGoal -> 8, 
  PrecisionGoal -> 8];

eps = 1.3*10^(-6);
isolateRootsInRectangle[func_, a_, b_, c_, d_, z_, eps_] := 
 Module[{numroots, avgx = (a + b)/2, avgy = (c + d)/2, r1, r2, r3, r4,
    len = 0, allrts}, 
  numroots = 
   Round[rectangleIntegrate[func, a - eps, b + eps, c - eps, d + eps, 
     z]];
  If[numroots <= 0, Return[{}]];
  If[numroots == 1, Return[{{a, b, c, d}}]];
  r1 = isolateRootsInRectangle[func, a, avgx, c, avgy, z, eps];
  r2 = isolateRootsInRectangle[func, a, avgx, avgy, d, z, eps];
  r3 = isolateRootsInRectangle[func, avgx, b, c, avgy, z, eps];
  r4 = isolateRootsInRectangle[func, avgx, b, avgy, d, z, eps];
  Join[r1, r2, r3, r4]]

I now find intervals for roots that lie in a box of half-length 1/10 centered at I.

r = 1/10;
Timing[intervals = 
   isolateRootsInRectangle[expr, -r, r, 1 - r, 1 + r, c, eps];]

(* Out[567]= {27.433912, Null} *)

This may have duplicates due to use of slight overlaps on region, and I did not attempt to remove them.

Now we refine them. Could use FindRoot with start points at interval centers, but this is more reliable (and a FindRoot based on these values might be better still).

Timing[rts = refineRoots[expr, intervals, c, eps];]

(* Out[573]= {6.487449, Null} *)

As for working with a circular rather than square region, there is more code at the link above. But that involves extra fuss and it's probably simpler just to remove those that are inside the box but outside the circle in question.

Answer 2

Solve also finds all complex roots and selects those that satisfy the inequality condition. Most of the timing difference is due to Solve evaluating the condition with exact values of roots and your code evaluating the condition with 10-digit approximations of roots. Exact evaluation of inequalities computes the expressions involved to at least 54 bits of precision. In this example computing the extra digits of precison is costly.

In[3]:= (rts=Solve[expr == 0, c]);//Timing
Out[3]= {126.159, Null}

In[4]:= (sol=Select[c /. rts, Abs[# - I] <= r &]);//Timing
Out[4]= {1229.03, Null}