5
$\begingroup$

Here is the equation I want to solve. [$c$ and $M$ have specified values, and I want the value of $p$.]:

$$(2p-c)\int_p^M\frac1{n^2}e^{p-n}\,\mathrm dn=1-\frac{c}{p}$$

For example with $c=0.1$ and $M=10$, I can plot the function in Mathematica and see that the solution is about 0.735. If I write the code

c=1

M=50

NSolve[(2 p - c)*
NIntegrate[(1/n^2)*Exp[p - n], {n, p, M}] - (1 - c/p) == 0, {p, 0.75}]

When I run the code in Mathematica I get

nlim: "n = p is not a valid limit of integration. 

I get it. p is in the limit of integration. But if I go from NIntegrate to regular Integrate, then after a few minutes I get

This system cannot be solved with the methods available to NSolve. >>


NSolve[ConditionalExpression[-1 + 1/
p + (-1 + 2 p) (-(1/50) E^(-50 + p) + 1/p - 
   E^p ExpIntegralEi[-50] + E^p ExpIntegralEi[-p]) == 0, 0 < Re[p] < 50 && Im[p] == 0], p]

I appreciate that the function I'm using is connected to the Exponential Integral (or Upper Incomplete Gamma Function). That's why I need Mathematica's help. Here a graph shows that the solution is about 1.75.

Ultimately, I'm looking to get $p(c, M)$ -- $p$ as a function of $c$ and $M$. Actually, I'd be happy with graphing $p$ as a function of $c$ and $M$.

Note also that with some other values of $c$ and $M$, I get a solution using Integrate but there is a different error. For example, with $c=0.1$, $M=10$ I get

Warning: NSolve was unable to prove that the solution set found is complete. >>

{{p -> 0.735228}}

Any guidance for how I can find a numerical solution for $p(c,M)$?

$\endgroup$
5
  • $\begingroup$ NSolve[] is not really intended for use with transcendental equations like yours. Use FindRoot[], but make sure you have a good starting guess. $\endgroup$ Aug 5, 2016 at 21:27
  • 1
    $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Michael E2
    Aug 5, 2016 at 21:50
  • $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this meta Q&A helpful $\endgroup$
    – Michael E2
    Aug 5, 2016 at 21:51
  • 2
    $\begingroup$ J.M is a god. Thank you! $\endgroup$ Aug 5, 2016 at 22:02
  • $\begingroup$ Barry Nalebuff here? I've heard about you here in Brazil. $\endgroup$
    – LCarvalho
    Aug 5, 2016 at 22:48

4 Answers 4

2
$\begingroup$

Starting with a plot from 0 to 100, one can narrow down the root:

Block[{c = 1, M = 50},
 NSolve[(2 p - c)*
      Integrate[(1/n^2)*Exp[p - n], {n, p, M}, 
       Assumptions -> 0 < p] - (1 - c/p) == 0 && 0 < p < 5, p]
 ]
(*  {{p -> 1.75081}}  *)
$\endgroup$
2
$\begingroup$

If I do:

i[p_, M_] := NIntegrate[(1/x^2)*Exp[p - x], {x, p, M}]
f[y_, c_] := (2 y - c)*i[y, 50] - (1 - 1/y)
Plot[f[y, 1], {y, 0, 1}]

I get a nice plot.

FindRoot[{f[y, 1] == 0}, {y, 0.75}]

complains like you describe, but eventually gives

 {y->1.7508107}

I'm not sure why there are the complaints about "x = y is not a valid limit of integration" ...

$\endgroup$
2
$\begingroup$

With J.M's help the FindRoot did the trick. Note that if one integrates by parts the result is the same as the ExpIntegralEi function. And now the solution is instantaneous:

FindRoot[Exp[-t]/(-ExpIntegralEi[-t] + ExpIntegralEi[-M]) == 2 t - c, {t, c + 1}]
$\endgroup$
1
  • $\begingroup$ FWIW: -ExpIntegralEi[-t] is the same as ExpIntegralE[1, t]. $\endgroup$ Aug 6, 2016 at 4:15
2
$\begingroup$

Actually, I'd be happy with graphing $p$ as a function of $c$ and $M$.

In that case, you can just use ContourPlot3D[] to generate your desired plot.

ContourPlot3D[(2 p - c) Exp[p] Gamma[-1, p, M] == 1 - c/p,
              {c, 0, 30}, {M, 0, 30}, {p, 0, 40}]

plot of implicit surface

(Note that I have taken the opportunity to use the three-argument form of the incomplete gamma function, which Mathematica has support for.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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