# Finding the roots of Hypergeometric1F1[]

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?

-

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

-
+1 Nice idea to use asymptotics – Vitaliy Kaurov Aug 3 '12 at 4:14
Can we get the plot with λ and β for each n, (n=0,...,10)? The y-intercepts should be {5.78306, 14.6819, 26.3744, 30.4715, 40.707, 49.2186, 57.5823, 70.8493, 74.8865, 76.9392, 95.2771}. – Sony Aug 3 '12 at 21:25
@Sony, "should be" - how did you compute these values, first of all? Where did these come from? – J. M. Aug 3 '12 at 21:50
These are the values when β approaches to 0 which are the same as the (BesselJZero[n,k])^2. For example, when n=0 & β=0.0001, we get λ=5.78306. – Sony Aug 3 '12 at 22:00
@Sony: Ah, so that's why the asymptotics worked particularly well here (see, you could have mentioned things like those in your question!); still, I'm wondering where you're getting 30.4715; ED[3, 10^-4] gives 40.70586582284382, and N[BesselJZero[3, 1]^2] gives 40.70646581820033. So, where did these values come from? Anyway, do check the output of Table[BesselJZero[n, 1]^2, {n, 0, 10}] // N, and you'll see that some of your values could not have come from that. – J. M. Aug 3 '12 at 22:06

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

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

-
It's a bit like what I did here, yes? :) Still, your solution can benefit greatly from using the Kummer-transformed version of the function. – J. M. Aug 4 '12 at 10:17
@J.M. oops, I did not remember that! (I did remember seeing this somewhere but not where). I've added a link! About the transformation, I agree, but I try to avoid specializing things unless absolutely necessary and, in this case, it's not; so I left it general. Plus, it's well described in your post. – acl Aug 4 '12 at 11:05
@J.M. as an aside, seeing that transformation reminds me of a few days of pain resulting from having to manipulate a complicated combination of 2F1 functions a few years ago. In the event, it took more computer time to explicitly evaluate the analytical solutions (which I obtained in the form of series) than to do numerical simulations of the system. That taught me a lesson! – acl Aug 4 '12 at 11:08
Heh. Still, sometimes analytical manipulations can yield a form that is eminently more suitable for numerics, so it still pays to try out an identity or two... – J. M. Aug 4 '12 at 11:12
@J.M. yes, that is pretty much the bulk of my working day, some days. – acl Aug 4 '12 at 11:13

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.

-
By using a better starting value, it is possible to get the needed roots. See my answer. – J. M. Aug 3 '12 at 3:49

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}}
-