Ask Your Question

Integration of dirac_delta

asked 2020-06-19 13:43:53 +0100

curios_mind gravatar image

updated 2020-06-19 13:44:31 +0100


When I tried to integrate the following


I get 1 (as expected)

But when I integrate


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 flag offensive close merge delete


Dirac delta is not very well handled by Maxima.

FrédéricC gravatar imageFrédéricC ( 2020-06-19 21:33:45 +0100 )edit

1 Answer

Sort by » oldest newest most voted

answered 2020-06-19 15:35:27 +0100

Florentin Jaffredo gravatar image

updated 2020-06-19 16:02:54 +0100


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)))):
        return E.integrate(x, a, b) # if not found, integrate normaly

    # If expression is a sum, call recursively
    if E.operator() is sage.symbolic.operators.add_vararg:
        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
        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)

Some tests:

sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, pi)
sage: integrate_delta(dirac_delta(cos(x)+0.5), x, 0, 6*pi)
sage: integrate_delta(dirac_delta(x^2-1), x, 0.5, 0.5)
sage: integrate_delta(dirac_delta(x^2-1), x, -1, 0.5)
sage: integrate_delta(dirac_delta(x^2-1), x, -2, 2)
edit flag offensive delete link 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.

curios_mind gravatar imagecurios_mind ( 2020-06-19 23:24:26 +0100 )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?

curios_mind gravatar imagecurios_mind ( 2020-06-20 00:05:19 +0100 )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.

Florentin Jaffredo gravatar imageFlorentin Jaffredo ( 2020-06-20 15:21:20 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools



Asked: 2020-06-19 13:43:53 +0100

Seen: 379 times

Last updated: Jun 19 '20