Finding the roots of Hypergeometric1F1[]

Question

I am trying to find the roots, λ, for this equation:

Hypergeometric1F1[1/4 (2 -  λ /β), n + 1, β] 

for certain β and n. Here is my Mathematica code.

eq[n_, β_, λ_] = Hypergeometric1F1[1/4 (2 - λ/β), n + 1, β]

Find the root near λ = β.

ED[n_, β_] := λ /. 
  FindRoot[eq[n, β, λ] == 0, {λ, β}] 

List all λ-values for each value of n when β = 0.00001

{ED[0, 0.00001], ED[1, 0.00001], ED[2, 0.00001], ED[3, 0.00001], ED[4, 0.00001],
 ED[5, 0.00001], ED[6, 0.00001], ED[7, 0.00001], ED[8, 0.00001], ED[9, 0.00001],
 ED[10, 0.00001]}

Then, I got {5.78319, 14.682, 26.3746, 40.7064, 57.5829, 76.9388, 98.7262, 122.907, 149.453, 178.337, 209.54}

However I'm supposed to get {5.78306, 14.6819, 26.3744, 30.4715, 40.707, 49.2186, 57.5823, 70.8493, 74.8865, 76.9392, 95.2771}

I guess Mathematica gives only one root. How can I also find the 2nd, 3rd, 4th roots, for each value of n?

Answer 1

I can see two possible improvements. As expected with a subject like this, considering these improvements require knowing a little something about the properties of the Kummer confluent hypergeometric function. Thankfully, such knowledge is easily available these days, at sites like the DLMF.

Here they are:

1) One can use the Kummer transformation

$${}_1 F_1\left({{a}\atop{b}}\mid z\right)=\exp(z)\;{}_1 F_1\left({{b-a}\atop{b}}\mid -z\right)$$

to factor out exponential behavior in β (which only makes the task of finding zeroes on the real axis more difficult). That is,

Hypergeometric1F1[(2 -  λ/β)/4, n + 1, β] 

is the same as

Exp[β] Hypergeometric1F1[n + 1 - (2 - λ/β)/4, n + 1, -β]

and one can remove the Exp[β] factor for root-finding purposes.

Thus, one can write a scaled version of the function whose roots are being sought, like so:

eqScaled[n_, β_, λ_] := Hypergeometric1F1[n + 1 - (2 - λ/β)/4, n + 1, -β]

2) Using the asymptotic behavior of Kummer's function for large values of the upper parameter, one can easily develop a better starting approximation for FindRoot[] to polish:

ED[n_, β_, k_Integer:1] := Block[{λ}, λ /. FindRoot[eqScaled[n, β, λ] == 0,
                                                    {λ, BesselJZero[n, k]^2 - 2 n β}]]

k is an index for the roots, e.g. k = 1 for the first root.

We can now try Vitaliy's example:

roots = Table[ED[1, 1, k], {k, 10}]
{13.0143, 47.552, 101.833, 175.854, 269.615, 383.115, 516.355, 669.334, 842.052, 1034.51}

Plot[eqScaled[1, 1, λ], {λ, 0, 1310}, 
 Epilog -> {AbsolutePointSize[4], Red, Point[{#, 0}] & /@ roots}]

Answer 2

Answer

Here is an alternative approach, leveraging Mathematica's intelligent plotting algorithms. The basic idea is to use the list of points evaluated when plotting the function.

Suppose I want the roots between 0 and 1400. Plot the function and assign the plot to a variable (I won't display the plot here):

eq[n_, β_, λ_] = Hypergeometric1F1[1/4 (2 - λ/β), n + 1, β];
p = Plot[eq[1, 1, λ], {λ, 0, 1400}];

Extract the actual list of points from the plot:

points = p[[1, 1, 3, 3, 1]];

This looks like this:

points // Shallow

so it is the actual list of points used by Plot. You could for instance ListLinePlot them and get something like Plot would give.

Now, the idea is to inspect this list, note when the $y$ coordinate changes sign between neighbouring points, and use the $x$ coordinate as the initial value for FindRoot. Clearly, this does not guarantee that we'll get the root nearest that point, but it usually works.

Here is how to locate where neighbouring points change sign: First pair up the neighbours:

pairs = Partition[points[[All, 2]], 2, 1];

then here is where the signs of neighbouring $y$-coordinates change:

flips = Position[Sign[Times @@ # & /@ pairs],-1];

and now find the roots:

roots=λ /. FindRoot[eq[1, 1, λ],
   {λ,
    First@First@points[[#]],
    First@First@points[[1 + #]]
    }] & /@ flips

And it works:

Show[
 p,
 Graphics@{Red, PointSize[Large], 
   Point[Partition[Riffle[roots, 0], 2]]}
 ]

(the last root is correctly calculated, but I was too lazy to plot it).

Comments

We didn't have to use Plot to extract the points; we could have just constructed a list of values at equidistant points using Table, say. I chose Plot because sometimes its adaptive stepsizes can be useful.

This method has the advantage that it requires less detailed knowledge than e.g. J.M.'s approach, and it may be more efficient than Vitaliy's method.

Please note that this method is also described here and, apparently, in a book by Stan Wagon (also mentioned in the answer linked to).

Answer 3

Kummer confluent hypergeometric function Hypergeometric1F1 has many zeros. Using your definitions:

Plot[eq[1, 1, λ], {λ, 0, 1310}, PlotRange -> {-.25, .15}, Filling -> 0]

You see that the 1st four roots are less than 200, give it a list of initial values in that range - large enough to get all the roots:

Union@Round[ED[1, 1, Range[1, 200, 5]], 10.^(-6)]

{13.0143, 47.552, 101.833, 175.854}

Some of the roots will be the same maybe with a small difference in the very last digit. This is why you need Union and Round.

Answer 4

Yet another approach to locating zeros of oscillatory functions is to find extrema via MedianFilter, then to useFindRoot with two adjacent extrema as starting values. I used this method to find roots of RiemannSiegelZ[t] before I learned about Gram points of the zeta function, and before Mathematica V6 came out with ZetaZero[k].

It requires a grid of function values $v$, so the method only works if the function is cheap to calculate. The grid may be linearly spaced, log spaced, or intelligently spaced via Plot as in the current answer from acl and the earlier one from J. M.

Next, run MedianFilter[v,1]. The function has a maximum at a point where its value exceeds the median of the three values centred on that point. A minimum occurs where the function value is less than the median. On slopes, the median equals the current, central point. Pick locations where the function value and its median are unequal, these are the extrema.

Finally, pass two adjacent extrema to FindRoot:

Block[{r=Range[2.,550.,5.], v, p, t, b=1, n=1},
       v = Transpose[{(v=Hypergeometric1F1[1/4(2-r/b),n+1,b]), MedianFilter[v,1]}];
       p = Pick[r, Apply[Unequal, v, 1]];
   Map[FindRoot[Hypergeometric1F1[1/4(2-t/b),n+1,b], Flatten[{t,#}]]&, Partition[p,2,1]]
]

{{t -> 13.0143}, {t -> 47.552}, {t -> 101.833}, {t -> 175.854},
 {t -> 269.615}, {t -> 383.115}, {t -> 516.355}}