# How to interpret these integrations with DiracDelta?

Let us consider the double integral of DiracDelta[y-x^2] over the triangle with vertices at $$(-1,0),\,(-1,2),\,(1,0)$$. Mathematica calculates it as an iterated integral

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

Let us trace that step by step.

Integrate[DiracDelta[y - x^2], {y, 0, 1 - x}] // InputForm
(*ConditionalExpression[DiscreteDelta[-1 + x] -  HeavisideTheta[1/2 - Sqrt[5]/2 + x] +
HeavisideTheta[(1 + Sqrt[5])/2 + x],Element[x, Reals]]*)
Integrate[%,{x,-1,1}]
(*1/2 (1 + Sqrt[5])*)

Now let us integrate in reverse order.

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

We see a different result from the previous one. Let us trace it.

Integrate[DiracDelta[y - x^2], {x, -1, 1 - y}] // InputForm
(*ConditionalExpression[Boole[-1 < -Sqrt[y] < 1 - y || 1 - y <  -Sqrt[y] < -1]/
(2*Sqrt[Abs[y]]),  Sqrt[y] < 1 || 1 + Sqrt[y] <  y]*)
Integrate[%,{y,0,2}]
(*Undefined*)

It's unclear whether the double integral exists in this case and how that double integral is defined. I have never seen definite multiple integrals of distributions over bounded regions in math literature. Constructive answers are welcome.

• Additionally, InverseFourierTransform[ FourierTransform[DiracDelta[y - x^2], {x, y}, {t, s}], {t, s}, {x, y}] produces 1/(\[Pi] (x^2 - y)) + DiracDelta[x^2 - y]. Commented Nov 8, 2020 at 14:36

The integral exists and the first result is ok.

Obviously Mathematica can't transform the nonlinear argument y-x^2 in your second approach Integrate[DiracDelta[y - x^2], {y, 0, 2}, {x, -1, 1 - y}]. The transformation by hand uses the identity DiracDelta[g[z]]=Sum[(z-zi)/Abs[g'[zi]]], g[zi]==0.

Integrate[DiracDelta[x + Sqrt[y] ]/( 2 Sqrt[y]) + DiracDelta[x - Sqrt[y] ]/(2 Sqrt[y]), {y, 0, 2}, {x, -1, 1 - y}]
(*1/2 (1 + Sqrt[5])*)
• Thank you. You calculate the second iterated integral in another way. However, your statement "The integral exists" is ungrounded. I repeat "How that double integral is defined?". I have never seen definite multiple integrals of distributions over bounded regions in math literature. Commented Nov 9, 2020 at 7:13
• Mathematica actually can't handle correctly the case DiracDelta[g[x,y]]. In your case the inner integration variable y is quadratic. My way is a well grounded workaround! Look for example for "Green's Function" to see multiple examples of integrals in bounde regions. Commented Nov 9, 2020 at 7:26
• @user64494 You might find this paper interesting. eej.aut.ac.ir/article_511_163e09592b261d64a85a3ca7673d9fbe.pdf Commented Nov 9, 2020 at 18:55
• @BillWats: Having looked in the article, I didn't see multiply integrals from distributions over bounded regions there. Commented Nov 9, 2020 at 23:26
• @user64494 The subclause in your comment is dispensable . Regardless , one link with multiple references: Weisstein mathworld Commented Nov 10, 2020 at 6:57

Let's take your first numerical integral

e = 10^-3; d = e/Pi*NIntegrate[1/(e^2 + (y - x^2)^2), {x, -1, 1}, {y, 0, 1 - x}, AccuracyGoal -> 3, PrecisionGoal -> 3]
(*1.57354*)

Change things up a little. Use Integrate instead of NIntegrate and make e smaller.

e = 10^-4; d = e/Pi*Integrate[1/(e^2 + (y - x^2)^2), {x, -1, 1}, {y, 0, 1 - x}] // N
(*1.60391 + 0. I*)

For e = 10^-5 we get d = 1.61356 +0. I and

for e = 10^-6 we get d = 1.6166 +0. I and

