Skip to main content

equation solving - 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.



Answer




...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]]

roots of a linear combination of PCFs



Comments

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...

list manipulation - Selecting multiple columns from a matrix?

Sample data: data = { {{2013, 1, 1}, 24.13, 167.67, 231.82}, {{2013, 1, 2}, 32.15, 170.92, 225.99}, {{2013, 1, 3}, 35.43, 172.68, 221.67}, {{2013, 1, 4}, 36.73, 173.05, 218.32}, {{2013, 1, 5}, 58.19, 165.96, 197.05}, {{2013, 1, 6}, 69.99, 163.50, 187.52}, {{2013, 1, 7}, 71.37, 154.21, 175.58}, {{2013, 1, 8}, 72.51, 149.66, 163.25}}; I want a DateListPlot with three graphs, so for a matrix formed by columns 1 and 2, one for columns 1 and 3, and 1 for columns 1 and 4. At the moment I'm using this code: data2 = Transpose[{data[[All, 1]], data[[All, 2]]}]; data3 = Transpose[{data[[All, 1]], data[[All, 3]]}]; data4 = Transpose[{data[[All, 1]], data[[All, 4]]}]; DateListPlot[{data2, data3, data4}, Joined -> True, Filling -> {3 -> {1}}] but I have a hunch that this can be done more efficiently. I don't like the Transpose s in particular. Any ideas? edit (for extra credit) What if I need to multiply the second column by 2, which in my solution is simp...

calculus and analysis - Integral of the Sinc product

Let us consider the following integral $$ B_n = \int_0^\infty \prod_{k=1,3,5,\dots}^n\frac{\sin (x/k)}{x/k}dx $$ By definition, Sinc[x] == Sin[x]/x , therefore B[n_?OddQ] := Integrate[Product[Sinc[x/k], {k, 1, n, 2}], {x, 0, Infinity}] B /@ Range[1, 13, 2] {Pi/2, Pi/2, Pi/2, Pi/2, Pi/2, Pi/2, Pi/2} OK, everything is fine. But... B[15] 467807924713440738696537864469 Pi/935615849440640907310521750000 What's going on? Answer As Eckhard wrote in comments B[n] is the n-th Borwein integral . (The letter B was not accidental :) ) This funny properties of Borwein integrals is related to the Fourier transform of Sinc function FourierTransform[Sinc[x], x, k] 1/2 Sqrt[Pi/2] (Sign[1 - k] + Sign[1 + k]) Plot[%, {k, -2, 2}, Filling -> 0] which is the box function. The result is $\pi/2$ while the sum $$ 1/3+1/5+\dots+1/n $$ B_n = \frac{\pi}{2} - \pi \bigg(\sum_{k=3,5,\ldots}^n\frac{1}{k} -1\biggr)^\frac{n-1}{2}\prod_{k=3,5,\dots}^n\frac{k}{k-1}. $$ For $n=15$ it is equal to $$ \frac{46780...