This is a follow-up to my question Root finding for holomorphic functions. I am trying to compute $10^6$ zeros of the derivative of the Riemann zeta function $\zeta(s)$ in the critical strip near height $t=5\cdot 10^6$. I've implemented the elegant solution of @azerbajdzan, but I need a slightly more sophisticated algorithm.
By what is known as 'the argument principle' for a holomorphic function like $\zeta^\prime(s)$, the change in Arg[$\zeta^\prime(s)$] around the boundary of a box, divided by $2\pi$, gives the number of zeros of $\zeta^\prime(s)$ inside the box. It is efficient to check this first before calling FindRoot: although each corner of the box requires a function evaluation, if we save the values, we can use them on an additional three boxes. Thus (neglecting boundary), we can determine if there's a zero in the box with about a single function evaluation. It will also be more accurate: if FindRoot returns an error FindRoot::reged because the algorithm converges to a root outside the box, I know not to ignore it but to recompute with a smaller step size.
Now, finally, to my question. I'm trying to write a recursive function to implement this idea, and I can't see what's wrong. Here's some code
ZetaPrime[s0_?NumericQ] :=
ZetaPrime[s0] =
N[-(100/21) (315 Zeta[s0] - 512 Zeta[1/800 + s0] +
224 Zeta[1/400 + s0] - 28 Zeta[1/200 + s0] + Zeta[1/100 + s0])];
This computes $\zeta^\prime(s)$ by Richardson extrapolation, just like ND in the NumericalCalculus package, for computing the arguments.
ZetaPrimeFine[
s0_?NumericQ] := -(100/21) (315 Zeta[s0] - 512 Zeta[1/800 + s0] +
224 Zeta[1/400 + s0] - 28 Zeta[1/200 + s0] + Zeta[1/100 + s0]);
The same but in symbolic form, not saved and without a specific number of digits, to pass to FindRoot.
argchange[a_, b_, zetaa_, zetab_] :=
Module[{argdif = Arg[zetab/zetaa], c, d},
If[Abs[argdif] < .75, argdif, c = (a + b)/2; d = ZetaPrime[c];
argchange[a, c, zetaa, d] + argchange[c, b, d, zetab]]];
This computes the change in Arg[$\zeta^\prime(s)$] from a to b.
SmartFindRoot[box_] :=
Module[{ll = box[[1]], ur = box[[2]],
lr = Re[box[[2]]] + I*Im[box[[1]]],
ul = Re[box[[1]]] + I*Im[box[[2]]]},
If[Round[(argchange[ll, lr, ZetaPrime[ll], ZetaPrime[lr]] +
argchange[lr, ur, ZetaPrime[lr], ZetaPrime[ur]] +
argchange[ur, ul, ZetaPrime[ur], ZetaPrime[ul]] +
argchange[ul, ll, ZetaPrime[ul], ZetaPrime[ll]])/(2 Pi)] > 0,
Check[
FindRoot[ZetaPrimeFine[s], {s, Mean[box], Sequence @@ box},
WorkingPrecision -> 30, AccuracyGoal -> 5] ,
Print[{Re[box[[1]]], Re[box[[2]]]},
" ", {Im[box[[1]]], Im[box[[2]]]}];
Flatten[
Table[{ii + I jj, ii + step/10 + I (jj + step/10)},
Evaluate@{ii,
Sequence @@ ({Re[box[[1]]], Re[box[[2]]]} - {0, step/10}),
step/10},
Evaluate@{jj,
Sequence @@ ({Im[box[[1]]], Im[box[[2]]]} - {0, step/10}),
step/10}], 1];
Map[SmartFindRoot, %] // Flatten, {FindRoot::reged}], {}]];
This function takes the lower left and upper right corners of the box, checks if there is a root in the box, and tries to find it with FindRoot. If a FindRoot::reged error occurs, Check calls SmartFindRoot with a smaller step size.
It works just fine when no FindRoot::reged error occurs. But when it does, the SmartRoot variable box_ does not seem to change in the recursion, and I go into an infinite loop:
step = 1/10;
re = {1/2, 3};
im = {45*10^5 + 30, 45*10^5 + 40};
Flatten[Table[{ii + I jj, ii + step + I (jj + step)},
Evaluate@{ii, Sequence @@ (re - {0, step}), step},
Evaluate@{jj, Sequence @@ (im - {0, step}), step}], 1];
Map[SmartFindRoot, %] // Flatten
results in
FindRoot::reged: The point {0.640762736317106198242023358404+4.50003870000000000000000000000\[CenterDot]10^6 I} is at the edge of the search region {0.60000000000000000000000+4.50003870000000000000000000000\[CenterDot]10^6 I,0.70000000000000000000000+4.50003880000000000000000000000\[CenterDot]10^6 I} in coordinate 1 and the computed search direction points outside the region.
{3/5,7/10} {45000387/10,22500194/5}
FindRoot::reged: The point {0.640762736317106198242023358404+4.50003870000000000000000000000\[CenterDot]10^6 I} is at the edge of the search region {0.60000000000000000000000+4.50003870000000000000000000000\[CenterDot]10^6 I,0.70000000000000000000000+4.50003880000000000000000000000\[CenterDot]10^6 I} in coordinate 1 and the computed search direction points outside the region.
{3/5,7/10} {45000387/10,22500194/5}
FindRoot::reged: The point {0.640762736317106198242023358404+4.50003870000000000000000000000\[CenterDot]10^6 I} is at the edge of the search region {0.60000000000000000000000+4.50003870000000000000000000000\[CenterDot]10^6 I,0.70000000000000000000000+4.50003880000000000000000000000\[CenterDot]10^6 I} in coordinate 1 and the computed search direction points outside the region.
General::stop: Further output of FindRoot::reged will be suppressed during this calculation.
{3/5,7/10} {45000387/10,22500194/5}
$Aborted
What am I doing wrong? I tried Block and With instead of Module but was unsuccessful.
Edit: In response to the comment, I can apply the SmartFindRoot function with a smaller step size to the the box causing the FindRoot::reged error, and it successfully finds the root:
step = 1/100;
re = {3/5, 7/10};
im = {45*10^5 + 38 + 7/10, 45*10^5 + 38 + 8/10};
Flatten[Table[{ii + I jj, ii + step + I (jj + step)},
Evaluate@{ii, Sequence @@ (re - {0, step}), step},
Evaluate@{jj, Sequence @@ (im - {0, step}), step}], 1];
Map[SmartFindRoot, %] // Flatten
{s -> 0.600152422567210509655436583395 +
4.50003870332967275482674049111\[CenterDot]10^6 I}
Second edit: I have a solution in my answer below; but I don't understand why it works or what the problem was. Happy to accept any answer that can enlighten me.
FindRootis converging towards a root just outside the box. From the documentation, I understand thatFindRootwill stop if condition 2 happens. Perhaps you need to restart from somewhere else in the box in this case. $\endgroup$