# Roots of Whittaker W function

I am interested in finding the roots $u$ of the equation $$W_{1,\imath b}(a)=0,$$ where $W_{\kappa,\mu}(z)$ denotes the Whittaker $W$ function, $a>0$ is a fixed parameter, $\imath=\sqrt{-1}$ and $b$ is the unknown. Due symmetry of the Whittaker $W$ with respect to $\mu$, it suffices to consider the above equation only for $b>0$. For any given $a>0$ the equation has countably many solutions, so numerically there is this issue of somehow telling Mathematica to ignore the solution it has already found.

I've written the following script

ClearAll[Evaluate[Context[] <> "*"]];

ϕ[x_, α_] := N[WhittakerW[1, (I α)/2, x]];
Subscript[α, opt][a_, β_] := α /.
FindRoot[ϕ[a, α], {α, β}];

a = 10;
β = 0;
For[i = 1, i <= 10, i += 1,
β =
Subscript[α, opt][
a, β*1.1]; (*tweak the starting point a bit so as to force \
Mathematica to look for zeros to the right of the previous one*)
Print["a=", a,
", \!$$\*SubscriptBox[\(α$$, $$opt$$]\)=", β,
", ϕ[a,\!$$\*SubscriptBox[\(α$$, $$opt$$]\)]=", ϕ[
a, β]];
]


I doubt though that the solutions the script outputs are correct. One issue is that they all have a nonzero imaginary component. Any suggestions as to how to improve the script?

• One issue is that they all have a nonzero imaginary component To clean, use Chop it is floating point computation noise. Commented May 20, 2015 at 23:26
• @Nasser, OP said "nonzero", but did not mention "tiny"; it may well be that his equation actually has complex roots. Commented May 20, 2015 at 23:38
• @J.M. I did try it before commenting. They are tiny and they all went away after Chop. V 10.1 Commented May 20, 2015 at 23:42
• They are actually supposed to be positive. That was the reason I explicitly took out the imaginary unit.
– Alex
Commented May 20, 2015 at 23:45
• I have some suggestions: convert your Whittaker function to a confluent hypergeometric function with FunctionExpand[], remove the exponential and power factors multiplying HypergeometricU[], and try working with that function instead. Commented May 20, 2015 at 23:46

This is a good example of the "dynamic range" problem in Boyd's CPR method (see also this answer by J.M.). As α moves toward -100, the OP's oscillatory Whittaker function decays exponentially. According to Boyd, where $ε_{mach}$ is the machine epsilon and $f_{max}$ is the maximum of $f(x)$ over an interval $[a,b]$, "If a function $f(x)$ has a magnitude as small as $|f (x)| \le O(ε_{mach}\, f_{max})$ on a subinterval within $[a,b]$ (excluding tiny subintervals around a zero), then $f(x)$ has, to borrow a term from photography, a large dynamic range." This can lead to numerical problems since the Chebyshev polynomials oscillate between ±1 throughout the interval. A Chebyshev series that is computed with floating-point numbers may have trouble representing a range greater than $ε_{mach}\, f_{max}$ to $f_{max}$. This is negligible if $f(x)$ passes quickly through the range $\pm ε_{mach}\, f_{max}$, such as near a zero of $f(x)$. There are two approaches to dealing with this issue. One can scale the function to alleviate the dynamic range, or one can subdivide the interval and find the roots piece-wise. We will show the scaling approach. The piecewise approach is feasible (with a interval length of 10 and WorkingPrecision -> 40), but we will show the scaling approach, which seems particularly appropriate here.

From the plot below, we can see that the rate of exponential decay is approximately

k = Rationalize[Log[10^0.3372], 0.01]
(*  7/9  *)

plot = Plot[RealExponent@WhittakerW[1, (I a)/2, 10], {a, -100, 0}, PlotPoints -> 100];
Manipulate[
Show[
plot,
Graphics[{Dynamic[Line[{{0, b}, {-100, -100 m + b}}]]}]],
{m, 0, 1, Appearance -> "Labeled"}, {b, 0, 2, Appearance -> "Labeled"}]


Alternatively, we can estimate the decay rate numerically from two distant peaks:

{m1, a1} = FindMaximum[Re@WhittakerW[1, (I a)/2, 10], {a, -40., -39.9}];
{m2, a2} = FindMaximum[Re@WhittakerW[1, (I a)/2, 10], {a, -100., -99.9}];
k = Rationalize[Log[m1/m2]/Subtract @@ (a /. {a1, a2}), 0.01]
(*  7/9  *)


