Numerical integration with Dirac delta

Question

I have some complicated function depending on many arguments $x,y,z$ and parameter $a$ multiplied by Dirac delta of another function,

$$ \tag 1 f(a,x,y,z) = g(a,x,y,z)\delta(t(a,x,y,z)) $$

I want to perform numerical integration over all the variables. Then, independently on the parameter a which has been chosen, the result is 0. However, evaluating the same integral analytically, I obtain non-zero result.

What is a reason that Mathematica doesn't evaluate the numerical integral correctly, and how to force it to do this? It seems that I can't use limit representations of the Dirac delta because of the numerical integration.

Edit. My example is

f[a_,x_,y_,z_] =(a^4 + 6400 a^2 - 81920000) Exp[-0.01 x^2] Sqrt[y^2 - a^2]Sin[z] UnitStep[11 Pi/45 - z] UnitStep[z - 11 Pi/90] DiracDelta[a^2 + 2 Sqrt[y^2 - a^2] x Cos[z] - 2 y Sqrt[x^2 + 6400] + 6400]

I want to integrate it over the region $x \in (0,3000), \ y\in (a,3000), \ z \in (0,\pi)$. I can perform the first step - evaluate one of the integrals by rewriting the Dirac delta as, say, $$ \delta (t(a,x,y,z)) = \frac{\delta(x-x_{1})}{|t'(x_{1})|}+\frac{\delta(x-x_{2})}{|t'(x_{2})|}, $$ where $x_{1,2} \equiv x_{1,2}(a,y,z)$ are solutions of $t(a,x,y,z) = 0$, and then to introduce reduced function depending only on $a,y,z$, $$ \tag 2 f_{1}(a,y,z) = \frac{g(a,x_{1},y,z)}{|t'(x_{1})|}+\frac{g(a,x_{2},y,z)}{|t'(x_{2})|}, $$ where $g$ is defined in $(1)$. However, this step introduces extra work which I would like to avoid.

Simple numerical integration

NIntegrate[f[2, x, y, z], {x, 0, 3000}, {y, 2, 3000}, {z, 0, Pi}]

gives zero, but the integration performed with $(2)$ is non-zero.

Answer 1

You can integrate symbolically of z and then use the result in NIntegrate.

newInt = Integrate[f[2, x, y, z], {z, 0, \[Pi]}]

(* -(1/x)
 40947192 E^(-(x^2/
   100)) (HeavisideTheta[
     3202 - Sqrt[6400 + x^2] y + 
      x Sqrt[-4 + y^2] Cos[(11 \[Pi])/90]] - 
    HeavisideTheta[
     3202 - Sqrt[6400 + x^2] y + x Sqrt[-4 + y^2] Cos[(11 \[Pi])/45]]) *)

NIntegrate[newInt, {x, 0, 3000}, {y, 2, 3000}, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4]

(* -4.85146*10^6 *)

NIntegrate[newInt, {x, 0, 3000}, {y, 2, 3000}, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4, 
 WorkingPrecision -> 30]

(* -4.19942269494560346742215774341*10^7 *)

Take note of NIntegrate's messages. Is this producing results you expect?

Answer 2

Here's an extension of Anton's approach:

newInt = Integrate[f[2, x, y, z], {z, 0, π}] /. 
  HeavisideTheta -> UnitStep; (* replacement by UnitStep is optional *)

The integrand has the form k * Exp[_] * (A-B) / x, where A and B are UnitStep (or HeavisideTheta) functions.

Block[{A, B},  (* there are two cases of UnitStep/HeavisideTheta *)
  {A, B} = Cases[newInt, UnitStep[u_] :> u, Infinity];
  reg = Reduce[  (* find reg where UnitStep[A]-UnitStep[B] is nonzero *)
    (A > 0 && B < 0 || A < 0 && B > 0) &&  
     0 < x < 3000 && 2 < y < 3000, {x, y}]
  ];

Cases[DeleteCases[reg, _Equal && _], (* delete zero-measure subsets *)
 _[a_, ___, x, ___, b_] && _[c_, ___, y, ___, d_] :> 
  NIntegrate[newInt, {x, a, b}, {y, c, d}]]
Total[%]

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions....

(*
  {-4.23643*10^7, 0.}
  -4.23643*10^7
*)

On the second part of reg, the exponential factor of newInt underflows.