Integration of dirac_delta

edit

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)))):
        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
            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