for e = 10^-10 we get 1.61801984663683206818852997777 + 0.*10^-30 I when I use N[#,30]& instead of N.

Clearly we are converging to Mathematica's exact answer.

• Have you looked my comment to my answer? How about g[u_]:=Exp[-u^2/e^2]/Sqrt[Pi]/e? Commented Nov 11, 2020 at 5:19

Too long for a comment.

Let us follow the primitive approach to the one-dimensional $$\delta$$-distribution and let us approximate DiracDelta under the consideration by a bell-shaped function and let us consider a usual two-dimensional integral. This can be done in many ways (compare Wiki and Encyclopedia of Mathematics). We start with a rational approximation f[u_}:=e/Pi/(e^2+u^2) in hopes to obtain a closed-form expression for the integral under consideration. Later we consider the case g[u_]:=Exp[-u^2/e^2]/Sqrt[Pi]/e too. Let's begin.

ClearAll[e]; a =  e/Pi*Integrate[1/(e^2 + (y - x^2)^2), {x, -1, 1}, {y, 0, 1 - x},
Assumptions -> e > 0 && e < 1]
(*1/(2 \[Pi]) e (1/e 4 (ArcCot[e] + (-1)^(1/4) Sqrt[
e] (-ArcTan[(-1)^(1/4)/Sqrt[e]] +
ArcTanh[(-1)^(1/4)/Sqrt[e]])) +
RootSum[1 + e^2 + 2 #1 - #1^2 - 2 #1^3 + #1^4 &, (
Log[1 - #1] #1)/(-1 - #1 + #1^2) &] -
RootSum[1 + e^2 + 2 #1 - #1^2 - 2 #1^3 + #1^4 &, (
Log[-#1] #1)/(-1 - #1 + #1^2) &] +
RootSum[1 + e^2 - 2 #1 - #1^2 + 2 #1^3 + #1^4 &, (
Log[1 - #1] #1)/(-1 + #1 + #1^2) &] -
RootSum[1 + e^2 - 2 #1 - #1^2 + 2 #1^3 + #1^4 &, (
Log[-#1] #1)/(-1 + #1 + #1^2) &])*)
N[a /. e -> 0.001, 20]
(*1.57354 + 0. I*)

I'd like to note N[1/2 + Sqrt[5]/2, 20] outputs 1.6180339887498948482. The above result is confirmed numerically by

e = 10^-3; d = e/Pi*NIntegrate[1/(e^2 + (y - x^2)^2), {x, -1, 1}, {y, 0, 1 - x},
AccuracyGoal -> 3, PrecisionGoal -> 3]
(*1.57353*)

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small

(A warning is not an error communication.) and, making use of the MonteCarlo method,

b = e/Pi*NIntegrate[1/(e^2 + (y - x^2)^2), {x, -1, 1}, {y, 0, 1 - x},
Method -> {"MonteCarlo", "MaxPoints" -> 10^7},AccuracyGoal -> 3,PrecisionGoal -> 2]
(*1.58846*)

We consider an integral of a continuous function over a tiangle which exists without shade of doubt and calculate it as an iterated integral. The consideration of the iterated integral in reverse order is not necessary, however, we treat that too. We split the iterated integral into three ones.

ClearAll[e]; e/Pi*Integrate[1/(e^2 + (y - x^2)^2), {y, 1, 2}, {x, -1, 1 - y},
Assumptions -> e > 0 && e < 1]
(*-((\[Pi] + I \[Sqrt](-1 - I e) \[Pi] - I \[Sqrt](-1 + I e) \[Pi] -
2 I \[Sqrt](-1 - I e) ArcTan[\[Sqrt](-1 - I e)] +
2 I \[Sqrt](-1 + I e) ArcTan[\[Sqrt](-1 + I e)] - (1 +
I) \[Sqrt](-10 I +
8 e) ArcTan[(1 + I)/(\[Sqrt](-10 I + 8 e))] + (1 +
I) \[Sqrt](10 I +
8 e) ArcTanh[(1 + I)/(\[Sqrt](10 I + 8 e))] -
I Log[-1 - I e] + I Log[-1 + I e])/(2 \[Pi]))*)
e/Pi*Integrate[1/(e^2 + (y - x^2)^2), {y, 0, 1}, {x, -1, 0},
Assumptions -> e > 0 && e < 1]
(*1/(2 \[Pi]) (\[Pi]+2 (-1)^(1/4) Sqrt[e] (ArcTanh[(-1)^(1/4)/Sqrt[e]]+
I ArcTanh[(-1)^(3/4)/Sqrt[e]])-2 (-1)^(3/4) Sqrt[-I+e] ArcTanh[(-1)^(3/4)/Sqrt[-I+e]]-
2 (-1)^(1/4) Sqrt[I+e] ArcTanh[(-1)^(1/4)/Sqrt[I+e]]+I (Log[I-e]-Log[I+e]))*)

Mathematica fails with the third piece

e/Pi*Integrate[1/(e^2 + (y - x^2)^2), {y, 0, 1}, {x, 0, 1 - y},
Assumptions -> e > 0 && e < 1]
(*(e Integrate[((-1)^(
1/4) (\[Sqrt](e +
I y) ArcTanh[((-1)^(3/4) (-1 + y))/(\[Sqrt](e - I y))] +
I \[Sqrt](e -
I y) ArcTanh[((-1)^(
1/4) (-1 + y))/(\[Sqrt](e + I y))]))/(2 e \[Sqrt](e^2 +
y^2)), {y, 0, 1}, Assumptions -> e > 0 && e < 1])/\[Pi]*)

It should be noticed, that the latest integral can be done numerically for concrete values of e. Let us continue.

-((\[Pi] + I \[Sqrt](-1 - I e) \[Pi] - I \[Sqrt](-1 + I e) \[Pi] -
2 I \[Sqrt](-1 - I e) ArcTan[\[Sqrt](-1 - I e)] +
2 I \[Sqrt](-1 + I e) ArcTan[\[Sqrt](-1 + I e)] - (1 +
I) \[Sqrt](-10 I +
8 e) ArcTan[(1 + I)/(\[Sqrt](-10 I + 8 e))] + (1 +
I) \[Sqrt](10 I +
8 e) ArcTanh[(1 + I)/(\[Sqrt](10 I + 8 e))] -
I Log[-1 - I e] + I Log[-1 + I e])/(2 \[Pi])) +
1/(2 \[Pi]) (\[Pi] +
2 (-1)^(1/4) Sqrt[
e] (ArcTanh[(-1)^(1/4)/Sqrt[e]] +
I ArcTanh[(-1)^(3/4)/Sqrt[e]]) -
2 (-1)^(3/4) Sqrt[-I + e] ArcTanh[(-1)^(3/4)/Sqrt[-I + e]] -
2 (-1)^(1/4) Sqrt[I + e] ArcTanh[(-1)^(1/4)/Sqrt[I + e]] +
I (Log[I - e] - Log[I + e])) /. e -> 0.001
(*0.977684 - 1.78111*10^-14 I*)

and

e = 0.001;(e *NIntegrate[((-1)^(
1/4) (\[Sqrt](e +
I y) ArcTanh[((-1)^(3/4) (-1 + y))/(\[Sqrt](e - I y))] +
I \[Sqrt](e -
I y) ArcTanh[((-1)^(
1/4) (-1 + y))/(\[Sqrt](e + I y))]))/(2 e \[Sqrt](e^2 +
y^2)), {y, 0, 1}])/\[Pi]
(*0.595856 - 1.37879*10^-17 I*)

This is in accordance with the results of completely numeric integration:

e=0.001;e/Pi*NIntegrate[1/(e^2 + (y - x^2)^2), {y, 0, 2}, {x, -1, 1 - y},
AccuracyGoal -> 3, PrecisionGoal -> 3]
(*1.57354*)

and

e = 10^-3; c = e/Pi*NIntegrate[1/(e^2 + (y - x^2)^2), {y, 0, 2}, {x, -1, 1 - y},
Method -> {"MonteCarlo", "MaxPoints" -> 10^7}, AccuracyGoal -> 3,PrecisionGoal -> 2]
(*1.56421*)

Now we turn to the case of g[u_]:=Exp[-u^2/e^2]/Sqrt[Pi]/e. Both iterated symbolic integrations fail, but both numeric integrations produce a different result from the previous ones.

e=0.001;e^-1/Sqrt[Pi]*NIntegrate[Exp[-(y - x^2)^2/e^2], {x, -1, 1}, {y, 0, 1 - x}]
(*0.10912*)

and

e=0.001;e^-1/Sqrt[Pi]*NIntegrate[Exp[-(y-x^2)^2/e^2],{y,0,2},{x,-1,1-y}]
(*0.111262*)

The executed codes as nb file on demand through Dropbox.

Addition. In this case one of the iterated integrals is reduced to an one-dimensional integral

ClearAll[e];e^-1/Sqrt[Pi]*Integrate[Exp[-(y - x^2)^2/e^2], {x, -1, 1}, {y, 0, 1 - x},
Assumptions -> e > 0 && e < 1]
(*(1/(e Sqrt[\[Pi]]))Integrate[1/2 e Sqrt[\[Pi]] (Erf[x^2/e] - Erf[(-1 + x + x^2)/e]),
{x, -1, 1},Assumptions -> e > 0 && e < 1]*)

Now calculate the latest integral numerically

e := 0.01;NIntegrate[1/2*e *Sqrt[\[Pi]]*(Erf[x^2/e]-Erf[(-1 + x + x^2)/e]), {x, -1, 1}]
(*0.0274534*)

e:=0.00100000000000000000000000000000000;
NIntegrate[1/2*e *Sqrt[\[Pi]] *(Erf[x^2/e] - Erf[(-1 + x + x^2)/e]), {x, -1, 1},
AccuracyGoal -> 8, PrecisionGoal -> 8, WorkingPrecision -> 30]
(*0.00282913941636254483746320250663*)
• N[1/(2 \[Pi]) e (1/ e 4 (ArcCot[ e] + (-1)^(1/4) Sqrt[ e] (-ArcTan[(-1)^(1/4)/Sqrt[e]] + ArcTanh[(-1)^(1/4)/Sqrt[e]])) + RootSum[1 + e^2 + 2 #1 - #1^2 - 2 #1^3 + #1^4 &, (Log[1 - #1] #1)/(-1 - #1 + #1^2) &] - RootSum[1 + e^2 + 2 #1 - #1^2 - 2 #1^3 + #1^4 &, (Log[-#1] #1)/(-1 - #1 + #1^2) &] + RootSum[1 + e^2 - 2 #1 - #1^2 + 2 #1^3 + #1^4 &, (Log[1 - #1] #1)/(-1 + #1 + #1^2) &] - RootSum[1 + e^2 - 2 #1 - #1^2 + 2 #1^3 + #1^4 &, (Log[-#1] #1)/(-1 + #1 + #1^2) &]) /. e -> 10^-5, 20] Commented Nov 10, 2020 at 13:27
• produces 1.6135641088989329018 + 0.*10^-20 I. Commented Nov 10, 2020 at 13:28