I am trying to solve the following equation numerically:

Equation to solve for $y$:

$$na(ay-y)^{n}\int_{ay}^{+\infty} \frac{(x-y)^{-n}}{x-ay}dx=b$$

with for example $a=4$, $n=1.25$ and $b=1.6$.

I have adapted an answer to a similar question that I found on StackExchange:

Pdf[x_, y_, a_, n_] = n (x - y)^(-1 - n) ((-1 + a) y)^n
eps[x_, y_, a_] = (x - y)*a/(x - y*a)
Func[x_, y_, a_, n_] = FullSimplify[eps[x, y, a]*Pdf[x, y, a, n]]
f[y_?NumericQ] := 
  NIntegrate[Func[x, y, 4, 1.25], {x, 4*y, Infinity}];
FindRoot[f[y] == 1.6, {y, 0}, WorkingPrecision -> 10000]

But I am getting a list of more than 10 error messages and the solution I get does not seem to be correct. My guess (from looking at plots of the integral NIntegrate[Func[x, y, 4, 1.25], {x, 4*y, Infinity}] which also generate loads of error messages) is that the solution is very small, that's why I am trying to search around 0. But I have no idea how to find it. Can anyone help please?


The integral is calculated exactly

f= Integrate[Func[x, y, a, n], x]

(*Out[]= (a n (x - y)^(
 1 - n) ((-1 + a) y)^n Hypergeometric2F1[1, 1 - n, 2 - n, (
  x - y)/((-1 + a) y)])/((-1 + a) (-1 + n) y)*)

But the integral diverges at the lower limit x=4ywith n=125/100,a=4. Check

f/. {x -> 4 y, a -> 4, n -> 125/100}

(*Out[]= -\[Infinity]*)
  • $\begingroup$ Thank you very much! One quick question regarding your answer: when I use your code for the limit evaluation in my Mathematica, I get two error messages (Power::infy: Infinite expression 1/0 encountered) and then Out[]=ComplexInfinity. Is that because I may be using a different version of Mathematica or am I forgetting to load some particular package...? Thanks again! $\endgroup$ – Lednacek Oct 15 '19 at 12:28
  • $\begingroup$ I have 12.0.0 for Microsoft Windows (64-bit) (April 6, 2019). What do you have? $\endgroup$ – Alex Trounev Oct 15 '19 at 12:33
  • $\begingroup$ 11.3.0 for Mac OS (32-bit). I guess that explains it. Thanks again! $\endgroup$ – Lednacek Oct 15 '19 at 12:34
  • $\begingroup$ @Lednacek You're welcome! $\endgroup$ – Alex Trounev Oct 15 '19 at 12:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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