2 replaced http://mathematica.stackexchange.com/ with https://mathematica.stackexchange.com/
source | link

This is a good example of the "dynamic range" problem in Boyd's CPR method (see also this answer by J.M.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.

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.

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.

1
source | link

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"}]

Mathematica graphics

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]

Mathematica graphics

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
* *)
Clear[adaptiveChebSeries];
Options[adaptiveChebSeries] = {PrecisionGoal -> 14, "Points" -> 32, 
   WorkingPrecision -> MachinePrecision, MaxIterations -> 5};
adaptiveChebSeries[f_, a_, b_, OptionsPattern[]] := 
 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, 
   Message[adaptiveChebSeries::cvmit, OptionValue[MaxIterations]]];
  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}]
   ];