2
$\begingroup$

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}]
$\endgroup$
  • 2
    $\begingroup$ To save other people the trouble of translating this, as you say, "monster", into Mathematica, can you supply the NIntegrate[] code you used? $\endgroup$ – J. M. is away Mar 13 '18 at 15:46
  • $\begingroup$ yes of course sorry $\endgroup$ – Riccardo Mar 13 '18 at 15:50
  • $\begingroup$ Please edit your question to put in the code. $\endgroup$ – J. M. is away Mar 13 '18 at 15:50
  • $\begingroup$ 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. $\endgroup$ – QuantumDot Mar 13 '18 at 16:13
  • $\begingroup$ 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? $\endgroup$ – Riccardo Mar 13 '18 at 17:23
1
$\begingroup$

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

Mathematica graphics

Note there are only quadratic terms:

ca

Mathematica graphics

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.

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

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

enter image description here

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

enter image description here

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}   *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.