Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
    
By the way, the error message I keep getting is: The line search decreased the step size to within tolerance specified \ by AccuracyGoal and PrecisionGoal but was unable to find a sufficient \ decrease in the merit function. You may need more than \ MachinePrecision digits of working precision to meet these tolerances –  misi Apr 17 '13 at 20:13
    
That's the message indicating that there probably isn't a root. Also, I vote for b = 2. –  wxffles Apr 17 '13 at 20:14
    
@wxffles So can I be confident that there is no solution after b~2 and that this is not a precision error? –  misi Apr 17 '13 at 20:18
    
I presume you've seen this already? –  J. M. Apr 17 '13 at 20:31
    
Actually you can have roots at negative M values for b > 2. –  wxffles Apr 17 '13 at 20:56

1 Answer 1

up vote 3 down vote accepted

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

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.

share|improve this answer
    
Great, that's really helpful. Thank you! By the way, what's your source for the gamma function identity? Just Mathematica Help? –  misi Apr 18 '13 at 16:05
    
The help tells you how to go from the generalised incomplete gamma function to the incomplete gamma function. Then you could just evaluate the integral directly to get the rest. Or many pages (such as mathworld) give the identity for integer arguments. –  wxffles Apr 18 '13 at 20:48

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.