# [✓] Compute a triple integral with NIntegrate?

Posted 1 year ago
1520 Views
|
4 Replies
|
7 Total Likes
|
 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?
4 Replies
Sort By:
Posted 1 year ago
 http://reference.wolfram.com/language/ref/NIntegrate.htmltake a look at Options - Methods - Integration Strategies
Posted 1 year ago
 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 
 Your integrand has a removable singularity at the origin, which causes the numerical problems. You can reduce its severity by using the undocumented InternalExpm1[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_] :> InternalExpm1[y] + 1 // Simplify (* -((8 InternalExpm1[-(z^2/2)])/( v (u + v) z^2 (v + 2 z) (u + v + 2 z))) + (8 (InternalExpm1[-(1/2) (v + z)^2]/(v (v + z)^2 (v + 2 z)) - InternalExpm1[-(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} *) 
 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}.