Find point at which equation stops having roots (if it exists)

Question

I am interested in the roots of this function:

f[M_, b_] := 1 - (2 M Gamma[2, 0, (1/M + b M)/Sqrt[b]])/(1/M + b M)

for fixed values of b. In particular I want the value of b at which the function of M stops having roots.

My code to do this is:

Do[s = FindRoot[f[m, b], {m, n}][[1, 2]]; Print[s]; n = s + .5;, {b, .005, 1000, 5}]

It seems to work but I'm not sure how to keep trying new values of b until there are no more solutions, instead of trying bigger and bigger b_max (in this case 1000).

Also I get some error messages from FindRoot; is this something to worry about? How accurate is this code?

Thanks!

Answer 1

It's best to try for an analytic solution rather than a numerical one. Just because you can't find a root numerically, doesn't mean there isn't one. And if you do find a root, it helps to check it by throwing it back into your original equation. You might not have noticed that Mathematica has suppressed error messages so as not to spam you.

We can rewrite your expression without the gamma function by using the identity:

grule = Gamma[2, 0, x_] -> 1 - (1 + x)/E^x

Mathematica won't simplify this by itself, but there is enough information in the help about Gamma to work this out.

If we are just concerned with roots, we only care if the numerator of a rational function is zero. So we can change your function as so:

h = FullSimplify@Numerator@Together[f[m, b] /. grule]

$\sqrt{b} e^{\frac{b m^2+1}{\sqrt{b} m}} \left((b-2) m^2+1\right)+2 m \left(b m^2+\sqrt{b} m+1\right)$

Now if $m > 0$ and $b > 2$ then this expression is always positive and hence has no roots. You can even get Mathematica to confirm this:

Reduce[{h > 0, b > 2, m > 0}] (* b > 2 && m > 0 *)

This is a reasonably robust answer to your initial problem. If you want to investigate further, then you can use Mathematica to gather empirical evidence as to what's happening. Plotting the function is always a good start. After playing around with values and PlotRanges, you might notice:

Plot[Evaluate[h /. b -> {5, 10}], {m, -0.5, 0}, PlotRange -> {-1, 1}]

So somewhere between $b = 5$ and $b = 10$ the peak of the curve will hit the axis. Let's find this numerically:

FindRoot[{h == 0, D[h, m] == 0}, {m, -0.3}, {b, 7}]

{m -> -0.256608, b -> 6.88164}

If you want a robust proof, then you might need to ask a mathematician. If eyeballing a graph is good enough, then ask Mathematica.