Trouble with the derivative of the delta function

Question

I am trying to evaluate the following integral with Mathematica:

\begin{align} I = \int_{0}^{\infty} da \, \frac{e^{-\frac{a ^2}{4s^2}} }{a^2} \mbox{sinc}\left(\tfrac{w}{2} a \right) \delta' \left( \frac{D^2}{a}- a \right), \end{align} where the prime on the delta function denotes differentiation with respect to the argument of the Delta function. When I evaluate this integral with Mathematica as:

Integrate[Exp[-a^2/(4 s^2)]/a^2 Sinc[w a / 2] Derivative[1][DiracDelta][D^2/a - a],{a,0,Infinity}, Assumptions -> s > 0 && w > 0 && D > 0]

I get the result: \begin{align} I_{Mathematica} = \frac{e^{-\frac{D^2}{4 s ^2}} }{4 D^4 s ^2 w } \left[\left(D^2+6 s ^2\right) \sin \left(\frac{D w }{2}\right)-D s ^2 w \cos \left(\frac{D w }{2}\right)\right]. \end{align}

However, if I evaluate this integral analytically, using the fact that \begin{align} \frac{d}{da} \delta\left( \frac{D^2}{a}- a \right) = - \delta' \left( \frac{D^2}{a}- a \right) \left(\frac{D^2}{a^2}+1\right) \implies \delta' \left( \frac{D^2}{a}- a \right) = - \left[\frac{d}{da} \delta\left( \frac{D^2}{a}- a \right) \right] \left(\frac{D^2}{a^2}+1\right)^{-1}, \end{align} I get the following result: \begin{align} I_{analytic} &= \int_{0}^{\infty} da \, \frac{e^{-\frac{a ^2}{4s^2}} }{a^2} \frac{\sin \left(\tfrac{w}{2} a \right) }{\tfrac{w}{2} a} \delta' \left( \frac{D^2}{a}- a \right) \\ &=- \int_{0}^{\infty} da \, \left[\frac{d}{da} \delta\left( \frac{D^2}{a}- a \right) \right] \left(\frac{D^2}{a^2}+1\right)^{-1} \frac{e^{-\frac{a ^2}{4s^2}} }{a^2} \frac{\sin \left(\tfrac{w}{2} a \right) }{\tfrac{w}{2} a} \\ &= \int_{0}^{\infty} da \, \delta\left( \frac{D^2}{a}- a \right) \left[\frac{d}{da} \left(\frac{D^2}{a^2}+1\right)^{-1} \frac{e^{-\frac{a ^2}{4s^2}} }{a^2} \frac{\sin \left(\tfrac{w}{2} a \right) }{\tfrac{w}{2} a} \right] \\ &= \int_{0}^{\infty} da \, \frac{\delta\left( D - a \right)}{2} \left[\frac{d}{da} \left(\frac{D^2}{a^2}+1\right)^{-1} \frac{e^{-\frac{a ^2}{4s^2}} }{a^2} \frac{\sin \left(\tfrac{w}{2} a \right) }{\tfrac{w}{2} a} \right] \\ &= - \frac{e^{-\frac{ D^2 }{ 4s^{2}}}}{4 D^{4} s^{2} w} \left[ \left(D^{2}+4 s^{2}\right)\sin\left( \frac{ Dw}{2} \right) - D s^{2} w \cos \left( \frac{ Dw}{2} \right) \right], \end{align} which differs from $I_{Mathematica}$ by an overall negative sign and the prefactor in front of $s^2$ in the first term.

I'm not sure if the issue is with the way Mathematica handles the derivative of the delta function or if I've made a mistake in my analytic calculation. Any help would be much appreciated, I've been staring at this for days!

Answer 1

Let's talk about the Dirac $\delta$-"function". Strictly speaking, it's a linear functional $$\delta:C^\infty(\mathbb R)\to\mathbb R\qquad\qquad\delta(f)=f(0).$$ However, we usually use the notation $$\int_{-\infty}^\infty\delta(x)f(x)dx$$ to denote the evaluation $\delta(f)$. The derivative of the $\delta$-"function" is computed via formal integration by parts: $$\delta'(f)=\int_{-\infty}^\infty\delta'(x)f(x)dx=-\int_{-\infty}^\infty\delta(x)f'(x)dx=-f'(0).$$ Your integral has the additional complications that there is a function inside the argument of $\delta'(x)$, and that the integral is not taken over all of $\mathbb R$. Composing distributions with functions is, in general, not possible, but in this case we can appeal to a theorem of Hormander:

Theorem: Suppose $f:M\to N$ is a smooth function whose differential is everywhere surjective. Then there is a linear map $f^*:\mathscr D(N)\to\mathscr D(M)$ such that $f^*u=u\circ f$ for all $u\in C(N)$.

For our purposes, this means $\int_{-\infty}^\infty\delta'(f(x))g(x)dx$ makes sense provided $f(x)$ is smooth and $f'(x)$ never vanishes. Similarly, reducing the domain of integration is, in general, not possible, but we have:

Theorem Suppose $E_1$ and $E_2$ are disjoint closed sets, and let $\mathscr D_{E_i}$ denote the set of distributions which coincide with a smooth function on $E_i^c$ for $i=1,2$. Then there is a bilinear map $$m:\mathscr D_{E_1}\times\mathscr D_{E_2}\to\mathscr D(\mathbb R^n)$$ such that $m(u,v)=uv$ when $u$ and $v$ are continuous.

