I am having difficulty computing with delta functions in Mathematica.

I want to do the following:

$\int_0^R dz \int_0^a d\rho e^{-k |z|}* \rho * J_1(k \rho)* \delta(\rho-a)[\delta(z - R/2) + \delta(z + R/2)]$

$J_1(k \rho)$ is just BesselJ[1, k*ρ].

I don't know how to correctly evaluate $\delta(\rho-a)$ for example.

I expected the result of the following to be 2;

Integrate[x*DiracDelta[x - 2], {x, 0, 2}]

Instead, I get:

 2 HeavisideTheta[0]

When I execute:

Integrate[ρ*BesselJ[1, k*ρ] *Exp[-k*Abs[z]] * 
  DiracDelta[ρ - a]*(DiracDelta[z - (R/2)] + 
    DiracDelta[z + (R/2)]), {z, 0, R}, {ρ, 0, a}]

I get:

 a E^(-(1/2) k Abs[R])
   BesselJ[1, a k] (-1 + 2 UnitStep[a]) (-1 + 2 UnitStep[R]), 
 a ∈ Reals && R ∈ Reals]

enter image description here

Can someone please help me with this?

  • 3
    $\begingroup$ Why would you expect that integral to be 2? Your limit of integration is exactly on top of the delta "peak", so it makes sense for the value to be somewhat undefined. (By the chain rule and using the delta functions distributional antiderivative HeavisideTheta, the integral is exactly HeavisideTheta[0], which is what is returned) $\endgroup$
    – Lukas Lang
    Apr 23, 2023 at 20:15
  • 1
    $\begingroup$ This integral makes no sense in traditional math. That was many times noticed at this forum. $\endgroup$
    – user64494
    Apr 24, 2023 at 5:26

3 Answers 3


(* "13.2.1 for Mac OS X ARM (64-bit) (January 27, 2023)" *)


Extend the limits of the integration by Infinitesimally small ϵ > 0

Assuming[ϵ > 0, 
 Integrate[x*DiracDelta[x - 2], {x, 0, 2 + ϵ}]]

(* 2 *)

Limit[Assuming[R > 0 && ϵ > 0, 
   Integrate[ρ*BesselJ[1, k*ρ]*Exp[-k*Abs[z]]*
     DiracDelta[ρ - a]*(DiracDelta[z - (R/2)] + 
       DiracDelta[z + (R/2)]), {z, 0, R}, {ρ, 0, 
     a + ϵ}]], ϵ -> 0] // Simplify

enter image description here

  • $\begingroup$ The result depends on the sign of $\varepsilon$, using $-\varepsilon$ instead of $+\varepsilon% leads to 0 for the first integral and \begin{cases} -a e^{-\frac{k R}{2}} J_1(a k) & a<0 \\ 0 & \text{True} \end{cases} for the second $\endgroup$ Apr 25, 2023 at 14:56
  • 2
    $\begingroup$ @userrandrand - He specifically indicated that he wanted the result 2 so I specially indicated that 𝜀 > 0 $\endgroup$
    – Bob Hanlon
    Apr 25, 2023 at 15:26

As stated by Lukas Lang, you have a limit of integration right on the DiracDelta "peak" and Mathematica does not know whether you mean the peak to be within or outside the limit of integration. In the limit HeavisideTheta[0] can be either $0$ or $1$.

Limit[HeavisideTheta[x], x -> 0, Direction -> "FromAbove"]
(* 1 *)

Limit[HeavisideTheta[x], x -> 0, Direction -> "FromBelow"]
(* 0 *)

so that HeavisideTheta[0] is essentially $0$ or$1$ at least in the limit and it is up to the user to tell Mathematica which he means.

easily done by:

Integrate[x*DiracDelta[x - 2], {x, 0, 2}] /. HeavisideTheta[0] -> 1
(* 2 *)

Older versions of Mathematica used to return $1/2$ instead of HeavisideTheta[0] which was always wrong.

For me, your first integral comes out clean just by making the

$Assumptions = a > 0 && R > 0
  • $\begingroup$ I can not seem to find find a good reference but I believe that in the path integral formulation of the Langevin equation, the choice HeavisideTheta[0]=0 corresponds to the Ito prescription as shown here: arxiv.org/pdf/1106.4129.pdf on page 8 and I believe (but I am not sure, it's a vague memory) that HeavisideTheta[0]=1/2 corresponds to the Stratonovich prescription. If that is the case then HeavisideTheta[0]=1/2 is not always wrong. $\endgroup$ Apr 25, 2023 at 15:58
  • $\begingroup$ There is also a discussion about why HeavisideTheta[0]=1/2 might seem natural in some cases here: physics.stackexchange.com/questions/61176/… $\endgroup$ Apr 25, 2023 at 16:01
  • 1
    $\begingroup$ @userrandrand I have seen those discussions of HeavisideTheta=1/2 and maybe someone can think of a real problem where that is appropriate. I'm guessing it comes from the idea that if you compute a Fourier series for the HeavisideTheta, you get 1/2 at the transition, and maybe that's where Mathematica got it. My point was that in old versions, Mathematica did not mention HeavisideTheta, but in integrals that had one of the limits on the DiracDelta peak, it just returned $1/2$, which for me was always wrong and I never got a choice. Now the user can choose an appropriate value. $\endgroup$
    – Bill Watts
    Apr 25, 2023 at 20:56

The integral factorizes

Integrate[exp[-Abs[k] z] (delta[z+R/2]+delta[z-R/2]),{z,0,R]]* 
Integrate[J[1, k r] delta[r-a],{r,0,a}]

The first integral extends over

(0,R > 0) -> Exp[ - Abs[k] R/2 ] 

or with direction reversed over

{R/2<0, 0}  ->  - Exp[- Abs[k] R/2] =Exp[ Abs[k  R/2 ] + I pi]

So the result is a simple exponential function multplied by Sign[R], that equals the Mathematica result

E^( R/2) DiscreteDelta[R] + (2 HeavisideTheta[R] - 1) E^(-( R/2))

== Sign[R] Exp[- |k| R/2 ]

Plot[ E^( R/2) DiscreteDelta[R] + (2 HeavisideTheta[R] - 1) E^(-( R/
2)), {R, -10^.4, 1}, Epilog -> {Disk[{0, 0}, 0.1]}]

sign *exp

Because of this ditochomy the point R=0 has to be left out, but since the exp-function is bounded, the integral has to be set to zero by zero length

I don't touch the question

Integrate[ r J[1, k r] delta[r-a] ,{r,0,a}]

because it represents the radial part of a 2d integral over the inner of a circle.

If thats the case, the coordinate DiracDelta is a wrong interpretation of the euclidean DeltaFunction on the circle in the plane as the limit of a normed Gaussian measure supported along the circle at radius a with length 2 pi a.

Integrate[f[x,y] 1/N[sigma] Exp[-(Sqrt[x^2+y^2]-a)^2/sigma^2 ]] 
, x^2 + y^2 < ? <= ? >  a^2 ]  

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.