Find all roots of a function with parabolic cylinder functions in a range of the variable

I want to find all roots of a function involving Parabolic Cylinder Functions. In what follows, I define 2 variables $\xi1$ and $\xi2$, which in turn depend on $\omega$. My function is then defined as f. I go on defining g and h (where I take specific values for my parameters $M$ and $\lambda$ which are real and positive. I then plot the real and imaginary part of h to locate the roots. I would like, however, to be able to find all roots of $Im[h]$ (The real part is essentially 0) in a range of $\omega$, say from 0 to 50.

ξ1[ ω_] := (-1 + I) (ω/Sqrt[λ] - Sqrt[λ]/2)

ξ2[ ω_] := (-1 + I) (ω/Sqrt[λ] + Sqrt[λ]/2)

f := (I ParabolicCylinderD[I M/(2 λ), I ξ1[ ω]] - Sqrt[M/λ]*(-I - 1)/2*
ParabolicCylinderD[I M/(2 λ) - 1, I ξ1[ ω]])*(
ParabolicCylinderD[-I M/(2 λ), ξ2[ ω]] + I *Sqrt[M/λ]*(I - 1)/2*
ParabolicCylinderD[-I M/(2 λ) - 1, ξ2[ ω]]) + (I ParabolicCylinderD[
I M/(2 λ), I ξ2[ ω]] + Sqrt[M/λ]*(-I - 1)/2*ParabolicCylinderD[I M/(2 λ) - 1,
I ξ2[ ω]])*(ParabolicCylinderD[-I M/(2 λ), ξ1[ ω]] - I *Sqrt[M/λ]*(I - 1)/2*
ParabolicCylinderD[-I M/(2 λ) - 1, ξ1[ ω]])
g:=FullSimplify[f, {ω>0&&M>0&&λ>0}]
h:=FullSimplify[g/.{M->2, λ->100}]
Plot[{ Re[h], Im[h]}, {ω, 0, 20}, PlotPoints -> 50, MaxRecursion -> 0]
FindRoot[Im[h],{ω,5}]


I have searched through some posts with the keyword "find all roots in a range"; however, most of the solutions are for simpler functions than this special parabolic cylinder functions, c.f. About multi-root search in Mathematica for transcendental equations and Find all roots of an interpolating function (solution to a differential equation).

I would appreciate any help. Thank you in advance.

• To help the reader, please provide values for all your constants. Also, do you mean f [ω_] rather than f? You could improve your question by providing the code to calculate at least one root. – bbgodfrey Jun 25 '15 at 18:30
• Thanks, I forgot some codes. I have updated. I mean f [ω_] as a function of ω, after assigning values for M and λ. Since I am not too familiar with Mathematica, so there is some anomaly in writing the codes. – Ambrose Chau Jun 25 '15 at 18:55

...most of the solutions are for simpler functions...

I'm not quite sure what gave OP that impression; certainly, FindAllCrossings[] is quite capable of handling transcendental equations, as long as all the roots being sought are simple.

But first: I slightly tidied up the definition of f[] (e.g. by using auxiliary variables for common subexpressions), as the original version brought tears to my sensitive eyes:

f[M_?NumericQ, λ_?NumericQ, ω_?NumericQ] := Module[{c, k, ξ1, ξ2},
c = Sqrt[M/λ]; k = I M/(2 λ);
ξ1 = (I - 1) (ω/Sqrt[λ] - Sqrt[λ]/2); ξ2 = (I - 1) (ω/Sqrt[λ] + Sqrt[λ]/2);
({1 + I, -I c}.ParabolicCylinderD[{-k, -k - 1}, ξ2]
{c, 1 + I}.ParabolicCylinderD[{k - 1, k}, I ξ1] -
{I c, 1 + I}.ParabolicCylinderD[{-k - 1, -k}, ξ1]
{c, -1 - I}.ParabolicCylinderD[{k - 1, k}, I ξ2])/2]


...and with that,

roots = FindAllCrossings[Im[f[2, 100, ω]], {ω, 0, 50}, WorkingPrecision -> 20]
{1.4210217375131208861, 4.9080677718060732317, 7.6276760758692264160,
11.242328271551264279, 14.025220377481373494, 17.188413671355074367,
20.686743750589305061, 23.568643080603343806, 26.490531437543067517,
29.849368653509459477, 33.222900929283185978, 36.429230282166527794,
39.466210718845193558, 42.459671861175573218, 45.512697669849073416,
48.625869297148536333}


As a graphical verification:

Plot[Im[f[2, 100, ω]], {ω, 0, 50},
Epilog -> {Red, AbsolutePointSize[4], Point[Thread[{roots, 0}]]},
Frame -> True, PlotStyle -> RGBColor[59/67, 11/18, 1/7]]


• Hi, thank you very much! I have made an extension to the problem. Please have a look. – Ambrose Chau Jun 26 '15 at 14:40
• As your addition is no longer about solving your transcendental equation, you probably should ask a new question. Link to this older question for context when you write your new question. – J. M.'s ennui Jun 26 '15 at 14:48

It is enormously faster to use

x = f /. {M -> 2, λ -> 100} // Simplify;
FindRoot[Im[x], {ω, 5}]


Then, given the space of roots for Im[x]

DeleteCases[Union[Table[ω /. FindRoot[Im[x], {ω, i}], {i, 50}],
SameTest -> (Abs[#1 - #2] < 10^-8 &)], z_ /; z < 0 || z > 50, Infinity]


finds all positive roots less than 50

(* {1.42102, 4.90807, 7.62768, 11.2423, 14.0252, 17.1884, 20.6867, 23.5686, 26.4905,
29.8494, 33.2229, 36.4292, 39.4662, 42.4597, 45.5127, 48.6259} *)


These results are consistent with a Plot of x.

Plot[{Re[x], Im[x]}, {ω, 0, 50}]


Introduction

For many transcendental functions, NSolve can solve for the roots, but not in this case. Since the roots are real we can apply the "Chebyshev-proxy rootfinder" method (CPR) which is based on the "colleague matrix" of a truncated Chebyshev series approximation to the function (see this answer by J.M. and the book by Boyd (2014)).

The first step is to form the Chebyshev series $f(\omega) = \sum a_n T_n(\zeta)$, where $\zeta = \frac{1}{2} (\omega (b-a)+(a+b))$ ranges over $-1 \le \zeta \le 1$ while $\omega$ ranges over $a \le \omega \le b$. We will take $[a,b]$ to be OP's example $[0,20]$. The series is truncated when the desired precision is achieved. By default, this is set to $10^{-14}$ in adaptiveChebSeries[]. For a function that is analytic in a complex neighborhood of $[a,b]$, the series will converge at least geometrically (or exponentially with an exponential rate of $1$). One can expect to get nearly machine precision for such a function. See Boyd (2014) for further discussion.

The next step is to get the eigenvalues of the colleague matrix. The roots will be the rescaled values $\omega_k$ in interval $[a,b]$ of the eigenvalues $\zeta_k$ in the interval $[-1,1]$. This is done by chebRoots[].

The success of this method depends primarily on the success of approximating $f(x)$. This has not been automated in adaptiveChebSeries[], so I'll show some hand checking at the end. One could also use this method to get a rough single-precision approximation of the roots and polish them off with FindRoot. (The option adaptiveChebSeries[..., PrecisionGoal->6] would accomplish this.)

Finding the roots

Code for the functions is given at the end of this answer and can also be found here.

h = f /. {M -> 2, λ -> 100} // Simplify;    (* Simplify seems sufficient *)
obj = Function @@ {ω, Im@h};                (* function equivalent of f[w] *)
dobj = Function @@ {ω, Im@D[h, ω]};         (* to be used in analyzing the roots *)

(a0 = 0; b0 = 20;                           (* specify interval *)
coeff = adaptiveChebSeries[obj, a0, b0];  (* approximate f(w) === obj[w] *)
roots = chebRoots[coeff, {a0, b0}]      (* get roots of approximant *)
) // AbsoluteTiming
Length[roots]  (* check how many roots found *)
(*
{1.57934,
{1.42102, 4.90807, 7.62768, 11.2423, 14.0252, 17.1884}}
6
*)


Checking the roots

We can see in the table below that the roots are calculated to roughly machine precision. (We approximate the error by a single Newton-Raphson step, which a fussier person could check with FindRoot at a higher precision than machine precision.)

The relative error shows the roots are almost all accurate to machine precision.

TableForm[Transpose@
{roots,
obj /@ roots,
(obj /@ roots/dobj /@ roots)/roots // Abs},
TableHeadings -> {None, {"root", "residual", "rel. err. in rt."}}
]


We can visualize the roots:

Plot showing the OP's function and roots (indicated by the red lines).

objPlot = Plot[obj[ω], {ω, 0, 20}];

Show[objPlot, GridLines -> {roots, None},
GridLinesStyle -> Directive[Red, Thickness[0.002]]]


Checking the Chebyshev approximant

We can see graphically that the approximant is close to the function.

Plot showing the OP's function (blue) and the Chebyshev approximant (light gray).

Show[
objPlot /. _AbsoluteThickness -> AbsoluteThickness[4.4],
Plot[obj2[ω], {ω, 0, 20}, PlotStyle -> LightGray]
]


Comparing the difference, the relative error shows the approximant is not quite down to machine precision. Except for near the roots where one expects the relative error to grow, the error stays below $10^{-14}$.

Plot showing the OP's function (blue) and the Chebyshev approximant (light gray).

LogPlot[(obj[ω] - obj2[ω])/obj2[ω] // Abs,
{ω, 0, 20},
PlotPoints -> 100, MaxRecursion -> 2,
Frame -> True, Axes -> False,
FrameLabel -> {HoldForm@ω, "relative error"},
GridLines -> {roots, None},
GridLinesStyle -> Directive[Red, Thickness[0.002]]]


Inspecting the Chebyshev coefficients, we can see that once $n >12$, the coefficients decrease roughly geometrically until $n=37$ when machine precision $\approx 10^{-16}$ is reached. After that, they bounce around at that level. (The function adaptiveChebSeries[] actually computed 65 coefficients and discarded about twenty using Boyd's heuristic with a cutoff of $10^{-14}$. Perhaps a few more should have been discarded.) One could recompute the coefficients at a higher working precision and see the geometric convergence continue.

Plot showing the the Chebyshev coefficients $a_n$. Once the degree of the approximation reaches roughly twice the number of roots (6), the series converges geometrically until it reaches machine precision. At that point rounding error dominates the coefficients.

ListLogPlot[Abs@coeff,
Frame -> True, Axes -> False,
FrameLabel -> {Row[{"Degree ", HoldForm@n, " of coefficient"}],
HoldForm@Abs[Subscript[a, n]]}]


Code dump

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
* *)
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 < 3 && # <= 2^OptionValue[MaxIterations] OptionValue["Points"] &
(* at least two coefficients dropped *)
];
If[len < 3,
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]
]

• This is great! A note: the type 1 DCT is conventionally done for lengths of the form $2^n +1$ (corresponding to a power-of-two FFT), so the traditional Clenshaw-Curtis hierarchy proceeds as $3, 5, 9, 17, 33, \dots$. Thankfully, Mathematica is not constrained to these lengths, so this works fine. – J. M.'s ennui May 6 '16 at 0:23
• @J.M. Probably I should rename "Points" -> n to "Order -> n", because it actually uses n+1 points. It iteratively doubles n until the desired precision is reached. So if you start with n = 2^k, it will always use 2^(k+m} + 1 points for m = 0, 1,...`. – Michael E2 May 6 '16 at 19:23