# Revision history [back]

Hello,

I don't know why this is the case, but here is a dirty workaround:

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.
"""

t = var('temporary_variable')
E *= t
# 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))):

match = E.match(SR.wild(0)*dirac_delta(SR.wild(1))) # match the pattern
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
roots_with_multiplicity = []

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
roots_with_multiplicity.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]
roots_with_multiplicity += [r[x].subs({z: i}) for i in IntegerRange(ceil(a0.lhs()), ceil(b0.rhs()))]

# returns sum(f(r)/abs(g'(r))) where the sum is over the roots of g.
return sum(f.subs({x: r, t: 1})/abs(g.diff(x).subs({x: r, t: 1})) for r in roots_with_multiplicity 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. But I tried it with some trigonometric functions, and it seems fine. 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


Hello,

I don't know why this is the case, but here is a dirty workaround:

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.
"""

t = var('temporary_variable')
E *= t
# 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))):
(E.find(SR.wild(0)*dirac_delta(SR.wild(1))) or E.find(dirac_delta(SR.wild(1)))):

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
roots_with_multiplicity = []

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
roots_with_multiplicity.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]
roots_with_multiplicity += [r[x].subs({z: i}) for i in IntegerRange(ceil(a0.lhs()), ceil(b0.rhs()))]

# returns sum(f(r)/abs(g'(r))) sum(f(r)*1/abs(g'(r))) where the sum is over the roots of g.
return sum(f.subs({x: r, t: 1})/abs(g.diff(x).subs({x: r, t: 1})) r})/abs(g.diff(x).subs({x: r})) for r in roots_with_multiplicity 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. But I tried it with some trigonometric functions, and it seems fine. 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


Hello,

I don't really know why this is the case, but here simplifying Dirac deltas is a dirty workaround: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)))):

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
roots_with_multiplicity = []

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
roots_with_multiplicity.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]
roots_with_multiplicity += [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 roots_with_multiplicity 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. But I tried it with some trigonometric functions, and it seems fine. 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


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
roots_with_multiplicity = []

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
roots_with_multiplicity.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]
roots_with_multiplicity += [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 roots_with_multiplicity 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


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
roots_with_multiplicity 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
roots_with_multiplicity.append(r[x])
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]
roots_with_multiplicity 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 roots_with_multiplicity 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