Approximating a function with a Chebyshev interpolation is a powerful way to find roots. See links at the beginning of the code dump. Best to use an analytic function with this method (i.e., drop Abs
, which is not differentiable). Takes about 5 sec. to find cc
and rts50
.
cc = adaptiveChebSeries[
x |-> Hypergeometric1F1[1/4 (3 - (2*x)/I), 15/10, I], 1, 20000,
"Points" -> 2^8, WorkingPrecision -> 50];
Length@cc
(* 130 *)
rts50 = chebRoots[cc, {1, 20000}];
Length@rts50
(* 63 *)
ReImPlot[
Hypergeometric1F1[1/4 (3 - (2*x)/I), 15/10, I], {x, 1, 20000},
Mesh -> {rts50}, MeshStyle -> Red,
ScalingFunctions -> {"Log", Automatic}]
Code dump
Code taken from my answer here; see my answer here for a discussion of the code.
(* 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 < 3 && # <= 2^OptionValue[MaxIterations] OptionValue["Points"] &
(* at least two coefficients dropped *)
];
If[len < 3,
Message[adaptiveChebSeries::cvmit, OptionValue[MaxIterations]]];
If[TrueQ[len > 1], Drop[cc, 1 - len], cc]
]
(* "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}]
];
(* Find the real roots of a truncated Chebyshev series
representing a function over an interval [a,b] *)
Options[chebRoots] = {(* TBD *)};
chebRoots::usage =
"chebRoots[c,{a,b}], c = {c0, c1,..., cn} Chebyshev coefficients, over the interval {a,b}
computes the (real) roots in the interval {a,b}";
chebRoots[coeff_, dom_: {-1, 1}, OptionsPattern[]] := Module[{eigs},
eigs = Eigenvalues@colleagueMatrix[coeff];
roots = Sort@Rescale[
Re@Select[
eigs,
Abs[Im[#]] < 1*^-15 && (* Im error tolerance *)
-1.0001 < Re[#] < 1.0001 &], (* Re error tolerance *)
{-1, 1}, dom]
]
Quiet@FindRoot[Abs[Hypergeometric1F1[1/4 (3 - (2*x)/I), 1.5, I]] == 0, {x, #, 0, a}] & /@ Range[1, a, 20] // Values // Flatten // Union[#, SameTest -> (Abs[#1 - #2]/Abs[#1] < 10^(-3) &)] &
. $\endgroup$