Root finding for holomorphic functions

Question

Mathematica's FindRoot has a lot of options, but the documentation is not helpful for what I would have expected is an important case, numerically finding complex roots of holomorphic functions. Via the Argument Principle, one can often identify a rectangle in which the holomorphic function has an exactly one zero. What is the best choice of options to return exactly that zero and no other?

It should not matter, but I'm interested in finding zeros of the derivative of the Riemann Zeta function, high in the critical strip. I'm calculating the derivative via:

<< NumericalCalculus`

ZetaPrime[s0_?NumericQ] := ND[Zeta[s], {s, 1}, s0, Method -> NIntegrate]

Answer 1

This seems to work:

d = 1/2;
re = {0, 5};
im = {10, 50};
Flatten[Table[{i + I j, i + d + I (j + d) }, 
   Evaluate@{i, Sequence @@ (re - {0, d}), d}, 
   Evaluate@{j, Sequence @@ (im - {0, d}), d}], 1];
Quiet@Check[
     FindRoot[
      Zeta'[s], {s, Mean[#], Sequence @@ #}], {}] & /@ % // Flatten
Clear[d, re, im]

(* {s -> 0.964686 + 48.8472 I, s -> 1.2865 + 31.7083 I, 
 s -> 1.38276 + 42.291 I, s -> 2.46316 + 23.2983 I, 
 s -> 2.30757 + 38.49 I} *)

I sampled region form 0 to 5 in real direction and from 10 to 50 in imaginary direction by d=1/2. And it found all five zeros in that region.

If there is a region where the zeros are more dense maybe smaller sampling value d than1/2 should be used.

Answer 2

FindRoot has a way of passing in a bounding rectangle:

FindRoot[f[s], {s, s0, smin, smax}]

If we know such a rectangle, we can use its centroid as the initial guess:

FindRoot[Zeta'[s], {s, 2.5 + 23.5 I, 2 + 23 I, 3 + 24 I}]

{s -> 2.46316 + 23.2983 I}

Answer 3

Somewhat surprisingly, roots of a holomorphic function over a bounded domain can be found exactly in terms of Root objects.

sol = Reduce[2 Sin[Exp[x]] - Cos[Pi x] == 1/3 && Abs[x] < 1, x]

   x == Root[{1/3 + Cos[\[Pi] #1] - 2 Sin[E^#1] &, -0.108292008603703700 - 0.243578228378785204 I}] 
|| x == Root[{1/3 + Cos[\[Pi] #1] - 2 Sin[E^#1] &, -0.108292008603703700 + 0.243578228378785204 I}]

N[sol, 30]

   x == -0.108292008603703700225204957273 - 0.243578228378785204495911435207 I 
|| x == -0.108292008603703700225204957273 + 0.243578228378785204495911435207 I

Here's a blog post about it: https://blog.wolfram.com/2008/12/18/mathematica-7-johannes-kepler-and-transcendental-roots/

---

A relevant excerpt from that post:

So how do we solve transcendental equations?

We start by looking at the class of holomorphic functions—essentially polynomials of infinite degree. A holomorphic function will often have infinitely many roots. But a key fact is that the roots are always countable—and, more importantly, there can only be finitely many roots in a given closed and bounded region.

Now we get to use some complex analysis. Given a region, we can tell how many roots are inside it either by directly computing the winding number around its boundary, or by using numerical integration and Cauchy’s theorem.

Once we have a count of the roots in a region we have a variety of methods for working out their rough specific locations. (Some methods are based on computational geometry and curve-curve intersections, some on complex interval arithmetic, and some on generalized eigenvalues derived from Cauchy’s theorem.) Given rough locations for roots, we then use Newton-like iterative methods to home in on the roots. And finally, we use Mathematica’s interval arithmetic to prove that there can only be one root inside each small region we’ve identified.

(There is some subtlety here. Proving zero equivalence is in general undecidable—and the way this shows up is that in pathological cases it can in principle require unboundedly much precision to distinguish multiple roots from closely spaced single roots.)