Integration analitycally and numerically

Question

I've been trying to solve the next integral :

$$\int_{0}^{\infty}{\frac{652rSin[kr]}{72248.9k+ke^{1.85185r}}dr}$$

But if I do it analytically , I obtain a function which I can not work with:

ConditionalExpression[ 1/k (((0. + 330.575 I) HypergeometricPFQ[{1., 1. - (0. + 0.54 I) k, 
     1. - (0. + 0.54 I) k}, {2. - (0. + 0.54 I) k, 
     2. - (0. + 0.54 I) k}, -72248.9])/((0. + 1.85185 I) + 
   1. k)^2 - ((0. + 330.575 I) HypergeometricPFQ[{1., 
     1. + (0. + 0.54 I) k, 
     1. + (0. + 0.54 I) k}, {2. + (0. + 0.54 I) k, 
     2. + (0. + 0.54 I) k}, -72248.9])/((0. + 1.85185 I) - 
   1. k)^2), Abs[Im[k]] < 1.85185]

And if I try to solve it numerically I got this :

"The relative error 7.769066430156647` is larger than expected for \the integrand (6611.51\(0. +r)\Sin[0.1\r])/(72248.9 +E^(1.85185 (0. \+r))) over {0,\[Infinity]} with DoubleExponentialOscillatory method \and automatic tuning parameters, TuningParameters -> {10,5}. The \integration will proceed with TuningParameters -> {1,5}. "

Several times!!! and I cannot fit a good function. The ideal thing here ( I guess is the goal) is to find a function because I have to put that function inside another integral but I don't know how to accomplish that. I really need help on this :D. Thank you in advance.

Answer 1

For this integral, you can use a trick where you replace a parameter with an exact numeric quantity, and then replace the exact numeric quantity with a variable. For example, lets replace your 72247.9 with Pi and 1.85185 with EulerGamma:

sint = Integrate[
    (r Sin[k r])/(Pi k + k Exp[EulerGamma r]),
    {r,0,Infinity},
    Assumptions -> k>0
]

(2 EulerGamma)/(EulerGamma^2 + k^2)^2 + ( I π (LerchPhi[-π, 2, 2 - (I k)/EulerGamma] - LerchPhi[-π, 2, 2 + (I k)/EulerGamma]))/(2 EulerGamma^2 k)

(it is possible to evaluate the integral using s and t directly, but it takes a much longer time). Now, we need to replace Pi with s and EulerGamma with t:

sint /. {EulerGamma -> t, Pi -> s}

(2 t)/(k^2 + t^2)^2 + ( I s (LerchPhi[-s, 2, 2 - (I k)/t] - LerchPhi[-s, 2, 2 + (I k)/t]))/(2 k t^2)

Based on the above, it's reasonable to assume that the integral is:

int[s_, t_, k_] = sint /. {EulerGamma -> t, Pi -> s};

Or, in TeX:

$$\int _0^{\infty }\frac{r \sin (k r)}{k e^{r t}+k s}dr=\frac{2 t}{\left(k^2+t^2\right)^2}+\frac{i s \left(\Phi \left(-s,2,2-\frac{i k}{t}\right)-\Phi \left(-s,2,\frac{i k}{t}+2\right)\right)}{2 k t^2}$$

The NIntegrate version is:

Options[nint] = {WorkingPrecision->50};
nint[s_, t_, k_, OptionsPattern[]] := NIntegrate[
    (r Sin[k r])/(s k + k Exp[t r]),
    {r,0,Infinity}, 
    WorkingPrecision->OptionValue[WorkingPrecision]
]

Here are some plots comparing int with nint:

Plot[int[1000, t, 2], {t, 0, 1}, WorkingPrecision->50,PlotRange->All]

Plot[nint[1000, t, 2], {t, 0, 1}, PlotRange->All, WorkingPrecision->50]

NIntegrate::precw: The precision of the argument function ((r Sin[2 r])/(2000+2 E^(0.0000204286 r))) is less than WorkingPrecision (50.`).

NIntegrate::deodiv: DoubleExponentialOscillatory returns a finite integral estimate, but the integral might be divergent.

These plots take a long time to evaluate, so I will only show one more. Here is a plot of the function at the point of interest for the OP:

Plot[int[72248.9`50, 1.85185`50, k], {k, 0, 2}, WorkingPrecision->40]

Plot[nint[72249.9`50, 1.85185`50, k], {k, 0, 2}, WorkingPrecision->40]

NIntegrate::precw: The precision of the argument function ((r Sin[0.0000408571 r])/(2.95192 +0.0000408571 E^(1.8518500000000000000000000000000000000000000000000 r))) is less than WorkingPrecision (50.`).

NIntegrate::deorela: The relative error 1.6273881429326276607971423877290476125687271847202`50. is larger than expected for the integrand (r Sin[0.000040857142857142848336223645011600069665291812270880 r])/(2.9519244857142847848763267393223941326141357421875+0.000040857142857142848336223645011600069665291812270880 2.7182818284590452353602874713526624977572470937000^(1.8518500000000000000000000000000000000000000000000 r)) over {0,[Infinity]} with DoubleExponentialOscillatory method and automatic tuning parameters, TuningParameters -> {10,5}. The integration will proceed with TuningParameters -> {1,5}.

Answer 2

This is your expression:

expr1 = (652*r*Sin[k r])/(72248.9*k + k*Exp[1.85*r]);

Let us replace k r->R so that there is Sin[R] instead of Sin[k r]:

expr2 = expr1 /. r -> R/k // Simplify

(*   (652. R Sin[R])/((7.22*10^4 + E^((1.85 R)/k)) k^2)  *)

Now one can integrate over R, not forget to divide once more by k. Try to use the Method->"LevinRule" option:

lst = Table[{k, 
    1/k^3 NIntegrate[(
      652.` R Sin[R])/(72248.9` + E^((1.85` R)/k)) , {R, 
       0, \[Infinity]}, Method -> "LevinRule"]}, {k, 
    Join[Range[0.01, 0.99, 0.01], Range[1, 5, 0.05]]}];

yielding this:

ListPlot[lst]

Have fun!