The problems with the method seems to arise from two sources, the objective function (and its components Re, Im, HarmonicNumber) and the working precision (MachinePrecision).
The objective function
The Jacobian: The functions Re and Im do not have derivatives. Inspecting the Jacobian at {56.9935, 171.987} shows unevaluated derivatives:
D[{Re[HarmonicNumber[n, x + I y]], Im[HarmonicNumber[n, x + I y]]}, {{x, y}}] /.
Thread[{x, y} -> {56.9935, 171.987}]
((* {{(-4.7632*10^-18 - 8.08894*10^-19 I) Derivative[1][Re][1. + 1.16699*10^-18 I],... *)
The derivatives of HarmonicNumber evaluate to a complex number but the derivatives of Re and Im do not. In such a case, I believe FindRoot will use a finite difference method for computing the Jacobian. Insufficient working precision can lead to issues in computing the Jacobian, which is not the hint provided by FindRoot::jsing message. In particular, function evaluations that mathematically should produce distinct numbers might be equal at MachinePrecision and cancel out.
Precision: Mathematically, HarmonicNumber[n, z] and Sum[1/k^z, {k, n}] are equivalent. One would like to think that HarmonicNumber would be at least as precise and fast as the sum. This may be true over some domains. For small values of n (certainly up to 200 and more), the sum is more precise. In fact, over the domain of this problem (n <= 200, z in the complex rectangle from -30 + 0 I to 1 + 6000 I), the sum is much more precise. For n = 20 and 1000 random, arbitrary-precision complex numbers with precision $MachinePrecision (15.95), the mean loss of precision for the sum was 4.3 digits and for HarmonicSum, it was 14.1 digits. For n = 20 and 1000 random, arbitrary-precision complex numbers with precision 100, the mean loss of precision for the sum was 4.3 digits and for HarmonicSum, it was 26.6 digits.
Speed: The sum is also considerable faster than HarmonicNumber for small values of n. At about n == 1800, the speed of the sum and HarmonicNumber are about equal; for larger values of n, HarmonicNumber is faster.
Compilability: HarmonicNumber is not compilable. This is not a serious problem, but the FindRoots2D code compiles the objective functions.
Working precision
MachinePrecision is good enough for many cases, but numerical errors can accumulate and (1) cause a Jacobian computed with finite differences to vanish erroneously or (2) confuse the root-finding algorithm with erratic function evaluations. These seem to be happening in this case. The first may be the reason for the FindRoot::jsing error; the second may be the reason for the FindRoot::lstol warning. For a function with isolated critical points, a FindRoot::jsing ought to be rare. A FindRoot::lstol warning can happen when FindRoot gets stuck near a local extremum, but I was not able to find such an example in the OP's problem after a (very) brief search.
Some improvements
Increase WorkingPrecision
Setting WorkingPrecision -> 100 took care of the FindRoot::jsing errors. It did not take care of all the FindRoot::lstol warnings. However, it is much slower to use arbitrary-precision arithmetic on all solutions.
Use the explicit sum instead HarmonicNumber
As described above, using an explicit sum (or Plus) expression results in greater accuracy and speed. Below myHN is my version of HarmonicNumber. If real and imaginary parts are desired, then myHNReIm can be used. It can be compiled, which will result in much faster computations at MachinePrecision. (For a specific integer n, each definition generates the sum to be computed and memoizes it. In this way the Sum is expanded only once.)
myHN[n_Integer, z0_] := Block[{z}, (* complex-valued *)
myHN[n, z_] := Evaluate@Sum[1/k^z, {k, n}];
myHN[n, z0]
];
myHNReIm[n_Integer, x0_, y0_] := (* real and imaginary parts *)
Block[{x, y, k},
With[{reHN = Simplify[Re[1/k^(x + I y)] // ComplexExpand,
x ∈ Reals && y ∈ Reals && k > 0],
imHN = Simplify[Im[1/k^(x + I y)] // ComplexExpand,
x ∈ Reals && y ∈ Reals && k > 0]},
myHNReIm[n, x_, y_] := Evaluate@{Sum[reHN, {k, n}], Sum[imHN, {k, n}]};
myHNReIm[n, x0, y0]
]];
myHNReImCx[n_Integer] := (* real part, compiled *)
myHNReImCx[n] =
Compile @@ Hold[{x, y}, Evaluate@First@myHNReIm[n, x, y]];
myHNReImCy[n_Integer] := (* imaginary part, compiled *)
myHNReImCy[n] =
Compile @@ Hold[{x, y}, Evaluate@Last@myHNReIm[n, x, y]];
myHNReImC[n_Integer, x_?NumericQ, y_?NumericQ] := (* interface to compiled parts *)
{myHNReImCx[20][x, y], myHNReImCy[20][x, y]};
Use the complex-valued function
Using myHN, which returns complex values is faster than using the individual components (such as in myHNReIm). This is true for HarmonicNumber, too. First of all, FindRoot can deal with such functions. For either myHN or HarmonicNumber, the Jacobian can be computed symbolically, which helps FindRoot to be efficient. And the way the real and imaginary parts are used in the call to FindRoots2D, the harmonic number has to be computed twice for each input.
Use Check to catch troublesome seeds
One can use Check to catch error/warning messages and re-execute FindRoot with adjusted settings. The following catches singular Jacobians, line-step tolerance warnings, and nonconvergence after the maximum number of interations and restarts FindRoot with a higher working precision. If messages are generated on the second try, an empty solution {} is returned on the assumption that a zero is not going to be found near the seed (starting point) x0. One could instead store them in some other variable and inspect them later. They might correspond to an actual root.
Check[{sol = FindRoot[f, {x, x0}]
Check[FindRoot[f, {x, x /. sol}, WorkingPrecision -> 100, AccuracyGoal -> 8],
{}],
{FindRoot::jsing, FindRoot::lstol, FindRoot::cvmit}]
Whether to restart FindRoot where it stopped x /. sol or from the initial starting point x0 is open to debate. If the solution sol is closer to the true zero, it will save some time. If round-off error has moved the solution far away from the true zero, the zero may be missed. In the OP's case, it seems best to start at the original seed (see below).
Use MeshFunctions to look for seed points
The following an arguable improvement. One can use ContourPlot in the form
ContourPlot[f1, {x, 0, 1}, {y, 0, 1}, MeshFunctions -> {f2}, Mesh -> {{0}}]
to find the seed points, which I've used before (see also 1, 2). Here f2 needs to be a pure function; it also ought to be continuous for good behavior. ContourPlot takes time and adding mesh functions increases the time it takes. It actually will search for the point on the contour line where f2 == 0.
Use the real and imaginary parts both ways in ContourPlot
Using one ContourPlot of the real part and searching every contour line for a change in sign of the imaginary part, as in FindRoots2D, is generally a good way to find zeros. One could also switch the roles of the real and imaginary parts. It turns out, more zeros are discovered if you do both. Since ContourPlot takes time, both plots could be computed in parallel, which won't take any more real time, if you have at least two cores.
Avoid Union
Unless you need the zeros sorted, Union tends to be slow. Here is a way to delete repeated roots that has two advantages, it is a bit faster and you can pick the root that gives the smallest value of HarmonicNumber. While HarmonicNumber is rather oscillatory, this probably gives the root closest to the actual zero. The function findClusters gathers, for each data point, the data points within a certain radius and then deletes the duplicate clusters:
findclusters[data_, r_: 10^-12] := With[{nf = Nearest[data]},
DeleteDuplicates[Sort[nf[#, {All, (* use Length[data] instead of All for pre-V10 *)
r}]] & /@ data]
];
One can then select the member of each cluster that produces the smallest function value like this:
reducedzeros = #[[First@Ordering[Abs @ myHN[20, #]]]] & /@ findclusters[zeros, 10^-6]
Overall this was about 6 times faster than Union on 27000 zeros.
Stopping runaway roots
If we use FindRoot like this
FindRoot[f, {z, seed, -30 + 0. I, 1 + 6000 I}]
we can catch the message FindRoot::reged with Check. This message indicates that the search for the root has left the region. If the region is slightly larger than {-30 + 0. I, 1 + 6000 I}, then perhaps such a seed can be safely discarded. It saves a lot of time, since FindRoot will pursue a root well beyond the region of interest.
Parallelize
There are a few things that can be parallelized, and since multi-core machines are now common, it makes sense to do so. This is particularly so for the FindRoot calls on the thousands of seeds.
Some code
I'll separate FindRoots2D into two parts, finding the initial seeds and finding the zeros.
Clear[findSeeds];
findSeeds[f_, {x_, a_, b_}, {y_, c_, d_}, opts : OptionsPattern[ContourPlot]] :=
Module[{plotf, mf1, mf2},
plotf = Function @@ Hold[{x, y}, Evaluate@f];
If[OptionValue[MeshFunctions] === {},
{mf1, mf2} = {
With[{mf =
Compile[{x0, y0}, Evaluate[Re[f] /. {x -> x0, y -> y0}]]},
mf[##] &],
With[{mf =
Compile[{x0, y0}, Evaluate[Im[f] /. {x -> x0, y -> y0}]]},
mf[##] &]
},
{mf1, mf2} = OptionValue[MeshFunctions]
];
Join @@ Parallelize[
{Cases[Normal@#, Point[p_] :> p, Infinity] &[
ContourPlot[
Re @ f == 0,
{x, a - (b - a)/97, b + (b - a)/103}, {y, c - (d - c)/98, d + (d - c)/102},
ContourStyle -> None, MeshFunctions -> {mf2}, Mesh -> {{0}}, opts]
],
Cases[Normal@#, Point[p_] :> p, Infinity] &[
ContourPlot[
Im @ f == 0,
{x, a - (b - a)/97, b + (b - a)/103}, {y, c - (d - c)/98, d + (d - c)/102},
ContourStyle -> None, MeshFunctions -> {mf1}, Mesh -> {{0}}, opts]
]}
]
]
The tolerance to reject a seed as a nonroot is set at 10^-6. This could be changed or made an option. It appears in two places.
Clear[myFRC];
myFRC[f_, {z_, a_, b_}, seeds_, opts___] := Module[{z1, z2},
z1 = Re[a - (b - a)/97] + Im[a - (b - a)/98] I;
z2 = Re[b + (b - a)/103] + Im[b + (b - a)/102] I;
z /. Flatten[
ParallelMap[
Module[{sol},
Quiet[ (* did not suppress FindRoot::jsing *)
Check[ (* for FindRoot::reged : no solution *)
Check[{ (* for FindRoot::jsing, FindRoot::lstol, FindRoot::cvmit : restart *)
sol = FindRoot[f, {z, #1, z1, z2},
Evaluate[FilterRules[{opts}, Options[FindRoot]]]]},
Check[ (* for any message : no solution *)
If[Abs[f /. sol] < 10^-6, (* tolerance 10^-6 for nonzero *)
{sol},
If[Abs[f /. #] < 10^-6, {#}, {}] &@
FindRoot[f, {z, SetPrecision[#1, 100], z1, z2},
WorkingPrecision -> 100, AccuracyGoal -> 8,
Evaluate[FilterRules[{opts}, Options[FindRoot]]]]],
{}],
{FindRoot::jsing, FindRoot::lstol, FindRoot::cvmit}],
{}, {FindRoot::reged}],
{FindRoot::reged, FindRoot::lstol, FindRoot::cvmit}
]
] &,
seeds
],
1]
]
Example
The steps in this example could replace the body of the code for FindRoots2D, if desired.
I used Pick with UnitStep instead of Select. In this case it is about ten times faster. Packing the array (Developer`ToPackedArray@N@... makes only a little difference, on my machine); there is also a small time savings from reducing the arbitrary precision zeros to MachinePrecision, whether the list is packed or not.
n = 20;
seeds = findSeeds[myHN[n, x + I y], {x, -30, 1}, {y, 0, 6000},
PlotPoints -> 25,
MeshFunctions -> {myHNReImCx[n][#1, #2] &,
myHNReImCy[n][#1, #2] &}]; // AbsoluteTiming
allfoundzeros =
Quiet @ myFRC[myHN[n, z], {z, -30, 1 + 6000 I},
seeds.{1, I}]; // AbsoluteTiming
uniqueszeros = #[[First @ Ordering[Abs @ myHN[20, #]]]] & /@
findclusters[allfoundzeros, 10^-6]; // AbsoluteTiming
zeros = Pick[
uniqueszeros,
With[{zz = Developer`ToPackedArray @ N @ uniqueszeros},
UnitStep[(#1 - (-30)) (1 - #1) #2 (6000 - #2)] &[Re[zz], Im[zz], 1]],
1]; // AbsoluteTiming
Length[zeros]
zeroplot = Graphics[{PointSize[Tiny], Red, Point @ Transpose@{Re[zeros], Im[zeros]}},
PlotRange -> {{-30, 1}, {0, 6000}}, AspectRatio -> 1, Frame -> True]
(*
{8.560900, Null}
{32.787730, Null}
{1.368914, Null}
{0.011737, Null}
{0.000582, Null}
2429
*)
Verification:
myHN[n, zeros] // Abs // Max
(* 5.11995*10^-7 *)
Caveat
Being sure you have found all the zeros seems problematic. Changing the PlotPoints option gave different numbers of seeds and zeros. There seem to be about ten times as many seeds generated as zeros, which strikes me as a very poor ratio. Increasing the WorkingPrecision on ContourPlot did not change the result.