We can scale the Whittaker function to have a fairly even variation over the interval {-100, 0}.

Plot[Exp[-k α] WhittakerW[1, (I α)/2, 10], {α, -100, 0},
WorkingPrecision -> 20, PlotRange -> All]


The function adaptiveChebSeries computes a Chebyshev series of a function $f(x)$ over $a \le x \le b$ to a given precision goal. The series $\sum c_j \, T_j(w)$ is in terms of $w = (2x - (a+b))/(b-a)$, where $-1 \le w \le 1$ for $a \le x \le b$. One can then compute companion matrix of the series, sometimes called the "colleague matrix," whose eigenvalues $w_j$ are the roots of the Chebyshev series. The roots $x_j$ of $f(x)$ will be the rescaled eigenvalues $x_j = ((b-a)w_j/2+(a+b)/2)$. The eigenvalues often have

(a0 = -100; b0 = 0;
coeff = adaptiveChebSeries[                 (* get Chebyshev coefficients *)
Exp[-k #] WhittakerW[1, (I #)/2, 10] &,
a0, b0, PrecisionGoal -> 14.7];
eigs = Eigenvalues@colleagueMatrix[coeff];  (* eigenvals of matrix contain the roots *)
seeds = Sort@Rescale[                       (* select roots in [-1, 1] and *)
Re@Select[                               (* rescale to [a0, b0] *)
eigs,
Abs[Im[#]] < 1*^-15 && -1.0001 < Re[#] < 1.0001 &],
{-1, 1}, {a0, b0}];) // AbsoluteTiming
Length[seeds]


Eigenvalues::arh: Because finding 126 out of the 126 eigenvalues and/or eigenvectors is likely to be faster with dense matrix methods, the sparse input matrix will be converted. If fewer eigenvalues and/or eigenvectors would be sufficient, consider restricting this number using the second argument to Eigenvalues. >>

(*
{0.87655, Null}
32
*)


The last coefficient gives an estimate of the error bound.

Abs@Last@coeff
(*  6.80667*10^-15  *)


One can polish the eigenvalue-roots with FindRoot.

roots = α /. FindRoot[WhittakerW[1, (I α)/2, 10], {α, #}] & /@ seeds; // AbsoluteTiming
(*  {1.44508, Null}  *)


The seeds are pretty close to the roots. The residuals of the seeds and the roots are appropriately small.

(seeds - roots)/roots // Abs // Max
WhittakerW[1, (I roots)/2, 10] // Abs // Max
WhittakerW[1, (I seeds)/2, 10] // Abs // Max
(*
3.54593*10^-15
9.20479*10^-21
1.09484*10^-18
*)


Since the Whittaker function has a fairly regular exponential decay, scaling by Exp[-k x] will give an approximate idea of the relative accuracy of the roots.

Exp[-k roots] WhittakerW[1, (I roots)/2, 10] // Abs // Max
Exp[-k seeds] WhittakerW[1, (I seeds)/2, 10] // Abs // Max
(*
1.47868*10^-13
6.84676*10^-12
*)


Code dump

The error in a Chebyshev series $\sum_{j=0}^N c_n T_n(x)$ for an analytic function is bounded by $\sum_{n > N} |c_n|$ and approximately equal to $|c_{N+1}| \le |c_{N}|$ if the series converges geometrically (or exponentially) and rapidly, since $|T_n(x)| \le 1$ for $-1 \le x \le 1$. If the function is analytic in a complex neighborhood of the interval $[-1,1]$, then the series will eventually converge geometrically. The coefficients are not guaranteed to decrease monotonically. For instance the Chebyshev polynomials are either even or odd functions, so for an even or odd function $f(x)$, every other Chebyshev coefficient will be zero. To help protect against mistakenly underestimating the error in such cases, we estimate the error by summing at least two terms of the tail of a series. Of course it is still possible for a function to be orthogonal to any number of successive Chebyshev polynomials, so such an error estimate can sometimes be fooled.

(* Chebyshev extreme points *)
chx[n_, prec_: MachinePrecision] := N[Cos[Range[0, n]/n Pi], prec];

(* Chebyshev series approximation to f *)
Clear[chebSeries];
chebSeries[f_, a_, b_, n_, prec_: MachinePrecision] := Module[{x, y, cc},
x = Rescale[chx[n, prec], {-1, 1}, {a, b}];
y = f /@ x; (* function values at Chebyshev points *)

cc = Sqrt[2/n] FourierDCT[y, 1]; (* get coeffs from values *)

cc[[{1, -1}]] /= 2;  (* adjust first & last coeffs *)
cc
];

(* recursively double the Chebyshev points
* until the PrecisionGoal is met
* The function values are memoized in f0
* *)
Options[adaptiveChebSeries] = {PrecisionGoal -> 14, "Points" -> 32,
WorkingPrecision -> MachinePrecision, MaxIterations -> 5};
Module[{cc, f0, max, len = 0, sum},
f0[x_] := f0[x] = f[x]; (* values reused in subsequent series *)
NestWhile[
(cc = chebSeries[f0, a, b, #, OptionValue[WorkingPrecision]];
(* check error estimate *)
max = Max@Abs@cc;   (*
sum the tail of the series *)
sum = 0;            (* relative to the max coefficient *)
len = LengthWhile[
Reverse@cc, (sum += Abs@#) < 10^-OptionValue[PrecisionGoal]*max &];
2 #) &,  (* double the number of points *)
OptionValue["Points"],
len < 2 && # <= 2^OptionValue[MaxIterations] OptionValue["Points"] &     (*
at least two coefficients dropped *)
];
If[len < 2,
Drop[cc, -len]
]

(* "Chebyshev companion matrix" (Boyd, 2014) /
"Colleague matrix" (Good, 1961) *)
colleagueMatrix[cc_] := With[{n = Length[cc] - 1},
SparseArray[{{i_, j_} /; i == j + 1 :> 1/2, {i_, j_} /; i + 1 == j :>
1/(2 - Boole[i == 1])}, {n, n}] -
SparseArray[{{n, i_} :> cc[[i]]/(2 cc[[n + 1]])}, {n, n}]
];

• One quick question: what's 1*^-15 in Abs[Im[#]] < 1*^-15? I don't understand the syntax.
– Alex
Commented Nov 15, 2016 at 4:37
• 1*^-15 is Mathematica's shorthand "scientific notation." It means 1 * 10^-15. Commented Nov 15, 2016 at 11:52
• I see, thank you. One more question: I would like the precision of FindRoot to be say 100 decimal places. I tried setting WorkingPrecision to 200, which should automatically set PrecisionGoal and AccuracyGoal to 100 each. However, at the zeros returned by FindRoot the target function does not evaluate to a number as small as $10^{-100}$. Is it fixable?
– Alex
Commented Nov 15, 2016 at 15:42
• @Alex I'm not sure whether you're asking in general, but, in general, that can happen if the derivative is very large, or if there are significant rounding errors in computing the function, or maybe some other things. I usually get an warning message. The first thing I try is setting the AccuracyGoal explicitly and increasing WorkingPrecision to 3 to 10 times the accuracy goal. If raising WP doesn't get me closer, then I wonder what's special about the problem. Commented Nov 15, 2016 at 16:07

You can use a plot to find the roots and then follow up with FindRoot to get precise solutions. Here I'll use ContourPlot to plot where the real part vanishes and use MeshFunctions to find where the imaginary part is simultaneously zero.

cp = Show[
Table[
ContourPlot[
Evaluate[Re[WhittakerW[1, (I α)/2, 10] /. α -> a + b I] == 0],
{a, a1, a1 + 10}, {b, -0.1, 0.1},
PlotPoints -> {21, 5}, Mesh -> {{0}},
MeshFunctions -> Function[{a, b}, Im[WhittakerW[1, (I α)/2, 10] /. α -> a + b I]],
MeshStyle -> Red],
{a1, -100, -10, 10}],
PlotRange -> All]


seeds = Sort@Cases[Normal@cp, Point[p_] :> p, Infinity];
seeds // Short[#, 3] &
(*  {{-98.9014,-0.0024671},<<30>>,{-14.6835,-1.11492*10^-7}}  *)

FindRoot[WhittakerW[1, (I α)/2, 10], {α, #[[1]]}] & /@ seeds


To verify that the imaginary parts do vanish, increase WorkingPrecision:

FindRoot[WhittakerW[1, (I α)/2, 10], {α, #[[1]]},
WorkingPrecision -> 30] & /@ seeds

• Using a plot to find roots has certainly been done before. I feel like there is a similar question on the site dealing with the roots of a complex-valued function. I couldn't find one, though. Commented May 21, 2015 at 2:09
• I've used FindAllCrossings2D[] plenty of times in the past for seeking complex roots; I just feed it the real and imaginary parts separately. Commented May 21, 2015 at 2:42
• @J.M. Do you mean like this? Commented May 21, 2015 at 11:27
• Heh, I forgot that I wrote an answer to that effect already. Thanks! Commented May 21, 2015 at 11:35