# Trouble with numerical evaluation of a four-fold integral

I am having trouble using NIntegrate in the following case:

$$\int_0^{100}\int_0^{100}\int_0^{0.5}\int_0^{0.5}\frac{378e^{-\frac{(x-y)^4}{10000}}e^{-625(x_3-y_3)^4}x^4}{\pi[0.02+\sqrt{((x-y)^2+(x_3-y_3)^2)^3}]}\mathrm dx\mathrm dy\mathrm dx_3\mathrm dy_3$$

When I evaluate my NIntegrate expression I get a message saying the integration is slow at converging. It also says it suspects that the cause is either a singularity or the value of integration being zero or highly oscillatory. I really need a way to integrate this monster, Any help?

In order to obtain this formula I wrote:

expr =
E^(-((x - y)^4/10000) - 625 (x3 -
y3)^4) ((e x + f x^2 + (m x + n x^2) x3 - e y -
f y^2 - (m y + n y^2) y3) ((
945 (x - y)^2 (e x + f x^2 + (m x + n x^2) x3 - e y -
f y^2 - (m y + n y^2) y3))/(
2 π (Sqrt[((x - y)^2 + (x3 - y3)^2)^3]) + 1/50) + (
945 (x - y) (x3 - y3) (g x + h x^2 + (o x + p x^2) x3 - g y -
h y^2 - (o y + p y^2) y3))/(
2 π (Sqrt[((x - y)^2 + (x3 - y3)^2)^3]) + 1/50)) + (g x +
h x^2 + (o x + p x^2) x3 - g y - h y^2 - (o y + p y^2) y3) ((
945 (x - y) (x3 - y3) (e x + f x^2 + (m x + n x^2) x3 - e y -
f y^2 - (m y + n y^2) y3))/(
2 π (Sqrt[((x - y)^2 + (x3 - y3)^2)^3]) + 1/50) + (
945 (x3 - y3)^2 (g x + h x^2 + (o x + p x^2) x3 - g y -
h y^2 - (o y + p y^2) y3))/(
2 π (Sqrt[((x - y)^2 + (x3 - y3)^2)^3]) + 1/50)));
kist = List @@ Expand[expr]


Then in order to evaluate the total integral I need to evaluate the integral of each term of the list. I am stuck at the first term which is the one I posted.

NIntegrate[
(945  E^(-((x - y)^4/10000) - 625 (x3 - y3)^4) x^4) /
(1/50 + 2 π Sqrt[((x - y)^2 + (x3 - y3)^2)^3]),
{x, 0, 100}, {y, 0, 100}, {x3, 0, 0.5}, {y3, 0, 0.5}]

• To save other people the trouble of translating this, as you say, "monster", into Mathematica, can you supply the NIntegrate[] code you used? – J. M. is away Mar 13 '18 at 15:46
• yes of course sorry – Riccardo Mar 13 '18 at 15:50
• Please edit your question to put in the code. – J. M. is away Mar 13 '18 at 15:50
• The integrand is strongly peaked at $x=y$, and essentially zero elsewhere. NIntegrate is likely the incorrect tool for this. Instead, obtain an asymptotic approximation using Laplace's method to get a very accurate (analytic) result. – QuantumDot Mar 13 '18 at 16:13
• I am not a mathematician so I might be hugely mistaken about this, but isn't the error for Laplace Method bounded for large exponents of "e"? If this were true, then I can apply the method to the x3,y3 variables (where I see as apex 625(x3-y3)^4)and not to the x,y, where I read as apex 10^-4 (x-y)^4; Did I get this right? – Riccardo Mar 13 '18 at 17:23

Given the partial symmetry and peaks along y == x and y3 == x3, it makes sense to sum the reflections of the integrand along these lines (or hyperplanes) and integrate the even part over half the x-y plane and half the x3-y3 plane:

# + Replace[#, {y -> x, x -> y}] &[expr];
# + Replace[#, {y3 -> x3, x3 -> y3}] &[%];
ca = Factor@ CoefficientArrays[%,
vars = Complement[Variables[%], {x, y, x3, y3}]];
coeffs = Flatten@Normal@ca // DeleteDuplicates // DeleteCases;
nint = 0;
Check[
++nint;
res = # -> Quiet[
NIntegrate[#, {x, 0, 100}, {y, x, 100}, {x3, 0, 1/2}, {y3, x3, 1/2},
PrecisionGoal -> 4, AccuracyGoal -> 10],
{NIntegrate::slwcon}],
Print[{nint, res}];
res
] & /@ coeffs];

Replace[Normal@ca[], subs, 2].vars.vars Note there are only quadratic terms:

ca Well, the coefficients are mostly smaller than @Bill's by orders of magnitude. There are no errors. (There are slow convergence warnings for most coefficients, but the numerical approximation converges nonetheless.) Maybe it's right.

• Regarding the result obtained, I see it's still a function of $x$ and $y$, but I don't understand why – Riccardo Mar 15 '18 at 11:22
• @Riccardo Bad substitution: replaced a factor, instead of the whole. I needed Replace instead of ReplaceAll (/.) Thanks for pointing it out. Fixed now. – Michael E2 Mar 15 '18 at 11:45

If want a fast answer with an error of about 1/1000, take into account, that this function peaks near x==y and x3==y3. Therefore limit the integration area.

Peaking is about Abs[x - y] < 20 and Abs[x3 - y3] < .30

Manipulate[
Plot3D[ (945 E^(-((x - y)^4/10000) - 625 (x3 - y3)^4) x^4)/(1/50 +
2 \[Pi] Sqrt[((x - y)^2 + (x3 - y3)^2)^3]), {x, 0, 100}, {y, 0,
100}, RegionFunction -> Function[{x, y, z}, Abs[x - y] < 20]], {x3,
0, 0.5}, {y3, 0, 0.5}] Manipulate[
Plot3D[ (945 E^(-((x - y)^4/10000) - 625 (x3 - y3)^4) x^4)/(1/50 +
2 \[Pi] Sqrt[((x - y)^2 + (x3 - y3)^2)^3]), {x3, 0, 0.5}, {y3, 0,
0.5}, RegionFunction ->
Function[{x3, y3, z3}, Abs[x3 - y3] < .30]], {{x, 50}, 0,
100}, {{y, 50}, 0, 100}] Integrate only this regions. Ignore error messages.

(nn = NIntegrate[
Boole[Abs[x - y] < 20] Boole[
Abs[x3 -
y3] < .30] (945 E^(-((x - y)^4/10000) -
625 (x3 - y3)^4) x^4)/(1/50 +
2 \[Pi] Sqrt[((x - y)^2 + (x3 - y3)^2)^3]), {x, 0, 100}, {y,
0, 100}, {x3, 0, 0.5}, {y3, 0, .5}]) // Timing

(*   {3.828, 3.79085*10^12}   *)