[?] Compute a triple integral with NIntegrate?

Questions
1

I am trying to numerically compute a triple integral. The code is as follows:

I61 =(8 (((-1+E^(-(1/2) (u+v+z)^2)) z^2)/(u+v+z)^2-((-1+E^(-(1/2) (u+v+z)^2)) (v+z)^2)/(u+v+z)^2+((-E^(-(z^2/2))+E^(-(1/2) (u+v+z)^2)) (v+z)^2)/((u+v) (u+v+2 z))+(E^(-(1/2) (u+v+z)^2) (-1+E^(1/2 u (u+2 (v+z)))) z^2)/(u (u+2 (v+z)))))/(v z^2 (v+z)^2 (v+2 z))
I62 = Simplify[Together[I61]]
NIntegrate[I62, {z, 0, Infinity}, {u, 0, Infinity}, {v, 0, Infinity}, PrecisionGoal -> 4]

This works if I set PrecisionGoal to be 4. But increasing it 5 causes an error:

"The integrand ... has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the 
region with boundaries {{3.,1.}, {0.,1.70984*10^14},{\[Infinity],7.}}"

Is there anyway to increase the precision a bit?

2

take a look at Options - Methods - Integration Strategies

3

I would suggest considering where the integrand evaluates to Indeterminate and why; perhaps the contributions from such regions might be negligible anyway. Ignoring these gives something like

In[4]:= f[z_?NumericQ, u_, v_] := I62 /. Indeterminate -> 0

In[5]:= N[NIntegrate[f[z, u, v], {z, 0, Infinity}, {u, 0, Infinity}, {v, 0, Infinity}, 
            PrecisionGoal -> 7, MinRecursion -> 3, WorkingPrecision -> 100], 7]

Out[5]= 1.374413
4

Your integrand has a removable singularity at the origin, which causes the numerical problems. You can reduce its severity by using the undocumented Internal`Expm1[x], which is equivalent to Exp[x] - 1:

I63 = Collect[I62 // Expand, Power[E, _], Simplify] /. Power[E, y_] :> Power[E, Simplify[y]]
(*
8/(z^2 (v + z)^2 (u + v + z)^2) -
 (8 E^(-(z^2/2))) / (v (u + v) z^2 (v + 2 z) (u + v + 2 z)) +
 (8 E^(-(1/2) (v + z)^2)) / (u v (v + z)^2 (v + 2 z) (u + 2 (v + z))) -
 (8 E^(-(1/2) (u + v + z)^2)) / (u (u + v) (u + v + z)^2 (u + v + 2 z) (u + 2 (v + z)))
*)

I64 = I63 /. Power[E, y_] :> Internal`Expm1[y] + 1 // Simplify
(*
-((8 Internal`Expm1[-(z^2/2)])/(
  v (u + v) z^2 (v + 2 z) (u + v + 2 z))) +
 (8 (Internal`Expm1[-(1/2) (v + z)^2]/(v (v + z)^2 (v + 2 z)) - 
      Internal`Expm1[-(1/2) (u + v + z)^2]/((u + v) (u + v + z)^2 (u +  v + 2 z)))) /
  (u (u + 2 (v + z)))
*)

MemoryConstrained[
  NIntegrate[I64, {z, 0, Infinity}, {u, 0, Infinity}, {v, 0, Infinity},
   Method -> {"GlobalAdaptive",  "SingularityHandler" -> "DuffyCoordinates"},
   PrecisionGoal -> 8],
  6000000000] // AbsoluteTiming
(*
NIntegrate::slwcon: Numerical integration converging too slowly....

{4.45456, 1.3744127474237668`}
*)
5

Actually the singularity handler "DuffyCoordinates" is enough.

NIntegrate[I62, {z, 0, Infinity}, {u, 0, Infinity}, {v, 0, Infinity},
  Method -> {"GlobalAdaptive", "SingularityHandler" -> "DuffyCoordinates"}, 
  PrecisionGoal -> 6] // AbsoluteTiming
(*
  NIntegrate::slwcon: Numerical integration converging too slowly....

  {3.011657`, 1.3744127589904662`} 
*)

But to get greater precision, a higher working precision and more time are needed:

NIntegrate[I62, {z, 0, Infinity}, {u, 0, Infinity}, {v, 0, Infinity},
  Method -> {"GlobalAdaptive", "SingularityHandler" -> "DuffyCoordinates"},
  WorkingPrecision -> 24, PrecisionGoal -> 8] // AbsoluteTiming
(*
  NIntegrate::slwcon: Numerical integration converging too slowly....

  {70.5418, 1.37441274742384630555707}
*)

One can save a small amount of time by separating out the singularity at z == 0 with {z, 0, 1, Infinity} in place of {z, 0, Infinity}.