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.