# 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. – Nasser May 20 '15 at 23:26
• @Nasser, OP said "nonzero", but did not mention "tiny"; it may well be that his equation actually has complex roots. – J. M. will be back soon May 20 '15 at 23:38
• @J.M. I did try it before commenting. They are tiny and they all went away after Chop. V 10.1 – Nasser May 20 '15 at 23:42
• They are actually supposed to be positive. That was the reason I explicitly took out the imaginary unit. – Alex May 20 '15 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. – J. M. will be back soon May 20 '15 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 Nov 15 '16 at 4:37
• 1*^-15 is Mathematica's shorthand "scientific notation." It means 1 * 10^-15. – Michael E2 Nov 15 '16 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 Nov 15 '16 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. – Michael E2 Nov 15 '16 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], {α, #[]}] & /@ seeds To verify that the imaginary parts do vanish, increase WorkingPrecision:

FindRoot[WhittakerW[1, (I α)/2, 10], {α, #[]},
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. – Michael E2 May 21 '15 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. – J. M. will be back soon May 21 '15 at 2:42
• @J.M. Do you mean like this? – Michael E2 May 21 '15 at 11:27
• Heh, I forgot that I wrote an answer to that effect already. Thanks! – J. M. will be back soon May 21 '15 at 11:35