Using NSolve in the complex plane
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
Verboseoption which controls whetherFindRoot::lstolmessages are displayed, and which works as described in this answer. In general you want to keepVerboseset toFalse, 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 toTruecan help you figure out what's going on.The number of
Seedsshould 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 theWorkingPrecision. These duplicates are removed, with the option to set theTolerancefor 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
seedsand then mapsFindRoot[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];
gridPointBasis = MapThread[
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.)
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.)
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.