# Same integral yielding to different results

I am currently working with the following integrals

$$$$\int_{0}^{\infty} dk\thinspace \frac{k^{3}e^{-2kd}}{\omega^{2}+k^{4}} = \frac{1}{\omega^{2}d^{4}}\int_{0}^{\infty}d\epsilon\thinspace\frac{\epsilon^{3}e^{-2\epsilon}}{1+\frac{\epsilon^{4}}{\omega^{2}d^{4}}}$$$$ which can be verified by the use of $$\epsilon = k d$$.

The problem is: depending on which side of the above equation I use, I obtain different results for its value as a function of $$d$$.

Here is how I define the integrals

SEND[\[Omega]_, d_] := NIntegrate[k*((k^2*Exp[-2*k*d])/((k^2)^2 + \[Omega]^2)), {k,0,Infinity}, AccuracyGoal -> 20, PrecisionGoal -> 20, MaxRecursion -> 20,   WorkingPrecision -> 50]

SEdimN[\[Omega]_, d_] := (1/(\[Omega]^2*d^4))*NIntegrate[(\[Epsilon]^3*Exp[-2*\[Epsilon]])/(1 + (1/(\[Omega]^2*d^4))*\[Epsilon]^4), {\[Epsilon], 0, Infinity},AccuracyGoal -> 20, PrecisionGoal -> 20,MaxRecursion -> 20, WorkingPrecision -> 50]


and here is how I plot them

LogLogPlot[{SEND[1, z], SEdimN[1, z]}, {z, 10^-10, 10^5}, LabelStyle -> Black, PlotLegends -> "Expressions", PlotRange -> {10^-6, 10^3}, PlotPoints -> 1000,WorkingPrecision -> 50]


The result I obtain is something like

from where we can clearly see that for most of the range both integrals yield the same value. So my question is: Whey they are yielding different results? How do I know which one is the correct one!?

• This integral can be calculated symbolically Integrate[(k^3 Exp[-2 d k])/(\[Omega]^2 + k^4) , {k, 0, \[Infinity]}, Assumptions -> d > 0 && \[Omega] > 0]. The shape of the blue curve suggests that the numerical integral is calculated with insufficient accuracy. For related issues see e.g. Numerical underflow for a scaled error function. Commented Sep 28, 2021 at 2:41
• The jump around $10^{-4}$ strongly suggests a numerical issue. Commented Sep 28, 2021 at 3:14

Better precision control is required.

Clear["Global*"]

SEND[ωv_?NumericQ, dv_?NumericQ] :=
Module[
{ω = SetPrecision[ωv, 60], d = SetPrecision[dv, 60]},
NIntegrate[
k*((k^2*Exp[-2*k*d])/((k^2)^2 + ω^2)),
{k, 0, Infinity},
MaxRecursion -> 100,
WorkingPrecision -> 50]]

SEdimN[ωv_?NumericQ, dv_?NumericQ] :=
Module[
{ω = SetPrecision[ωv, 60], d = SetPrecision[dv, 60]},
(1/(ω^2*d^4))*
NIntegrate[
(ϵ^3*
Exp[-2*ϵ])/(1 + (1/(ω^2*d^4))*ϵ^4),
{ϵ, 0, Infinity},
MaxRecursion -> 100,
WorkingPrecision -> 50]]

{SEND[1, z], SEdimN[1, z]} /. z -> 10^-7

(* {14.847733027640972818058721332928476020439557891260,
14.847733027640972818058721332928476020439557891260} *)

LogLogPlot[
{SEND[1, z], SEdimN[1, z]},
{z, 10^-10, 100},
LabelStyle -> Black,
PlotStyle -> {Automatic, Dashed},
PlotLegends -> Placed["Expressions", {0.35, 0.4}],
PlotRange -> {10^-6, 10^2},
PlotPoints -> 75,
MaxRecursion -> 5,
WorkingPrecision -> 50]


Remove the AccuracyGoal -> 20 from SEdimN[]` or from both (from both shown):