# Using NSolve in the complex plane

I am trying to find numerical values for $\Im\ e^{\frac{1}{2} +i\ y} = \Im\ \zeta(\frac{1}{2} +i\ y)$ for a given range. I have tried:

NSolve[Im[Exp[0.5 + I x]] == Im[Zeta[0.5 + I x]] && 0 <= x <= 20, x]


I can see that:

Plot[{Im[Exp[0.5 + I x]], Im[Zeta[0.5 + I x]]}, {x, 0, 20}]


gives:

... I am sure it is a problem with my syntax again, but can't seem to fix it :/

• Why don't you use FindRoot? We discussed the same issue here. Now you might try this: FindRoot[ Im[Zeta[1/2 + I t] - Exp[1/2 + I t]], {t, #, 0, 20}] & /@ {3, 5}. It's hard to say definitely what's wrong with NSolve here, but similar problems one encounters here Issue with NSolve. – Artes Nov 12 '13 at 18:07
• @Artes, I was just looking at that post, trying to figure it out. Is there any way (rather than manually going back & plugging in values) to output full coordinates? Also, is it possible to output all values? FindRoot only outputs two values for the given range. – martin Nov 12 '13 at 18:14
• ... what does the {3, 5} stand for? – martin Nov 12 '13 at 18:17
• It appears that entering {3, 6, 5} gives the first three values, but I have no idea why. – martin Nov 12 '13 at 18:31
• Another problem: entering t = 3.1701346526064005; NSolve[Im[Exp[1/2 + I t]]] gives error message NSolve::precsm: Requested precision -0.0470514 is smaller than MinPrecision. Using MinPrecision instead.  :/ – martin Nov 12 '13 at 18:37

This solution will only work for zeros of odd multiplicity that are spaced far enough apart. We can just look at the sign changes from Plot.

