# Integration of dirac_delta

Hello,

When I tried to integrate the following

(dirac_delta(x)).integrate(x,-pi,pi)


I get 1 (as expected)

But when I integrate

(dirac_delta(cos(x))).integrate(x,-pi,pi)


I get integrate(dirac_delta(cos(x)), x, -pi, pi) back. I would expect to get 2 since there are two points which are zero between the upper and lower bounds of the integral. So by any means of approaching these points with limit during the integration, the integral of dirac_delta should have returned 1 for each point which adds up to 2

If I try the same with Mathematica,

Integrate[DiracDelta[x], {x, -Pi, Pi}]


I get 1 and when I try

Integrate[DiracDelta[Cos[x]], {x, -Pi, Pi}]


I get 2 as I would expect.

What is the reason of this behaviour of Sage?

edit retag close merge delete

Dirac delta is not very well handled by Maxima.

( 2020-06-19 21:33:45 +0200 )edit

Sort by » oldest newest most voted

Hello,

I don't really know why this is the case, but simplifying Dirac deltas is a hard problem, since solving an equation is always needed.

Here is a version designed to work with trigonometric functions. It probably breaks if no solution is found, so be careful.

from operator import eq
def integrate_delta(E, x, a, b):
"""
Integrate E over [a, b] with respect to x, but tries to find a dirac delta first.
"""
# try to find a pattern of the shape f(x)*delta(g(x))
if not (E.find(SR.wild(0)*dirac_delta(SR.wild(1))) or E.find(dirac_delta(SR.wild(1)))):

# If expression is a sum, call recursively
return sum(integrate_delta(op, x, a, b) for op in E.operands())

match = E.match(SR.wild(0)*dirac_delta(SR.wild(1))) # match the pattern
if match is None:
match = E.match(dirac_delta(SR.wild(1)))
match[SR.wild(0)] = 1

f, g = match[SR.wild(0)], match[SR.wild(1)] # extract the functions f and g

roots = list(g.solve(x, solution_dict=True, to_poly_solve="force")) # Solve g(x)==0
explicit_roots = []

for r in roots:
# sometimes there is an unknown integer in the solutions
extra_variables = set(r[x].variables()).difference(g.variables())
if len(extra_variables)==0: # if not, keep the root
explicit_roots.append(r[x])
continue
if len(extra_variables)>1: # if there are more, we don't know how to solve, fall back.
return E.integrate(x, a, b)
z = extra_variables.pop()

# if there is, we must make sure to keep only the values that stay inside the bounds
a0, b0 = [t1 for t1 in solve_ineq([r[x]>=a, r[x]<=b], (z,))
if all(t2.operator() is not eq for t2 in t1)][0]
explicit_roots += [r[x].subs({z: i}) for i in IntegerRange(ceil(a0.lhs()), ceil(b0.rhs()))]

# returns sum(f(r)*1/abs(g'(r))) where the sum is over the roots of g.
return sum(f.subs({x: r})/abs(g.diff(x).subs({x: r})) for r in explicit_roots if r>=a and r<=b)


This function will in most case call the usual integrate method, but will first try to find a Dirac delta. The code is absolutely ugly, and there is no guarantee whatsoever that it always works. Use it like this:

sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2


Some tests:

sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, pi)
2/3*sqrt(3)
sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, 6*pi)
4*sqrt(3)
sage: integrate_delta(dirac_delta(x^2-1), x, 0.5, 0.5)
0
sage: integrate_delta(dirac_delta(x^2-1), x, -1, 0.5)
1/2
sage: integrate_delta(dirac_delta(x^2-1), x, -2, 2)
1

more

Thank you very much. The code works with the trigonometric expressions (at least for the ones I tested).

I am surprised that this was not implemented yet in sage.

( 2020-06-19 23:24:26 +0200 )edit

Unfortunately, I think there are more bugs in Sage's integral with dirac_delta. I opened up a new thread for them. Another: dirac_delta integration: possible bug?

( 2020-06-20 00:05:19 +0200 )edit

If you start using symbolic variable inside the Dirac delta, then the problem is A LOT harder, because now the roots can be anywhere, depending on the value of the free parameter. I don't think maxima is powerful enough to handle it. But I think there should at least be a NotImplemented error, instead of a silent wrong answer.

( 2020-06-20 15:21:20 +0200 )edit