In our case, we would like to compute the integral $$\int_0^\infty\delta'\left(\frac{D^2}{x}-x\right)g(x)dx=\int_{-\infty}^\infty\chi_{(0,\infty)}(x)\delta'\left(\frac{D^2}{a}-a\right)g(x)dx,$$ where $\chi_{(0,\infty)}$ is the characteristic function of the half-line $(0,\infty)$. The theorem says that the product $$\chi_{(0,\infty)}(x)\delta'\left(\frac{D^2}{x}-x\right)$$ makes sense whenever the singular support of $\chi_{(0,\infty)}$, namely $\{0\}$, does not intersect the singular support of $\delta'\left(\frac{D^2}{x}-x\right)$, namely $\{D,-D\}$. Thus when $D\neq 0$, our integral makes sense and $$\int_0^\infty\delta'\left(\frac{D^2}{x}-x\right)g(x)dx=\begin{cases}g'(D),&D>0\\g(-D),&D<0\end{cases}.$$ To compute your integral, just plug in your particular function $g(x)$. When you're working with distributions (like $\delta$) you need to be very careful about what you do with them. I don't know how Mathematica conceptualizes the $\delta$-distribution, but I wouldn't be inclined trust that it would go through the necessary analytical reasoning and get the right answer.

TL;DR: Do your distributional calculus by hand.

Answer 2

Here my attempt to solve the integral Integrate[f[a] Derivative[1][DiracDelta][d^2/a - a],{a,0,Infinity}]:

f[a_] := Exp[-a^2/(4 s^2)]/a^2 Sinc[w a/2] 

Substitution u[a]=d^2/a-a (integrationlimits change to u[0]=Infinity],u[Infinity]=-Infinity)

u[a_] := d^2/a - a

sola = Solve[u == d^2/a - a, a][[2]] (*solution a>0*)

Now Mathematica is able to solve the integral

int=Integrate[f[a/.sola] Derivative[1][DiracDelta][u]/u'[a]/.sola ,{u, Infinity,-Infinity}]

(*(E^(-(d^2/(4 s^2))) (d s^2 w Cos[(d w)/2] - (d^2 + 4 s^2) Sin[(d w)/2]))/(4 d^3 Sqrt[d^2] s^2 w)*)

Hope it helps solving your problem!

Answer 3

My previous answer and comments were wrong. I didn't notice the argument of the δ function was not linear in the integration variable (and I wasn't even drunk).

In the Wikipedia page, there is this paragraph

In the integral form the generalized scaling property may be written as $∫_{-∞}^∞ f ( x ) δ ( g ( x ) ) d x = ∑_i f ( x_i ) / | g ′ ( x_i ) | $.

The Jacobian of the transformation is 1/g'(x). Please note the absolute value in the denominator.

Basically, find the zeroes of the argument of the δ, and integrate around them (by parts if necessary). Also,

The distributional derivative of the Dirac delta distribution is the distribution δ′ defined on compactly supported smooth test functions φ by $δ ′ [ φ ] = − δ [ φ ′ ] = − φ ′ ( 0 )$ .

(1) Finding the zeroes:

Solve[-a + Δ^2/a == 0, a]

(2) Finding the Jacobian:

jac = Solve[Dt[-a + Δ^2/a == u[a]], u'[a]] 
  /. Dt[Δ] → 0 /. a → Δ // FullSimplify

(3) Evaluating the integral by parts (don't forget the minus sign in front):

v1 = -D[(E^(-(a^2/(4 s^2))) Sinc[(a w)/2])/a^2, a] / Abs[jac]
  /. a → Δ // FullSimplify

(4) Divide the integral by the Jacobian (the previous division was because of the integration by parts, this one because of the scaling):

v1 / Abs[ jac ]

The answer is the same as $I_{MMA}$. By the way, MMA is simply using

Integrate[ f[a] DiracDelta'[2 a], {a, -∞, ∞}]
(*-(1/4) f'[0]*)

Prove that analytically and you will find the error in your analytic calculation.

Answer 4

speculation: Mathematica cannot handle Derivative[1][DiracDelta][1/x-x]in a right way?

Here I'll give a simplified example which perhaps shows that Mathematica gives a wrong result, when applied to Derivative[1][DiracDelta][1/x-x]!

Let's consider the integral

Integrate[Derivative[1][DiracDelta][1/x - x], {x, 0, Infinity} ]
(*0*)

which MMA (v12) evaluates to zero!

Alternatively integration with substitution u=1/x-x, x=-u/2+Sqrt[1+(u/2)^2] (see my first answer)

us=D[1/x-x,x]/. x->-u/2+Sqrt[1+(u/2)^2];
Integrate[ Derivative[1][DiracDelta][u]/us, {u, Infinity,-Infinity} ]
(*1/4*)

To "proof" the last result I'll consider the deltadistribution as a well known limit

dirac = Function[x, Exp[-(x^2/(2 eps))]/Sqrt[2 Pi eps]] (* eps->0 *)

int=Integrate[dirac'[1/x - x], {x, 0, Infinity} ]

(*(E^(1/eps) (-BesselK[0, 1/eps] + BesselK[1, 1/eps]))/(eps^(3/2) Sqrt[2 \[Pi]])*)     

eps->0

Simplify[ Normal[Series[int, {eps, 0, 0}]], eps > 0]

(*1/4*)

Why can't Mathematica find this result? What's wrong here?