Options[FindAllRootsInRange] = Options[FindRoot];
FindAllRootsInRange[e1_ == e2_, {x_, a_, b_}, ops:OptionsPattern[]] := Module[{p, s, g},
p = Plot[e1 - e2, {x, a, b}];
s = SplitBy[
p[[1, 1, 3, 2, 1]],
Sign[Last[#]] &
];
g = s[[All, 1, 1]];
FindRoot[e1 == e2, {x, #}, ops] & /@ g
]

FindAllRootsInRange[Im[Exp[0.5 + I x]] == Im[Zeta[0.5 + I x]], {x, 0, 20},
WorkingPrecision -> 20]

(* {{x -> -4.8148248609680896326*10^-35},
{x -> 3.1701346526063997557},
{x -> 6.5144026505762818508},
{x -> 9.3682393343876903168},
{x -> 12.091009965923843486},
{x -> 15.176592892192292174},
{x -> 18.423560419608622067}} *)


To be more thorough, one could catch all zeros of even multiplicity by applying this function to the derivative too, then choosing which ones are roots to the original function. (That's not necessary for this example though.)

• @ RiemannZeta, thanks for your input. My nachine is outputting x -> 1.2037062152420224082*10^-33 for the first value. What is this discrepancy in accuracy about? – martin Nov 12 '13 at 21:39
• @martin perhaps we're using different version? One could fix this by putting Chop around FindRoot. – Chip Hurst Nov 12 '13 at 22:23

The current answer by Chip Hurst deals well with the case of solving for the roots of a (possibly) complex-valued function of a real variable, which NSolve struggles with. However, the more general case of finding roots of a complex-valued equation of a complex variable is a good bit more complicated, and it requires a stronger approach.

One thing to note is that NSolve is often fine with complex-valued equations of complex variables (such as NSolve[Sin[t] == 2, t]) but the class of functions it can solve for is limited, and does not include many interesting trascendental equations, like NSolve[Sin[t] == 2 + 0.01 t, t]. On the other hand, FindRoot works well in those cases (so try e.g. FindRoot[Sin[t] == 2 + 0.01 t, {t, π/2 - I}]) but it does not handle domain constraints, it only returns a single solution, and it requires a seed to get going.

The obvious way forward is therefore to start with a bunch of random seeds in the desired patch of the complex plane, and then to feed those seeds to FindRoot. If one tries with enough seeds, one can be reasonably certain that the procedure will find all the findable roots.

The function FindComplexRoots described below does exactly that. The function and some more documentation can be found in this GitHub repository.

The calling syntax is of the form

FindComplexRoots[e1==e2, {z, zmin, zmax}]


and it will attempt to find roots of the equation in the complex rectangle whose corners are zmin and zmax; its output is identical to that of NSolve and FindRoot. To do multiple equations, you then do

FindComplexRoots[{e1==e2, e3==e4, …}, {z1, z1min, z1max}, {z2, z2min, z2max}, …]


It comes with a Seeds option to control how many random seeds are used (50 by default), a SeedGenerator option to control how the seeds are created, and a few other niceties.

Let me start with a bunch of options and documentation. This includes syntax information so you'll get pretty colouring and a usage message so autocomplete templates will work in versions that have them.

FindComplexRoots::usage="FindComplexRoots[e1==e2, {z, zmin, zmax}] attempts to find complex roots of the equation e1==e2 in the complex rectangle with corners zmin and zmax.

FindComplexRoots[{e1==e2, e3==e4, …}, {z1, z1min, z1max}, {z2, z2min, z2max}, …] attempts to find complex roots of the given system of equations in the multidimensional complex rectangle with corners z1min, z1max, z2min, z2max, ….";
Seeds::usage="Seeds is an option for FindComplexRoots which determines how many initial seeds are used to attempt to find roots of the given equation.";
SeedGenerator::usage="SeedGenerator is an option for FindComplexRoots which determines the function   used to generate the seeds for the internal FindRoot call. Its value can be RandomComplex, RandomNiederreiterComplexes, RandomSobolComplexes, DeterministicComplexGrid, or any function f such that f[{zmin, zmax}, n] returns n complex numbers in the rectancle with corners zmin and zmax.";

Options[FindComplexRoots] = Join[Options[FindRoot], {Seeds -> 50, SeedGenerator -> RandomComplex, Tolerance -> Automatic, Verbose -> False}];
SyntaxInformation[FindComplexRoots] = {"ArgumentsPattern" -> {_, {_, _, _}, OptionsPattern[]},    "LocalVariables" -> {"Table", {2, ∞}}};
FindComplexRoots::seeds = "Value of option Seeds -> 1 is not a positive integer.";
FindComplexRoots::tol = "Value of option Tolerance -> 1 is not Automatic or a number in [0,∞).";
$MessageGroups=Join[$MessageGroups,{"FindComplexRoots":>{FindRoot::lstol}}]

Protect[Seeds];
Protect[SeedGenerator];


The function itself is

FindComplexRoots[equations_List, domainSpecifiers__, ops : OptionsPattern[]] :=
Block[{seeds,tolerances},

If[
! IntegerQ[Rationalize[OptionValue[Seeds]]] || OptionValue[Seeds]<=0,
Message[FindComplexRoots::seeds, OptionValue[Seeds]]
];
If[
! (OptionValue[Tolerance] === Automatic || OptionValue[Tolerance]>=0),
Message[FindComplexRoots::tol, OptionValue[Seeds]]
];

seeds=OptionValue[SeedGenerator][
{domainSpecifiers}[[All,{2,3}]],OptionValue[Seeds]
];
tolerances=Which[
ListQ[OptionValue[Tolerance]],OptionValue[Tolerance],
True,ConstantArray[
Which[
NumberQ[OptionValue[Tolerance]],OptionValue[Tolerance],
True,10^If[NumberQ[OptionValue[WorkingPrecision]], 2-OptionValue[WorkingPrecision],2-\$MachinePrecision]
]
,Length[{domainSpecifiers}]]
];

If[OptionValue[Verbose],Hold[], Hold[FindRoot::lstol]] /. {
Hold[messageSequence___] :> Quiet[
DeleteDuplicates[
Select[
Check[
FindRoot[
equations
,Evaluate[Sequence@@Table[{{domainSpecifiers}[[j,1]],#[[j]]},{j,Length[{domainSpecifiers}]}]]
,Evaluate[Sequence @@ FilterRules[{ops}, Options[FindRoot]]]
],
##&[]
]&/@seeds,
Function[
repList,
ReplaceAll[
Evaluate[And@@Table[
And[
Re[{domainSpecifiers}[[j,2]]]<=Re[{domainSpecifiers}[[j,1]]]<=Re[{domainSpecifiers}[[j,3]]],
Im[{domainSpecifiers}[[j,2]]]<=Im[{domainSpecifiers}[[j,1]]]<=Im[{domainSpecifiers}[[j,3]]]
]
,{j,Length[{domainSpecifiers}]}]]
,repList]
]
],
Function[{repList1,repList2},
And@@Table[
Abs[({domainSpecifiers}[[j,1]]/.repList1)-({domainSpecifiers}[[j,1]]/.repList2)]<tolerances[[j]]
,{j,Length[{domainSpecifiers}]}]
]
]
, {messageSequence}]}
]
FindComplexRoots[e1_==e2_,{z_,zmin_,zmax_},ops:OptionsPattern[]]:=FindComplexRoots[{e1==e2},{z,zmin,zmax},ops]


A few notes on the implementation:

• There is a Verbose option which controls whether FindRoot::lstol messages are displayed, and which works as described in this answer. In general you want to keep Verbose set to False, because if you try many times then you are likely to hit some ill-conditioned seeds or generally some bad stuff, but if you try enough times then this should not be a problem. If you don't get any roots out, or you suspect the output of being fishy, then setting it to True can help you figure out what's going on.

• The number of Seeds should be set high enough that you're reasonably certain to get all the roots you need. When that happens, though, you will get many duplicates, which will be off from each other by one or two digits less than the WorkingPrecision. These duplicates are removed, with the option to set the Tolerance for how different they're allowed to be.

• Other than a bunch of machinery to make it nicer, the function is quite simple. In essence, it sets a bunch of random seeds and then maps FindRoot[e1 == e2, {z, #}]& over that list.

One thing that has a large impact on how well this works is the choice of initial seeds. Random seeds tend to work well, but they do require a relatively large number of iterations to work well consistently. Another alternative would be to use an evenly distributed grid, but that does mean that if you hit a bad grid then you're going to have a hard time getting off it. This can be fixed by using a randomly-perturbed uniform grid, which you can implement if you want.

A nice middle ground between the two, on the other hand, is to use a low-discrepancy quasirandom sequence. These are nice because they tend to cluster a lot less and therefore require fewer seeds to cover a given square patch of complex plane uniformly.

Mathematica comes with two types of low-discrepancy sequences: Niederreiter and Sobol sequences. Unfortunately these are only implemented for Intel chips, and they are described under the "MKL" method in the Random Number Generation tutorial. If you are not running an Intel chip the functions below won't work (presumably), so you will need to roll your own as per the specs in the SeedGenerator::usage message, or try to obtain the QRSToolbox package described in this paper.

RandomSobolComplexes::usage="RandomSobolComplexes[{zmin, zmax}, n] generates a low-discrepancy Sobol sequence of n quasirandom complex numbers in the rectangle with corners zmin and zmax.

RandomSobolComplexes[{{z1min,z1max},{z2min,z2max},…},n] generates a low-discrepancy Sobol sequence of n quasirandom complex numbers in the multi-dimensional rectangle with corners {z1min,z1max},{z2min,z2max},….";

RandomSobolComplexes[pairsList__, number_] :=Map[
Function[randomsList,
pairsList[[All,1]]+Complex@@@Times[
ReIm[pairsList[[All,2]]-pairsList[[All,1]]],
randomsList
]
],
BlockRandom[
SeedRandom[Method->{"MKL",Method->{"Sobol", "Dimension" -> 2Length[pairsList]}}];
SeedRandom[];
RandomReal[{0, 1}, {number,Length[pairsList],2}]
]
]
RandomSobolComplexes[{zmin_?NumericQ,zmax_?NumericQ},number_]:=RandomSobolComplexes[{{zmin,zmax}},number][[All,1]]


RandomNiederreiterComplexes::usage="RandomNiederreiterComplexes[{zmin, zmax}, n] generates a low-discrepancy Niederreiter sequence of n quasirandom complex numbers in the rectangle with corners zmin and zmax.

RandomNiederreiterComplexes[{{z1min,z1max},{z2min,z2max},…},n] generates a low-discrepancy Niederreiter sequence of n quasirandom complex numbers in the multi-dimensional rectangle with corners {z1min,z1max},{z2min,z2max},….";

RandomNiederreiterComplexes[pairsList__, number_] :=Map[
Function[randomsList,
pairsList[[All,1]]+Complex@@@Times[
ReIm[pairsList[[All,2]]-pairsList[[All,1]]],
randomsList
]
],
BlockRandom[
SeedRandom[Method->{"MKL",Method->{"Niederreiter", "Dimension" -> 2Length[pairsList]}}];
SeedRandom[];
RandomReal[{0, 1}, {number,Length[pairsList],2}]
]
]
RandomNiederreiterComplexes[{zmin_?NumericQ,zmax_?NumericQ},number_]:=RandomNiederreiterComplexes[{{zmin,zmax}},number][[All,1]]


To keep things on a more even footing, and to allow for a deterministic option, here is a corresponding deterministic uniform grid with the same syntax:

DeterministicComplexGrid::usage="DeterministicComplexGrid[{zmin, zmax}, n] generates a grid of about n equally spaced complex numbers in the rectangle with corners zmin and zmax.

DeterministicComplexGrid[{{z1min,z1max},{z2min,z2max},…},n] generates a regular grid of about n equally spaced complex numbers in the multi-dimensional rectangle with corners {z1min,z1max},{z2min,z2max},….";

DeterministicComplexGrid[pairsList_, number_] :=
Block[{sep, separationsList, gridPointBasis, k},
sep = NestWhile[0.99 # &,
Min[Flatten[ReIm[pairsList[[All, 2]] - pairsList[[All, 1]]]]],
Times @@ (
Floor[Flatten[ReIm[pairsList[[All, 2]] - pairsList[[All, 1]]]],
0.99 #]/(0.99 #)) <= number &];
separationsList =
Round[Floor[
Flatten[ReIm[pairsList[[All, 2]] - pairsList[[All, 1]]]], sep]/
sep];
Function[{l, n},
Range[l[[1]], l[[2]], (l[[2]] - l[[1]])/(n + 1)][[2 ;; -2]]],
{Flatten[Transpose[ReIm[pairsList], {1, 3, 2}], 1],
separationsList}
];
Flatten[Table[
Table[k[2 j - 1] + I k[2 j], {j, 1, Length[pairsList]}],
Evaluate[
Sequence @@
Table[{k[j], gridPointBasis[[j]]}, {j, 1,
2 Length[pairsList]}]]
], Evaluate[Range[1, 2 Length[pairsList]]]]
]
DeterministicComplexGrid[{zmin_?NumericQ, zmax_?NumericQ}, number_] :=
DeterministicComplexGrid[{{zmin, zmax}}, number][[All, 1]]


Note that the number of points returned by DeterministicComplexGrid is always ≤number but it need not match it. Because of the uniformity requirement, this scales pretty badly with the number of points, but of course you can roll your own to suit your needs and feed the pure function to the SeedGenerator option of FindComplexRoots.

Finally, the built-in random number generator needs some slight tweaking to allow for the same syntax in the multi-variable case:

Unprotect[RandomComplex];
RandomComplex[{range1_List, moreRanges___}, number_] :=
Transpose[RandomComplex[#, number] & /@ {range1, moreRanges}]
Protect[RandomComplex];


When you're faced with a single equation, the best procedure is to choose some pretty high value for Seeds and run it a few times to make sure you're good to go. If you're looping over many similar equations, on the other hand, you do want to choose a value that's big enough to ensure you find all the roots of all the equations, but not a bigger one so that it doesn't clog up performance.

To help with choosing that value, here is a benchmarking suite for this function. You first calculate for a bunch of different seed numbers and generators,

Table[
{seedGenerator, AbsoluteTiming[
benchmark[seedGenerator] = Flatten[Table[
{seedNumber, Length[#[[2]]], #[[1]]} &[AbsoluteTiming[
FindComplexRoots[
(1 + (1 - Sin[t])^2) (1 + (1 + Sin[t])^2) == 0.01 t, {t, -2 I, 2 \[Pi] + 2 I},
Seeds -> seedNumber, SeedGenerator -> seedGenerator
]
]]
, {seedNumber, 1, 30}, {repetition, 100}
], 1];
][[1]]}
, {seedGenerator, {RandomComplex, RandomSobolComplexes,
RandomNiederreiterComplexes, DeterministicComplexGrid}}
]


and then you graph how many roots were found for each number of seeds, with the dot size representing how many iterations found that number.

GraphicsColumn[
Table[
Show[
Graphics[
{PointSize[0.0002 #2], Point[#1]}] & @@@ Tally[
benchmark[seedGenerator][[All, {1, 2}]]
]
, Frame -> True
, ImageSize -> 600
, AspectRatio -> 1/3
, PlotRange -> {{.5, 30.5}, {-0.5, 9}}
, PlotLabel -> ToString[seedGenerator]
]
, {seedGenerator, {RandomComplex, RandomSobolComplexes,
RandomNiederreiterComplexes, DeterministicComplexGrid}}]
]


The results are about what you'd expect. You need about three to five times more random seeds than roots to ensure you get them all with a ~99% certainty, and the quasirandom sequences slightly outperform the random ones. When the deterministic grid works, it tends to do so well, but again if you hit a bad spot then you can't get off it. (Note also that DeterministicComplexGrid[{zmin, zmax}, n] need not return a list of length n.)

• Just seeing this, very nice! +1 – Chip Hurst Aug 8 '15 at 16:53

It has been noted a fair number of times on this site that one can use the MeshFunctions option of Plot[] to help in root-finding. Applied to this case:

ie[x_] := Im[Exp[1/2 + I x]]
iz[x_] := Im[Zeta[1/2 + I x]]

pic = Plot[{ie[x], iz[x]}, {x, 0, 20}, Mesh -> {{0.}},
MeshFunctions -> {(ie[#] - iz[#]) &},
MeshStyle -> Directive[Red, AbsolutePointSize[6]]]


(* needed to remove duplicate meshes given to each function *)
seeds = DeleteDuplicates[Cases[Normal[pic], Point[pt_] :> pt, ∞],
EuclideanDistance[##] < 0.1 &][[All, 1]];

With[{prec = 20},
(\[FormalT] /. First[FindRoot[ie[\[FormalT]] - iz[\[FormalT]],
{\[FormalT], SetPrecision[#, prec]}, WorkingPrecision -> prec]]) &
/@ seeds]
{3.1701346526063997557, 6.5144026505762811360, 9.3682393343876916520,
12.091009965923843486, 15.176592892192292174, 18.423560419608621761}


although the intersection at the origin is missed here.

In fact, the function FindAllCrossings[] internally performs a similar strategy:

FindAllCrossings[ie[x] - iz[x], {x, 0, 20}, WorkingPrecision -> 20]


which should give the same results as above.

For finding the complex roots of a complex-valued function, as was done in episanty's answer, you can use FindAllCrossings2D[], as was done here.

• nice! neat solution :) – martin Jul 4 '15 at 8:57