7

I'm currently using Reduce to find zeros of a function in a region on a complex plane:

Reduce[q[s]==0 && 0 < Re[s] < 1 && 0 < Im[s] < 50]

q[s] is a large Dirichlet polynomial:

q[s_] := Sum[a[n]/n^s, {n, 1, 70}]

(the sum in my case is explicit with explicit coefficients, but it does not matter here, I believe)

This function is entire, and so I expected Mathematica to find its zeros rather quickly - as far as I know, there exist special algorithms for finding zeros of a holomorphic function in a region. However, Reduce is very slow - on my machine it could easily take 2-3 hours or even more to locate zeros of one polynomial.

Does Mathematica have some special built-in algorithm for this? For example, Maple has RootFinding[Analytic] function which does indeed find zeros of my functions relatively quickly (usually in less than a minute or two), but it fails for unknown reason when the degree of a polynomial is above 60-70 and when the region is larger than approximately 1x50 (and anyway I prefer using Mathematica for my computations). Numeric approximation of roots is perfectly acceptable for me, absolute accuracy is unnecessary - however, if I use NSolve it works even slower :( Maybe there is a way to tweak the performance of Reduce somehow? I couldn't find any relevant options in its documentation.

I believe it is possible to run FindRoot over some grid of points in the region but I don't want to do it because it is easy to miss some zeros this way.

Here is an example of the Dirichlet polynomials I'm working with.

Vladimir Matveev
  • 245
  • 1
  • 6
  • Is the variable s? If so, this seems like it might be difficult since it's not polynomial. – Daniel Lichtblau Apr 20 '15 at 15:01
  • Yes, the variable is s. I searched the internet and found some papers which suggest methods of numerical approximations of zeros of holomorphic functions - and it seems that one of them is used in Maple. I wonder if Mathematica also has something like it, maybe in form of parameters to generic functions like with NIntegrate and its methods of integration. – Vladimir Matveev Apr 20 '15 at 15:06
  • 1
    I think numeric integration is sometimes used by Reduce although I am not absolutely certain it is used in this example. That of course does not address the speed issue under consideration. One thing I will suggest is that you post an explicit example (no symbolic a[n], that is). That way readers might have something to test against. – Daniel Lichtblau Apr 20 '15 at 15:09
  • For a numerical solution, you could try one of the answers to Updating Wagon's FindAllCrossings2D function, using the real and imaginary parts of your function as input. – Jens Apr 20 '15 at 16:04
  • @DanielLichtblau, thanks, I should have done it since the beginning but all this time Mathematica was computing the next portion of zeros and I didn't want to interrupt it. I've added a link to an example function. – Vladimir Matveev Apr 20 '15 at 16:58
  • @Jens, yes, I saw that some time ago and unless I find something more appropriate I think I will go that way. Still hoping for something less "hacky". – Vladimir Matveev Apr 20 '15 at 17:01

1 Answers1

11

It turns out that Reduce finds candidate solutions relatively quickly and spends the vast majority of time proving correctness and completeness of the result. NSolve didn't have its own code for handling such problems, and was ending up using the same code as Reduce, finding symbolic solutions, and then numericizing them. I have implemented an NSolve version of the code, which only verifies candidate solutions numerically. It will be a part of the next Mathematica release. Unfortunately, I don't see an easy way of making the current Mathematica skip the symbolic solution verification... You can try the following code which directly calls the internal function finding candidate solutions.

nsol[f_, {x_, a_, b_}] :=
   Module[{pf, cand, rts},
      pf=Evaluate[f/.x->#]&;
      Reduce`AnalyticRootIsolation;
      cand=System`TRootsDump`GuessRoots[pf, {a, b}];
      rts=Quiet[x/.(FindRoot[f, {x, #}]&/@cand)];
      rts=Select[rts, 
            Chop[pf[#]]==0 && Re[a]<=Re[#]<=Re[b] && Im[a]<=Im[#]<=Im[b]&]; 
      Union[rts, SameTest->(Chop[#1-#2]==0&)]]

In[3]:= (rts=nsol[f, {s, 0, 1+50 I}]);//Timing
Out[3]= {1.471776, Null}

In[4]:= Length[rts]
Out[4]= 11
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Adam Strzebonski
  • 3,510
  • 23
  • 17
  • 1
    Wow, that's very interesting. I did notice that Reduce said something about proving correctness. For now I ended up using the approach with analyzing ContourPlot as suggested in comments above, but this does answer the question. Looking forward for the next version then :) Thanks! – Vladimir Matveev Apr 20 '15 at 22:22
  • Does Mathematica 10.4 implement the process described in your answer? – bbgodfrey May 09 '16 at 04:13
  • 1
    System`TRootsDump`GuessRoots seems now to require a third argument, True/False, that indicates whether to try harder to find roots. – Michael E2 Oct 08 '17 at 17:23