Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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:

enter image description here

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

share|improve this question
1  
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

2 Answers 2

up vote 5 down vote accepted

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.)

share|improve this answer
    
@ 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 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. It comes with a Seeds option to control how many random seeds are used (50 by default), as well as 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.

Options[FindComplexRoots] = Join[Options[FindRoot], 
  {Seeds -> 50, SeedGenerator -> RandomComplex, Tolerance -> Automatic, Verbose -> False}];
Quiet[Seeds = Seeds, {Set::wrsym}]; Protect[Seeds];
Quiet[SeedGenerator = SeedGenerator, {Set::wrsym}]; Protect[SeedGenerator];
FindComplexRoots::usage = "FindComplexRoots[e1==e2, {z, zmin, zmax}] attempts to find 
  complex roots of the equation e1==e2 in the rectangle with corners zmin and zmax.";
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,∞).";
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, or any function f
  such that f[{zmin, zmax}, n] returns n complex numbers in the rectancle with corners 
  zmin and zmax."
$MessageGroups = Join[$MessageGroups, {"FindComplexRoots" :> {FindRoot::lstol}}]

The function itself is

FindComplexRoots[e1_ == e2_, {z_, zmin_, zmax_}, ops : OptionsPattern[]] := Module[
  {seeds, msglist},
  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][{zmin, zmax}, OptionValue[Seeds]];

  If[OptionValue[Verbose], Hold[], Hold[FindRoot::lstol]] /. {
   Hold[messageSequence___] :> Quiet[
   DeleteDuplicates[
    Select[
     Check[
        FindRoot[e1 == e2, {z, #},
         Evaluate[Sequence @@ FilterRules[{ops}, Options[FindRoot]]]
         ]
        , ## &[]
        ] & /@ seeds
     , (Re[zmin] < (Re[z] /. #) < Re[zmax] && 
        Im[zmin] < (Im[z] /. #) < Im[zmax]) &]
    , Abs[(z /. #1) - (z /. #2)] < If[
       NumberQ[OptionValue[Tolerance]],
       OptionValue[Tolerance],
       10^If[NumberQ[OptionValue[WorkingPrecision]], 
           2 - OptionValue[WorkingPrecision], 2 - $MachinePrecision]
       ] &] 
  , {messageSequence}]}
  ]

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, but if you try enough times then this should not be a problem. If you don't get any roots out, 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 graphics

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[{zmin_, zmax_}, number_] := (
  ((Re[zmax] - Re[zmin]) #1 + Re[zmin] + 
      I ((Im[zmax] - Im[zmin]) #2 + Im[zmin])) & @@@ BlockRandom[
    SeedRandom[Method -> {"MKL", Method -> {"Sobol", "Dimension" -> 2}}];
    SeedRandom[];
    RandomReal[{0, 1}, {number, 2}]
    ]
  )

RandomNiederreiterComplexes[{zmin_, zmax_}, number_] := (
  ((Re[zmax] - Re[zmin]) #1 + Re[zmin] + 
      I ((Im[zmax] - Im[zmin]) #2 + Im[zmin])) & @@@ BlockRandom[
    SeedRandom[ Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> 2}}];
    SeedRandom[];
    RandomReal[{0, 1}, {number, 2}]
    ]
  )

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

DeterministicComplexGrid[{zmin_, zmax_}, number_] := Module[
  {p = Round[ Sqrt[(Re[zmax] - Re[zmin])/(Im[zmax] - Im[zmin]) number]], 
   q = Round[ Sqrt[(Im[zmax] - Im[zmin])/(Re[zmax] - Re[zmin]) number]]},
  Flatten[Outer[#1 + I #2 &,
    Most@Rest@N[Range[Re[zmin], Re[zmax], (Re[zmax] - Re[zmin])/(p + 2)]],
    Most@Rest@N[Range[Im[zmin], Im[zmax], (Im[zmax] - Im[zmin])/(q + 2)]]
    ], 2]
  ]

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.)

Mathematica graphics

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.