# Integrating a difficult function (Numerically ?)

I've been trying to make the following integral but I get no results :

Integrate[ (3/(16 + r^2)- 4/(6.25 + r^2)) E^(-2r^2) Sin[k r] r, {r, 0, Infinity}]

I am new in Mathematica and I have tried with the function Integrate and NIntegrate. Is there any command so that I can get a result of the above integral? Thanks :D

• try Nintegrate but you need to specify numerical value for k Commented Aug 17, 2017 at 6:26
• This does not look like it has a closed form known to Mathematica. But, if you specify a value for k, NIntegrate[] should be able to handle it. Commented Aug 17, 2017 at 6:41

The exact integration fails in this case. One can do this numerically. As much as I understand, you need to have the result as a function of k. If you do this numerically, the result may only be an approximate function. If this is what you agree to have, the following approach may be undertaken.

First, note that because of the factor Sin[k r] your integral is highly oscillating at large values of k. In this case, it can be easily circumvented by the replacement k r->R. Let us introduce your expression staying under the integral as

Clear[expr];
expr[k_, R_] :=
Simplify /@ (((3/(16 + r^2) - 4/(6.25 + r^2)) Exp[(-2 r^2)] Sin[k r]*
r) /. r -> R/k // Expand)


The expression takes the following form:

Now let us integrate, not forgetting to multiply the expression by 1/k coming from dr=dR/k:

lst = Table[{k,
NIntegrate[expr[k, R]*1/k, {R, 0, \[Infinity]}]}, {k,
0.1, 10, 0.1}]


giving you a list of pairs with the structure {k, int[k]}. You can plot it, if you wish.

Further, if you need to have an approximate analytical function, in the end, you can fit this result to some reasonable function:

model = -a*k*Exp[-c*k^2];
ff = FindFit[lst, model, {a, c}, k]
Show[{
ListPlot[lst, PlotStyle -> Blue,
AxesLabel -> {Style["k", Italic, 16], Style["J(k)", Italic, 16]}],
Plot[model /. ff, {k, 0.1, 10}, PlotStyle -> Red]
}]
(*  {a -> 0.0619, c -> 0.114} *)


where the plot is drawn to control the fitting quality:

Here the blue points show the results of the numerical integration, while the red solid line that of the fitting.

However, I guess you may try to do the integral analytically, but by hand. To do this, integrate the terms of the expression above separately. Replace in each of them Sin[R] by Exp[I R] and then, observing that the both integrands are even, integrate them from -infinity to +infinity and divide by 2, rather than from 0 to +infinity. Further, use the method of residuals for calculating Fourier integrals. You may have a look here for the method description: V. I. Smirnov, A Course of Higher Mathematics., 16 ed. (Pergamon Pr., Oxford, 1964), Vol. 3 Chapter 3, or somewhere else. In the end take the imaginary part of the solution.

Have fun!

• Thank you very much for your help Alexei, I will check it and reproduce it. As you correctly infer I need the result of the integral as a function of k. If it's not much to ask, how did you think in the function to fit ?. I need to make the next integral too: Integrate[ (2k)/(1 + E^(4 (r - 1))) Sin[k r] r, {r, 0, Infinity}]. Should I use the same process you did ? Thank you very much for your help :D Commented Aug 17, 2017 at 18:31
• @Jhoan Perez It is not much to ask, but I do not have much to answer. I just guessed it and tried. The more that it is easy to vary the model a bit and see what it gives. The third or fourth attempt gave a good fitting. With your new function, the process might be the same, but I would first try to get the answer as Integrate, and only if it does not work I would apply the procedure I proposed above. Concerning the approximation of your integral, everything depends on how the ListPlot[lst] looks like. If it is close to the dependence shown above, probably the same approximation is OK. Commented Aug 18, 2017 at 7:02
• Continuation. May be with some slight variations. If not, it is exactly the place, where the science turns into art. I wish you success. Commented Aug 18, 2017 at 7:03

You can do this integral analytically. One method is to recast the integral as an ODE, and then solve the ODE. Let's consider the integral:

int[k_] := Inactive[Integrate][(r Exp[-2 r^2] Sin[k r])/(s^2+r^2), {r, 0, Infinity}]


The trick is to realize that:

LHS = s^2 int[k] - int''[k]


s^2 Inactive[Integrate][(E^(-2 r^2) r Sin[k r])/( r^2 + s^2), {r, 0, [Infinity]}] - Inactive[Integrate][-((E^(-2 r^2) r^3 Sin[k r])/(r^2 + s^2)), {r, 0, [Infinity]}]

can be simplified as follows:

e = Replace[
LHS,
a_ h_[i1_, lim_] + b_ h_[i2_, lim_] :> h[a i1 + b i2, lim]
] //Simplify


Inactive[Integrate][E^(-2 r^2) r Sin[k r], {r, 0, [Infinity]}]

and this integral is one that Mathematica can handle:

RHS = e //Activate


1/8 E^(-(k^2/8)) k Sqrt[π/2]

Finally, we need to know 2 initial conditions:

i0 = int[0] //Activate
i0p = Assuming[s>0, int'[0] //Activate]


0

1/4 (Sqrt[2 π] - 2 E^(2 s^2) π s Erfc[Sqrt[2] s])

Putting the above pieces together, we have the following ODE to solve:

soln = DSolve[
{
s^2 i[k] - i''[k] == RHS,
i[0] == 0,
i'[0] == i0p
},
i,
k
]

res = FullSimplify[i[k] /. First @ soln, k>0 && s>0]


{{i -> Function[{k}, 1/4 E^(-k s + 2 s^2) π (Erf[(k - 4 s)/(2 Sqrt[2])] + Erf[Sqrt[2] s] - E^(2 k s) Erf[Sqrt[2] s] + E^(2 k s) Erf[(k + 4 s)/(2 Sqrt[2])] + Erfc[Sqrt[2] s] - E^(2 k s) Erfc[Sqrt[2] s])]}}

1/4 E^(-(k - 2 s) s) π (1 + Erf[(k - 4 s)/(2 Sqrt[2])] - E^(2 k s) Erfc[(k + 4 s)/(2 Sqrt[2])])

Summarizing, we have:

integral[s_, k_] = res


1/4 E^(-(k - 2 s) s) π (1 + Erf[(k - 4 s)/(2 Sqrt[2])] - E^(2 k s) Erfc[(k + 4 s)/(2 Sqrt[2])])

or in TeX:

$$\int _0^{\infty }\frac{e^{-2 r^2} r \sin (k r)}{r^2+s^2}dr=\frac{1}{4} \pi e^{s (-(k-2 s))} \left(\operatorname{erf}\left(\frac{k-4 s}{2 \sqrt{2}}\right)-e^{2 k s} \operatorname{erfc}\left(\frac{k+4 s}{2 \sqrt{2}}\right)+1\right)$$

The integral in question is:

3 integral[4, k] - 4 integral[2.5, k]


-E^(-2.5 (-5. + k)) π (1 + Erf[(-10. + k)/(2 Sqrt[2])] - E^(5. k) Erfc[(10. + k)/(2 Sqrt[2])]) + 3/4 E^(-4 (-8 + k)) π (1 + Erf[(-16 + k)/(2 Sqrt[2])] - E^(8 k) Erfc[(16 + k)/(2 Sqrt[2])])

and the plot is:

Plot[3 integral[4, k] - 4 integral[2.5, k], {k, 0, 10}, WorkingPrecision->30]


in agreement with the version given by @AlexeiBoulbitch