# Why are the results of integrations of a DiracDelta related integrand and its approximations strikingly different?

Let us consider a double integral, treating it as a functional which corresponds a number to a function,

Integrate[ Exp[2 x + 3 y]*DiracDelta[x - UnitStep[x + y]], {x, -3/2, 3/2}, {y, -3/2, 3/2}]


1/E^(9/2) - 2 E^2 + E^(13/2)

In fact, Mathematica calculates the iterated integral, so we also consider another iterated integral

Integrate[Exp[2 x + 3 y]*DiracDelta[x - UnitStep[x + y]], {y, -3/2,3/2}, {x, -3/2, 3/2}]


-((2 DiracDelta)/E^(9/2))

As we see, the results are strikingly different. Let us consider approximations of DiracDelta[x - UnitStep[x + y]] in the weak topology (Exactly saying, the functions which are associated with that approximations. I'd like to recall that the $$\delta$$-distribution is not associated with any usual function.) ant the limits when those approximations tend to DiracDelta[x - UnitStep[x + y]] in the weak topology:

Integrate[ Exp[2 x + 3 y]*UnitBox[(x - UnitStep[x + y])/eps]/eps, {x, -3/2, 3/2},
{y, -3/2, 3/2}, Assumptions -> eps > 0 && eps < 1]


(-2 E^(7/2) Sinh[eps/2] + 2 E^(9/2) Sinh[eps/2] - Sinh[eps] + E^11 Sinh[eps])/(3 E^(9/2) eps)

Limit[%, eps -> 0, Direction -> "FromAbove"]


(-1 - E^(7/2) + E^(9/2) + E^11)/(3 E^(9/2))

and

Integrate[Exp[2 x + 3 y]*UnitBox[(x - UnitStep[x + y])/eps]/eps, {y, -3/2, 3/2},
{x, -3/2, 3/2}, Assumptions -> eps > 0 && eps < 1]


1/(3 eps) E^(-(9/2) - (3 eps)/2) (-2 E^(7/2 + (3 eps)/2) Sinh[eps/2] + 3 E^(9/2 + (3 eps)/2) Sinh[eps/2] + E^(7/2 + (5 eps)/2) Sinh[eps/2] + E^(7/2 + (7 eps)/2) Sinh[eps/2] + E^(3/2) Sinh[eps] + E^(9/2) Sinh[eps] + E^(13/2) Sinh[eps] - E^(3 eps/2) Sinh[eps] + E^(11 + (3 eps)/2) Sinh[eps] - E^(3/2 + 3 eps) Sinh[eps] - E^(7/2 + 3 eps) Sinh[eps] - E^(13/2 + 3 eps) Sinh[eps] - E^(9/2 + eps/2) Sinh[(3 eps)/2] + 2 E^(3/2 + (3 eps)/2) Sinh[eps] Sinh[(3 eps)/2] + 2 E^(13/2 + (3 eps)/2) Sinh[eps] Sinh[(3 eps)/2])

Limit[%, eps -> 0, Direction -> "FromAbove"]


(-1 - E^(7/2) + E^(9/2) + E^11)/(3 E^(9/2)).

To be sure, let us also calculate those with non-symmetric approximations:

Integrate[Exp[2 x + 3 y]* UnitBox[(x - UnitStep[x + y] - Sqrt[eps])/eps]/eps, {x, -3/2,3/2},
{y, -3/2, 3/2}, Assumptions -> eps > 0 && eps < 1];
Limit[%, eps -> 0, Direction -> "FromAbove"]


(-1 - E^(7/2) + E^(9/2) + E^11)/(3 E^(9/2))

and

eps = 0.005; NIntegrate[Exp[2*x + 3*y]*eps/((x - UnitStep[x + y] - eps)^2 + eps^2)/Pi,
{y, -3/2, 3/2}, {x, -3/2, 3/2}]


223.53

(compare the latter with N[(-1 - E^(7/2) + E^(9/2) + E^11)/(3 E^(9/2))] which results in 221.921 and with N[1/E^(9/2) - 2 E^2 + E^(13/2)] which results in 650.375).

It should be noticed the case under consideration is substantially different from that case: two iterated integrals of  Exp[2 x + 3 y]*DiracDelta[x - UnitStep[x + y]] and limits of the integrals of its approximation in the weak topology Exp[2 x + 3 y]*UnitBox[(x - UnitStep[x + y])/eps]/eps as eps tends to zero from above produce three different results; if the limits of the integration are slightly perturbed, then the intersection of the support of DiracDelta[x - UnitStep[x + y]] with the boundary of the set of the integration does not qualitativily change as

ContourPlot[x - UnitStep[x + y] == 0, {x, -2, 2}, {y, -2, 2}] shows.

The question is open: what does the result of the integration of Exp[2 x + 3 y]*DiracDelta[x - UnitStep[x + y]] over the square {x, -3/2, 3/2}, {y, -3/2, 3/2} mean? The only reference I know where the integration of distributions over bounded sets is defined is Antosik, P., Mikusinski, J., Sikorski, R. Theory of distributions. The sequential approach. Reprints. (English) Zbl 0267.46028 Amsterdam: Elsevier Scientific Publishing Company; Warszawa: PWN-Polish Scientific Publishers. XIV, 273 p.(1973), but a smooth integral introduced there does not include the integral under consideration.

• Somebody suggests to close the question as based on options. Which options do the work? Also a downvote is not motivated. What is wrong in my question? – user64494 Mar 10 at 13:40
• DiracDelta is a tool of mathematics as applied to physical situations. What real-world object are you attempting to model here? – John Doty Mar 10 at 17:38
• @JohnDogy: Thank you for your constant interest to my questions. I don't know such an object. I simply study two syntactically correct commands of Mathematica which produce doubtful results and ask about it in the forum devoted to Mathematica. Constructive replies of you on the topic are welcome. – user64494 Mar 10 at 18:04
• Doubtful results from ill-posed problems are normal. It may be syntactically correct to attempt to pick up a chain saw by its blade, but the outcome will not be pretty. – John Doty Mar 10 at 18:41
• I don't think Fubini's theorem holds on the original singular integrand. One order gives a wrong result, possibly GIGO or possibly a bug. – Daniel Lichtblau Mar 10 at 19:40

"Integral[ Exp[2 x + 3 y]*DiracDelta[x - UnitStep[x + y]], {x, -3/2, 3/2}, {y, -3/2, 3/2}]"
==Integrate[Exp[2 0 + 3 y], {y, -3/2, 0}]
+Integrate[Exp[2 1 + 3 y], {y, -1, 3/2}]


This can be seen by examining the argument of DiracDelta in your example:

ContourPlot[x - UnitStep[x + y] == 0, {x, -3/2, 3/2},{y, -3/2, 3/2}, FrameLabel -> {x, y} ] It's obvious to first integrate "over x" using the rule Integrate[DiracDelta[x-p] f[x],{x,p-eps,x+eps}]==f[p]

The value of the integral follows to

Integrate[Exp[2 0 + 3 y], {y, -3/2, 0}]+Integrate[Exp[2 1 + 3 y], {y, -1, 3/2}]
(*1/3 - 1/(3 E^(9/2)) + (-1 + E^(15/2))/(3 E)==221.921*)


It might be "confirmed" numerical with

dirac = Function[{x, eps}, Which[-eps < x < eps, 1/(2eps), True, 0]]
NIntegrate[Exp[2 x + 3 y]*dirac[x - UnitStep[x + y], .001], {x,-3/2,3/2}, {y, -3/2, 3/2}]
(*221.921*)


final note: Direct evaluation with Mathematica

Integrate[ Exp[2 x + 3 y]*DiracDelta[x -UnitStep[x + y]], {x, -3/2, 3/2}, {y, -3/2,3/2}]
(*1/E^(9/2) - 2 E^2 + E^(13/2)==650.375*)


gives wrong result. Perhaps because argument of DiracDelta also vanishs on the boundary of the integration range...

• Thank you. However, the equality Integrate[ Exp[2 x + 3 y]*DiracDelta[x - UnitStep[x + y]], {x, -3/2, 3/2}, {y, -3/2, 3/2}] ==Integrate[Exp[2 0 + 3 y], {y, -3/2, 0}] +Integrate[Exp[2 1 + 3 y], {y, -1, 3/2}] is not grounded: Integrate[ Exp[2 x + 3 y]*DiracDelta[x - UnitStep[x + y]], {x, -3/2, 3/2}, {y, -3/2, 3/2}] performs a different value. – user64494 Mar 10 at 10:32
• If you read my complete answer you might found the "ground" – Ulrich Neumann Mar 10 at 10:33
• No, I don't find arguments there though there's something about it. Thank you anyway. – user64494 Mar 10 at 10:44
• It seems surprising that my analytical result gives the same value as your numerical approximations! Let me know if someone else except you doesn't find arguments in my answer... – Ulrich Neumann Mar 10 at 10:44
• UlrichNeumann(@ does not work.): Yuriy Golovatyi. You may contact him with his kind permission. – user64494 Mar 10 at 14:01