ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sat, 20 Jun 2020 15:21:20 +0200Integration of dirac_deltahttps://ask.sagemath.org/question/52088/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?Fri, 19 Jun 2020 13:43:53 +0200https://ask.sagemath.org/question/52088/integration-of-dirac_delta/Comment by FrédéricC for <p>Hello,</p>
<p>When I tried to integrate the following</p>
<pre><code>(dirac_delta(x)).integrate(x,-pi,pi)
</code></pre>
<p>I get <code>1</code> (as expected)</p>
<p>But when I integrate</p>
<pre><code>(dirac_delta(cos(x))).integrate(x,-pi,pi)
</code></pre>
<p>I get <code>integrate(dirac_delta(cos(x)), x, -pi, pi)</code> back. I would expect to get <code>2</code> 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</p>
<p>If I try the same with Mathematica,</p>
<pre><code>Integrate[DiracDelta[x], {x, -Pi, Pi}]
</code></pre>
<p>I get <code>1</code>
and when I try</p>
<pre><code>Integrate[DiracDelta[Cos[x]], {x, -Pi, Pi}]
</code></pre>
<p>I get <code>2</code> as I would expect.</p>
<p>What is the reason of this behaviour of Sage?</p>
https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?comment=52106#post-id-52106Dirac delta is not very well handled by Maxima.Fri, 19 Jun 2020 21:33:45 +0200https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?comment=52106#post-id-52106Answer by Florentin Jaffredo for <p>Hello,</p>
<p>When I tried to integrate the following</p>
<pre><code>(dirac_delta(x)).integrate(x,-pi,pi)
</code></pre>
<p>I get <code>1</code> (as expected)</p>
<p>But when I integrate</p>
<pre><code>(dirac_delta(cos(x))).integrate(x,-pi,pi)
</code></pre>
<p>I get <code>integrate(dirac_delta(cos(x)), x, -pi, pi)</code> back. I would expect to get <code>2</code> 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</p>
<p>If I try the same with Mathematica,</p>
<pre><code>Integrate[DiracDelta[x], {x, -Pi, Pi}]
</code></pre>
<p>I get <code>1</code>
and when I try</p>
<pre><code>Integrate[DiracDelta[Cos[x]], {x, -Pi, Pi}]
</code></pre>
<p>I get <code>2</code> as I would expect.</p>
<p>What is the reason of this behaviour of Sage?</p>
https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?answer=52093#post-id-52093Hello,
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)
1Fri, 19 Jun 2020 15:35:27 +0200https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?answer=52093#post-id-52093Comment by curios_mind for <p>Hello,</p>
<p>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.</p>
<p>Here is a version designed to work with trigonometric functions. It probably breaks if no solution is found, so be careful.</p>
<pre><code>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)
</code></pre>
<p>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:</p>
<pre><code>sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2
</code></pre>
<p>Some tests:</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?comment=52113#post-id-52113Thank 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.Fri, 19 Jun 2020 23:24:26 +0200https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?comment=52113#post-id-52113Comment by curios_mind for <p>Hello,</p>
<p>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.</p>
<p>Here is a version designed to work with trigonometric functions. It probably breaks if no solution is found, so be careful.</p>
<pre><code>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)
</code></pre>
<p>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:</p>
<pre><code>sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2
</code></pre>
<p>Some tests:</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?comment=52118#post-id-52118Unfortunately, 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?](https://ask.sagemath.org/question/52117/dirac_delta-integration-possible-bug/)Sat, 20 Jun 2020 00:05:19 +0200https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?comment=52118#post-id-52118Comment by Florentin Jaffredo for <p>Hello,</p>
<p>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.</p>
<p>Here is a version designed to work with trigonometric functions. It probably breaks if no solution is found, so be careful.</p>
<pre><code>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)
</code></pre>
<p>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:</p>
<pre><code>sage: integrate_delta(dirac_delta(cos(x)), x, -pi, pi)
2
</code></pre>
<p>Some tests:</p>
<pre><code>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
</code></pre>
https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?comment=52127#post-id-52127If 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.Sat, 20 Jun 2020 15:21:20 +0200https://ask.sagemath.org/question/52088/integration-of-dirac_delta/?comment=52127#post-id-52127