# How to explain these results of integration of DiracDelta?

Let us consider

ClearAll["Global*"];
Integrate[DiracDelta[1 - x^2 - y^2], {x, -1, 1}, {y, -1, 1}]
(*\[Pi]/2*)


Let us verify it, making use of approximations of delta-distribution in weak topology (more exactly, making use of usual functions associated with these approximations) and usual double integrals.

eps = 0.005; NIntegrate[ eps/((1 - x^2 - y^2)^2 + eps^2)/Pi, {x, -1, 1}, {y, -1, 1}]
(*2.947*)
NIntegrate[1/Sqrt[Pi]/eps*Exp[-(1 - x^2 - y^2)^2/eps^2], {x, -1, 1}, {y, -1, 1}]
(*0.468889*)


I have never seen double integrals of distributions over bounded sets in math literature and don't know any definitions of such integrals. I would be very thankful for accessible and serious references. How to explain the difference between the exact calculation (BTW, Maple produces $$\pi$$) and the numeric results?

• Interestingly Integrate[DiracDelta[1 - Sqrt[x^2 + y^2]], {x, -1, 1}, {y, -1, 1}]  yields π – chris Dec 28 '20 at 9:36
• The general idea of going from infinite ranges to finite has been explained a few times I think. In this example more is required since the integrand has four singular points on the integration boundary. A limit argument allows to work around that detail. – Daniel Lichtblau Dec 28 '20 at 14:01
• (1) Four singular points on the integration boundary are x=+-1 and y=+-1. What are the other four? – Daniel Lichtblau Dec 29 '20 at 1:42
• (2) It is not my job to fill gaps in your knowledge. If you want a hint for how to extend to finite bounds, ask nicely. Making annoying remarks like "I prefer formulas... over nice, but ungrounded, words" is simply impolite. And, worse, it seems to be your modus operandi. Bad habit. Really bad habit. Try to lose it before people just start to ignore you. Me, I am on that verge. – Daniel Lichtblau Dec 29 '20 at 1:45
• (3) To state the utterly obvious, Mathematica might be giving an erroneous result. – Daniel Lichtblau Dec 29 '20 at 1:46

A very simple limiting representation of the Dirac $$\delta$$-function is

f[ε_, x_] = Piecewise[{{1/(2*ε), -ε < x < ε}}];


with $$\lim_{\epsilon\to0^+} f_{\epsilon}(x)=\delta(x)$$ in the sense of a distribution.

The integral in question is then

