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
* *)
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]
]
f [ω_]rather thanf? You could improve your question by providing the code to calculate at least one root. $\endgroup$ – bbgodfrey Jun 25 '15 at 18:30