Mathematica multi-dimensional numerical integration default method

Question

I'm performing multidimensional Numerical integrations with mathematica I was wondering what was the default method that mathematica was using. Also i'm changing some parameter inside the integration, and for some of those (which are not known before) the integrand is very small $(\approx 0)$ and the method becomes very slow; is there an automatic way to spot this?

here is an example of the code i'm using:

Gslip33timesxintegrated[h_, λ_, θ_, l_] := 
 NIntegrate[
 x*distance[s, h, θ]*
 Exp[-u/2]*(distance[s, h, θ] + λ*u)*
 GDz[surfacecoord[x, 0, θ, ν, h][[1]] - 
    imagelinedensitycoord[s, θ, h][[1]], 
   surfacecoord[x, 0, θ, ν, h][[2]] - 
    imagelinedensitycoord[s, θ, h][[2]], 
   surfacecoord[x, 0, θ, ν, h][[3]] - 
    imagelinedensitycoord[s, θ, h][[3]] + λ*u][[
  3]] -x*distance[s, h, θ]*Exp[-u/2]*
 GSDzz[surfacecoord[x, 0, θ, ν, h][[1]] - 
    imagelinedensitycoord[s, θ, h][[1]], 
   surfacecoord[x, 0, θ, ν, h][[2]] - 
    imagelinedensitycoord[s, θ, h][[2]], 
   surfacecoord[x, 0, θ, ν, h][[3]] - 
    imagelinedensitycoord[s, θ, h][[3]] + λ*u][[
  3]], {x, -l, l}, {s, -l, l}, {u, 0, Infinity}];

Where the function GDz and GSDzz are defined in the following way:

GDz[x_, y_, z_] = -D[sourcesing[x, y, z], z];
GSDzz[x_, y_, z_] = -D[stokesletz[x, y, z], z];

In turn the functions sourcesing and stokesletz are defined as:

sourcesing[x_, y_, z_] =
 {(x/(Sqrt[x^2 + y^2 + z^2])^(3)), (y/(Sqrt[
    x^2 + y^2 + z^2])^(3)), (z/(Sqrt[x^2 + y^2 + z^2])^(3))};

stokesletz[x_, y_, z_] = 
  {{1/(Sqrt[
     x^2 + y^2 + z^2]) + (x^2)/((Sqrt[x^2 + y^2 + z^2])^3), (x*
    y)/((Sqrt[x^2 + y^2 + z^2])^3), (x*
    z)/((Sqrt[x^2 + y^2 + z^2])^3)}, {(x*
    y)/((Sqrt[x^2 + y^2 + z^2])^3), 
   1/(Sqrt[x^2 + y^2 + 
      z^2]) + (y^2)/((Sqrt[x^2 + y^2 + z^2])^3), (y*
    z)/((Sqrt[x^2 + y^2 + z^2])^3)}, {(x*
    z)/((Sqrt[x^2 + y^2 + z^2])^3), (y*
    z)/((Sqrt[x^2 + y^2 + z^2])^3), 
 1/(Sqrt[x^2 + y^2 + 
      z^2]) + (z^2)/((Sqrt[x^2 + y^2 + z^2])^3)}}.{0, 0, 1};

the argument of the functions GDz and GSDzz is the difference between the components of the following two vectors:

surfacecoord[x_, r_, θ_, ν_, h_] = 
 {x*Cos[θ] + r*Sin[θ] Sin[ν], r*Cos[ν], 
  x*Sin[θ] + h - r*Cos[θ]*Sin[ν]};

imagelinedensitycoord[s_, θ_, h_] = {s*Cos[θ], 0, 
   s*Sin[θ] - h};

Finally the function distance is defined as:

distance[s_, h_, θ_] = h - s*Sin[θ];

what typically happens is that when $\theta$ is close to zero the integrand is nearly zero...

Welcome! Supplying more detail and code should help others answering your question. As it is we can only speculate. — Yves Klett, Oct 22 '14 at 19:10
Indeed i can supply the code, but the functions i have to integrate are very long and complicated... — SSC Napoli, Oct 22 '14 at 19:14
@YvesKlett just edited my question to include the code, here the functions however are implicit... — SSC Napoli, Oct 22 '14 at 19:28
sure! i just thought that there could have been a general answer — SSC Napoli, Oct 22 '14 at 19:46
What emerged from the question that you linked to really scares me, i'm doing an integration over a parallelepipedon domain so the integration order should not matter but the guy in the question got different results by changing the integration order...do mathematica at least returns some warning when this happens? — SSC Napoli, Oct 22 '14 at 20:48

Acus · Accepted Answer · 2014-10-23T12:02:08.313

Probably You have to play with AccyracyGoal and PrecisionGoal. However in this case I usually use universal trick by shifting integration function value from zero and then subtracting corresponding integral.

In you case, for example, the following works:

Gslip33timesxintegratedSubtract =
NIntegrate[ 1/(u + 1)^2, {x, -l, l}, {s, -l, l}, {u, 0, Infinity}];




Gslip33timesxintegrated2[h_, λ_, θ_, l_] := 
  NIntegrate[ 1/(u + 1)^2 +
    x*distance[s, h, θ]*
     Exp[-u/2]*(distance[s, h, θ] + λ*u)*
     GDz[surfacecoord[x, 0, θ, ν, h][[1]] - 
        imagelinedensitycoord[s, θ, h][[1]], 
       surfacecoord[x, 0, θ, ν, h][[2]] - 
        imagelinedensitycoord[s, θ, h][[2]], 
       surfacecoord[x, 0, θ, ν, h][[3]] - 
        imagelinedensitycoord[s, θ, h][[3]] + λ*
         u][[3]] - 
    x*distance[s, h, θ]*Exp[-u/2]*
     GSDzz[surfacecoord[x, 0, θ, ν, h][[1]] - 
        imagelinedensitycoord[s, θ, h][[1]], 
       surfacecoord[x, 0, θ, ν, h][[2]] - 
        imagelinedensitycoord[s, θ, h][[2]], 
       surfacecoord[x, 0, θ, ν, h][[3]] - 
        imagelinedensitycoord[s, θ, h][[3]] + λ*
         u][[3]], {x, -l, l}, {s, -l, l}, {u, 0, Infinity}];

Then the fast integration answer would be, for example,

Gslip33timesxintegrated2[1., 0.4, 0.001, 1.] - 
 Gslip33timesxintegratedSubtract

It is a trick and use it at your own risk.

Nice trick, I'll check how faster the code becomes and check if it's worth the risk and some extra verification! Maybe the optimization can be done in terms of which function to add/subtract is the fastest to integrate... — SSC Napoli, Oct 24 '14 at 15:53

Mathematica multi-dimensional numerical integration default method

1 Answers1