Assuming[0 < ε < 1/2,
Integrate[f[ε, 1-x^2-y^2], {x, -1, 1}, {y, -1, 1}] // FullSimplify]

(*    (2 Sqrt[ε] + (1+ε)*ArcCot[Sqrt[ε]] + (ε-1)*ArcCsc[Sqrt[1/ε-1]] +
(ε-1)*ArcTan[Sqrt[1/ε-2]] - ArcTan[Sqrt[ε]] - ε*ArcTan[Sqrt[ε]])/ε    *)


The limit is $$\pi$$:

Limit[%, ε -> 0, Direction -> "FromAbove"]
(*    π    *)

• +1. Thank you for your approach. Do you mean "in the sense of distributions"? – user64494 Dec 28 '20 at 15:35
• @user64494 isn't that precisely what I wrote? – Roman Dec 28 '20 at 15:41
• The questions are still open: why the difference with the result ofIntegrate[DiracDelta[1 - x^2 - y^2], {x, -1, 1}, {y, -1, 1}]? why the difference with the result of eps = 0.005; NIntegrate[1/Sqrt[Pi]/eps*Exp[-(1 - x^2 - y^2)^2/eps^2], {x, -1, 1}, {y, -1, 1}]? – user64494 Dec 28 '20 at 15:45
• The difference with your Integrate looks like a bug. The difference with your NIntegrate is, as @UlrichNeumann says, a numerical artefact: try f[eps_?NumericQ] := NIntegrate[1/Sqrt[Pi]/eps*Exp[-(1 - x^2 - y^2)^2/eps^2], {x, -1, 1}, {y, -1, 1}] and then Plot[f[eps], {eps, 0, 1}] and you'll see that it converges to $\pi$ as $\text{eps}\to0^+$ apart from numerical errors for small eps. – Roman Dec 28 '20 at 16:17
• @user64494 Haha of course it's not a proof. There's a whole world of mathematical intuition that lies outside of proofs; numerical analysis is a crazy mix of solid math and shaky heuristics. For proofs I think you should go to the math stackexchange. – Roman Dec 28 '20 at 18:09

Here a way to get the correct symbolic result in Mathematica

Integrate[DiracDelta[1 - x^2 - y^2], {x, -1 - delta,1 + delta}, {y, -1 - delta, 1 + delta}
, Assumptions -> delta > 0]
(*Pi*)


numerical verification

eps = 0.005;
delta = .1
NIntegrate[1/Sqrt[2 Pi eps] Exp[-((1 - x^2 - y^2)^2/(2 eps))], {x, -1 - delta,1 + delta}, {y, -1 - delta, 1 + delta} ]
(*3.14091*)

• Why do you think Pi is the correct symbolic result? Can you kindly ground that claim? Concerning the numeric verification. First, you changed the notations. Second, putting your eps=0.0001 in your code, one obtains 1.86765. – user64494 Dec 28 '20 at 12:25
• Cut and paste error, I modified my answer. My claim is grounded , see my comment concerning maple result. By the way you shouldn't ignore the NIntegrate messages concerning convergence – Ulrich Neumann Dec 28 '20 at 12:31
• Warning about a slow convergence is not any error. In any case you don't explain the defference between integration over {x, -1 - delta,1 + delta}, {y, -1 - delta, 1 + delta} and {x, -1, 1}, {y, -1, 1}. The heart of the problem is not explained and the definition of the double integral of DiracDelte[1-x^2-y^2] is no done. – user64494 Dec 28 '20 at 12:36
• Do you know why the fact that the contour touches the edge yields a factor of Pi/2 drop? Is this a known mathematical result or a kind of bug of the mathematica implementation? – chris Dec 28 '20 at 15:36
• @chris: Sorry, don't understand you. The question is that different approaches produce different results. I don't know why and how to explain it. – user64494 Dec 28 '20 at 16:17

The Mathematica expression

Integrate[
DiracDelta[(x^2 + y^2) - 1], {x, -1, 1}, {y, -1, 1}]


evaluates (wrongly) to Pi/2 since your integration limits coincide with (x, y) points such that the argument of DiracDelta is equal zero (correct value: Pi, since we integrate over unit circle).

Thus, increasing the integration limit solves the problem:

Integrate[
DiracDelta[(x^2 + y^2) - 1], {x, -2, 2}, {y, -2, 2}]


evaluates to Pi.

On the other hand, also integration over the upper half-plane yields the correct result:

Integrate[
DiracDelta[(x^2 + y^2) - 1], {x, -1, 1}, {y, 0, 1}]


evaluates correctly to Pi/2.

Furthermore, by changing from cartesian to polar coordinates, Mathematica gives you a hint on why it's acting strange:

Integrate[r*DiracDelta[r^2 - 1], r, {\[CurlyPhi], 0, 2 \[Pi]}]


evaluates to Pi HeavisideTheta[r^2 - 1]. Note that HeavisideTheta != UnitStep, i.e. undefined behavior for r=1.

As a conclusion, your integration limit may not be element of your region of integration, which I think is a wrong behavior of mathematica.

• Your statement "evaluates (wrongly) to Pi/2 since your integration limits coincide with (x, y) points such that the argument of DiracDelta is equal zero" is built on the sand. An integration delta = 1/100; Integrate[ DiracDelta[1 - x^2 - y^2], {x, -1 + delta, 1 - delta}, {y, -1 + delta, 1 - delta}] which results in -\[Pi] + 4 ArcTan[99/Sqrt[199]] contradicts your cllaim since there are 8 points on the boundary of the square $[-1+\delta,1-\delta]^2$ where $1-x^2-y^2=0$ . Two of them are $(1-\delta,\pm\sqrt{1-(1-\delta^2})$. I leave to find other points on your own. – user64494 Dec 29 '20 at 16:25
• You may find the "integral" Integrate[ DiracDelta[1 - x^2 - y^2], {x, -1 + delta, 1 - delta}, {y, -1 + delta, 1 - delta}, Assumptions->delta>0&&delta<1/4] in closed form and verify that its limit as delta tends to zero from above equals Pi. Good luck! Don't hesitate to ask for further explanation in need. – user64494 Dec 29 '20 at 17:26
• This HeavisideTheta appearance might well be the issue. Pretty much as @chris noted in a comment to a different response. (If I told my colleagues their debugging was built on sand, they'd save enough of it to bury me. While the poster might not appreciate this effort to locate the problem source, I certainly do.) – Daniel Lichtblau Dec 30 '20 at 2:46

I will focus on only one issue raised in the OP, namely, why MA produces the $$\pi/2$$ result. I avoid on purpose any discussion on the mathematical justification of the integrals as I think this

I have never seen double integrals of distributions over bounded sets in math literature and don't know any definitions of such integrals. I would be very thankful for accessible and serious references.

goes beyond the scope of this forum. Coming back to this integral:

Integrate[DiracDelta[1 - x^2 - y^2], {x, -1, 1}, {y, -1, 1}]


1. Mathematica treats it as an iterated integral. Thus, it is sufficient to consider only the inner integral

Integrate[DiracDelta[1 - x^2 - y^2], {x, -1, 1}]


(* ConditionalExpression[Boole[-1 < -Sqrt[1 - y^2] < 1]/(2 Sqrt[Abs[1 - y^2]]), -1 < Sqrt[1 - y^2] < 1] *)


2. The correct answer can obtained by the integration in the limits $$[-a,a]$$

 z=Integrate[DiracDelta[1 - x^2 - y^2], {x, -a, a}]
(*(Boole[-a < -Sqrt[1 - y^2] < a || a < -Sqrt[1 - y^2] < -a] + Boole[-a < Sqrt[1 - y^2] < a || a < Sqrt[1 - y^2] < -a])/(2 Sqrt[Abs[1 - y^2]])*)


and setting $$a=1$$ (no limit is necessary)

z1=z/.{a->1};
Integrate[z1, {y, -1, 1}]
(*    π    *)


3. Alternatively we can FullSimplify the intermediate result

za=FullSimplify[z, Assumptions -> a > 0 && -1 <= y <= 1]
(* Boole[Sqrt[1 - y^2] < a]/Sqrt[1 - y^2] *)


Integration over y can be performed in full generality

Integrate[Boole[Sqrt[1 - y^2] < a]/Sqrt[1 - y^2], {y, -1, 1}]


$$\begin{array}{cc} \Bigg\{ & \begin{array}{cc} \pi & a>1 \\ 2 \sin ^{-1}(a) & 0

Conclusion The reason for the $$\pi/2$$ answer seems to be that mathematica forgot one root of the

$$1-x^2-y^2=0 \tag{2}$$

equation. It certainly knows how to perform integrals of the $$\delta$$-function of a function

$$\int\delta(f(x))\,dx=\int \sum_k\frac{\delta(x-\alpha_k)}{\left|f'(\alpha_k)\right|}\,dx\tag{3}$$

where $$\alpha_k$$ are the real roots of $$f(x)=0$$, however, in the OP example $$a=1$$ it misses one root and does not miss it otherwise.

It has been suggested by other answers that because the circle touches the boundaries of the integration domain this should somehow influence the result. It should not because these are points of measure zero: integration domain is $$2D$$, $$\delta$$-function cuts out a $$1D$$ manifold out of it, whereas ill-defined integrand is on the $$0D$$ manifold. This seems to be just a coincidence that for $$a=1$$ mathematica misses one root of (3).

• First of all, thank you for your work and answer. +1. Indeed, Mathematica treats the integral under consideration as an iterated integral. The question arises: how MMA produces Integrate[DiracDelta[1 - x^2 - y^2], {x, -a, a}]==(Boole[-a < -Sqrt[1 - y^2] < a || a < -Sqrt[1 - y^2] < -a] + Boole[-a < Sqrt[1 - y^2] < a || a < Sqrt[1 - y^2] < -a])/(2 Sqrt[Abs[1 - y^2]])? In a usual integral the Newton-Leibniz formula is applied or the limit of an integral sum is found. What does Mathematica apply here? I want to know the truth, don't cavil about trifles. – user64494 Dec 31 '20 at 9:28
• @user64494 Thank you. Good question, I do not know actually. Only developers can say how they programmed Integrate`. I would say that eq. (3) is converted to a sum over roots by adding a Heaviside function to assure that these roots fall into the integration interval. If you look into the Rubi package, there is huge number of rules for all possible cases of life (it cannot deal with deltas). I would guess that MA too first reduces the integral to known forms. Newton-Leibniz is probably the last rule that is being tried for definite integrals. – yarchik Dec 31 '20 at 10:09