Trouble with numerical evaluation of a four-fold integral

Question

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}]

Answer 1

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[0];
nint = 0;
subs = Thread[
   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[[3]], 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.

Answer 2

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}   *)