# Root finding for holomorphic functions

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]


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}

• Unfortunately, FindRoot[Zeta'[s], {s, 2.5 + 25.5 I}] produces "FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations." and {s -> 146.243 + 23.8385 I} Sep 20 at 18:36
• So the convention for the complex rectangle is to name the lower left and upper right corners. The documentations says "FindRoot[lhs==rhs,{x,xstart,xmin,xmax}] searches for a solution, stopping the search if x ever gets outside the range xmin to xmax." May I assume this holds for complex x as well? Sep 20 at 18:46
• @stopple I'm not sure, but it appears to work this way. Sep 20 at 18:54
• FindRoot[Zeta'[s], {s, 2.20 + 19* I, 2 + 18 I, 3 + 24 I}] produces "FindRoot::reged: The point {3. +19.0696 I} is at the edge of the search region {2. +18. I,3. +24. I} in coordinate 1 and the computed search direction points outside the region." and {s -> 3. + 19.0696 I}. Sep 21 at 3:50
• D[Zeta[s] ,s]/. {s -> 3. + 19.069576455342865* I} results in -0.0591024 + 0.101857 I. Sep 21 at 19:16

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.)

• Unfortunately, Reduce[D[Zeta[s], {s, 1}] == 0 && Re[s] < 2 && Re[s] > -2 && Im[s] > 0 && Abs[s] > 10^5 && Abs[s] < 10^5 + 100, s] returns "Reduce::nsmet: This system cannot be solved with the methods available to Reduce.". Sep 20 at 17:05
• This method takes 158 seconds to fail to find even the first non-real zero, at height about Im[s]=23. Sep 20 at 17:18
• Hmm, indeed: Reduce[D[Zeta[s], {s, 1}] == 0 && 2 < Re[s] < 3 && 23 < Im[s] < 24, s] comes back unevaluated for me. Sep 20 at 17:43
• Just to compare. The command of Maple 2023 RootFinding:-Analytic(diff(Zeta(z), z) = 0, z = 300*I .. 1 + 300*I + 100*I) produces 23 roots 0.8723871165 + 349.8295184*I, 0.9943529705 + 324.4687500*I,..., 0.6857226375 + 395.9738587*I, 0.9041376510 + 392.9718602*I, 0.6691537285 + 391.7936220*I in few minutes. Sep 20 at 20:03

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.

• This is really helpful; thank you. Sep 21 at 20:26