Numerical integration with Dirac delta

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.

• The problem with DiracDelta is that it will evaluate to 0 when fed by a nonzero numerical argument. This way, $\delta(t(a,x,y,z))$ will quite likely be interpreted as function that is zero almost everywhere (but in reality, $\delta$ is not a function). It is better not to use it for numerical code. If I understand you correctly, you want to intergrate over the hypersurface defined by the equation $t(a,x,y,z) = 0$. Please, provide the concrete equations. – Henrik Schumacher Sep 30 '18 at 10:08
• In some applications it shoudl also be noted, that KroneckerDelta and DiscreteDelta and DiracDelta are not the same thing. – Johu Sep 30 '18 at 11:32
• Maybe something like NIntegrate[ g[a, x, y, z], {a, x, y, z} \[Element] ImplicitRegion[t[a, x, y, z] == 0, {a, x, y, z}]] – Szabolcs Sep 30 '18 at 11:56
• @HenrikSchumacher : I've added my example. Please take a look on it. – John Taylor Sep 30 '18 at 12:33
• @Szabolcs : applied to my example (see the update of my question please) it doesn't give the result (only displays the code in the output). – John Taylor Sep 30 '18 at 12:38

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?

• See my comment to the question. – user64494 Sep 30 '18 at 15:52

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.

• See my comment to the question. – user64494 Sep 30 '18 at 15:51
• @user64494 Why? how is your comment in any way relevant? either to the OP or this answer? It doesn't add any information related to either the question or the answer, and it is completely unnecessary. – AccidentalFourierTransform Sep 30 '18 at 17:19
• @AccidentalFourierTransform: I am sure that my comment adds useful information: in the Properties in n dimensions section of the cited Wiki article such type integrals are treated – user64494 Sep 30 '18 at 19:46