202 lines
7.3 KiB
Python
202 lines
7.3 KiB
Python
|
from sympy.core.mul import Mul
|
||
|
|
from sympy.core.singleton import S
|
||
|
|
from sympy.core.sorting import default_sort_key
|
||
|
|
from sympy.functions import DiracDelta, Heaviside
|
||
|
|
from .integrals import Integral, integrate
|
||
|
|
|
||
|
|
|
||
|
|
def change_mul(node, x):
|
||
|
|
"""change_mul(node, x)
|
||
|
|
|
||
|
|
Rearranges the operands of a product, bringing to front any simple
|
||
|
|
DiracDelta expression.
|
||
|
|
|
||
|
|
Explanation
|
||
|
|
===========
|
||
|
|
|
||
|
|
If no simple DiracDelta expression was found, then all the DiracDelta
|
||
|
|
expressions are simplified (using DiracDelta.expand(diracdelta=True, wrt=x)).
|
||
|
|
|
||
|
|
Return: (dirac, new node)
|
||
|
|
Where:
|
||
|
|
o dirac is either a simple DiracDelta expression or None (if no simple
|
||
|
|
expression was found);
|
||
|
|
o new node is either a simplified DiracDelta expressions or None (if it
|
||
|
|
could not be simplified).
|
||
|
|
|
||
|
|
Examples
|
||
|
|
========
|
||
|
|
|
||
|
|
>>> from sympy import DiracDelta, cos
|
||
|
|
>>> from sympy.integrals.deltafunctions import change_mul
|
||
|
|
>>> from sympy.abc import x, y
|
||
|
|
>>> change_mul(x*y*DiracDelta(x)*cos(x), x)
|
||
|
|
(DiracDelta(x), x*y*cos(x))
|
||
|
|
>>> change_mul(x*y*DiracDelta(x**2 - 1)*cos(x), x)
|
||
|
|
(None, x*y*cos(x)*DiracDelta(x - 1)/2 + x*y*cos(x)*DiracDelta(x + 1)/2)
|
||
|
|
>>> change_mul(x*y*DiracDelta(cos(x))*cos(x), x)
|
||
|
|
(None, None)
|
||
|
|
|
||
|
|
See Also
|
||
|
|
========
|
||
|
|
|
||
|
|
sympy.functions.special.delta_functions.DiracDelta
|
||
|
|
deltaintegrate
|
||
|
|
"""
|
||
|
|
|
||
|
|
new_args = []
|
||
|
|
dirac = None
|
||
|
|
|
||
|
|
#Sorting is needed so that we consistently collapse the same delta;
|
||
|
|
#However, we must preserve the ordering of non-commutative terms
|
||
|
|
c, nc = node.args_cnc()
|
||
|
|
sorted_args = sorted(c, key=default_sort_key)
|
||
|
|
sorted_args.extend(nc)
|
||
|
|
|
||
|
|
for arg in sorted_args:
|
||
|
|
if arg.is_Pow and isinstance(arg.base, DiracDelta):
|
||
|
|
new_args.append(arg.func(arg.base, arg.exp - 1))
|
||
|
|
arg = arg.base
|
||
|
|
if dirac is None and (isinstance(arg, DiracDelta) and arg.is_simple(x)):
|
||
|
|
dirac = arg
|
||
|
|
else:
|
||
|
|
new_args.append(arg)
|
||
|
|
if not dirac: # there was no simple dirac
|
||
|
|
new_args = []
|
||
|
|
for arg in sorted_args:
|
||
|
|
if isinstance(arg, DiracDelta):
|
||
|
|
new_args.append(arg.expand(diracdelta=True, wrt=x))
|
||
|
|
elif arg.is_Pow and isinstance(arg.base, DiracDelta):
|
||
|
|
new_args.append(arg.func(arg.base.expand(diracdelta=True, wrt=x), arg.exp))
|
||
|
|
else:
|
||
|
|
new_args.append(arg)
|
||
|
|
if new_args != sorted_args:
|
||
|
|
nnode = Mul(*new_args).expand()
|
||
|
|
else: # if the node didn't change there is nothing to do
|
||
|
|
nnode = None
|
||
|
|
return (None, nnode)
|
||
|
|
return (dirac, Mul(*new_args))
|
||
|
|
|
||
|
|
|
||
|
|
def deltaintegrate(f, x):
|
||
|
|
"""
|
||
|
|
deltaintegrate(f, x)
|
||
|
|
|
||
|
|
Explanation
|
||
|
|
===========
|
||
|
|
|
||
|
|
The idea for integration is the following:
|
||
|
|
|
||
|
|
- If we are dealing with a DiracDelta expression, i.e. DiracDelta(g(x)),
|
||
|
|
we try to simplify it.
|
||
|
|
|
||
|
|
If we could simplify it, then we integrate the resulting expression.
|
||
|
|
We already know we can integrate a simplified expression, because only
|
||
|
|
simple DiracDelta expressions are involved.
|
||
|
|
|
||
|
|
If we couldn't simplify it, there are two cases:
|
||
|
|
|
||
|
|
1) The expression is a simple expression: we return the integral,
|
||
|
|
taking care if we are dealing with a Derivative or with a proper
|
||
|
|
DiracDelta.
|
||
|
|
|
||
|
|
2) The expression is not simple (i.e. DiracDelta(cos(x))): we can do
|
||
|
|
nothing at all.
|
||
|
|
|
||
|
|
- If the node is a multiplication node having a DiracDelta term:
|
||
|
|
|
||
|
|
First we expand it.
|
||
|
|
|
||
|
|
If the expansion did work, then we try to integrate the expansion.
|
||
|
|
|
||
|
|
If not, we try to extract a simple DiracDelta term, then we have two
|
||
|
|
cases:
|
||
|
|
|
||
|
|
1) We have a simple DiracDelta term, so we return the integral.
|
||
|
|
|
||
|
|
2) We didn't have a simple term, but we do have an expression with
|
||
|
|
simplified DiracDelta terms, so we integrate this expression.
|
||
|
|
|
||
|
|
Examples
|
||
|
|
========
|
||
|
|
|
||
|
|
>>> from sympy.abc import x, y, z
|
||
|
|
>>> from sympy.integrals.deltafunctions import deltaintegrate
|
||
|
|
>>> from sympy import sin, cos, DiracDelta
|
||
|
|
>>> deltaintegrate(x*sin(x)*cos(x)*DiracDelta(x - 1), x)
|
||
|
|
sin(1)*cos(1)*Heaviside(x - 1)
|
||
|
|
>>> deltaintegrate(y**2*DiracDelta(x - z)*DiracDelta(y - z), y)
|
||
|
|
z**2*DiracDelta(x - z)*Heaviside(y - z)
|
||
|
|
|
||
|
|
See Also
|
||
|
|
========
|
||
|
|
|
||
|
|
sympy.functions.special.delta_functions.DiracDelta
|
||
|
|
sympy.integrals.integrals.Integral
|
||
|
|
"""
|
||
|
|
if not f.has(DiracDelta):
|
||
|
|
return None
|
||
|
|
|
||
|
|
# g(x) = DiracDelta(h(x))
|
||
|
|
if f.func == DiracDelta:
|
||
|
|
h = f.expand(diracdelta=True, wrt=x)
|
||
|
|
if h == f: # can't simplify the expression
|
||
|
|
#FIXME: the second term tells whether is DeltaDirac or Derivative
|
||
|
|
#For integrating derivatives of DiracDelta we need the chain rule
|
||
|
|
if f.is_simple(x):
|
||
|
|
if (len(f.args) <= 1 or f.args[1] == 0):
|
||
|
|
return Heaviside(f.args[0])
|
||
|
|
else:
|
||
|
|
return (DiracDelta(f.args[0], f.args[1] - 1) /
|
||
|
|
f.args[0].as_poly().LC())
|
||
|
|
else: # let's try to integrate the simplified expression
|
||
|
|
fh = integrate(h, x)
|
||
|
|
return fh
|
||
|
|
elif f.is_Mul or f.is_Pow: # g(x) = a*b*c*f(DiracDelta(h(x)))*d*e
|
||
|
|
g = f.expand()
|
||
|
|
if f != g: # the expansion worked
|
||
|
|
fh = integrate(g, x)
|
||
|
|
if fh is not None and not isinstance(fh, Integral):
|
||
|
|
return fh
|
||
|
|
else:
|
||
|
|
# no expansion performed, try to extract a simple DiracDelta term
|
||
|
|
deltaterm, rest_mult = change_mul(f, x)
|
||
|
|
|
||
|
|
if not deltaterm:
|
||
|
|
if rest_mult:
|
||
|
|
fh = integrate(rest_mult, x)
|
||
|
|
return fh
|
||
|
|
else:
|
||
|
|
from sympy.solvers import solve
|
||
|
|
deltaterm = deltaterm.expand(diracdelta=True, wrt=x)
|
||
|
|
if deltaterm.is_Mul: # Take out any extracted factors
|
||
|
|
deltaterm, rest_mult_2 = change_mul(deltaterm, x)
|
||
|
|
rest_mult = rest_mult*rest_mult_2
|
||
|
|
point = solve(deltaterm.args[0], x)[0]
|
||
|
|
|
||
|
|
# Return the largest hyperreal term left after
|
||
|
|
# repeated integration by parts. For example,
|
||
|
|
#
|
||
|
|
# integrate(y*DiracDelta(x, 1),x) == y*DiracDelta(x,0), not 0
|
||
|
|
#
|
||
|
|
# This is so Integral(y*DiracDelta(x).diff(x),x).doit()
|
||
|
|
# will return y*DiracDelta(x) instead of 0 or DiracDelta(x),
|
||
|
|
# both of which are correct everywhere the value is defined
|
||
|
|
# but give wrong answers for nested integration.
|
||
|
|
n = (0 if len(deltaterm.args)==1 else deltaterm.args[1])
|
||
|
|
m = 0
|
||
|
|
while n >= 0:
|
||
|
|
r = S.NegativeOne**n*rest_mult.diff(x, n).subs(x, point)
|
||
|
|
if r.is_zero:
|
||
|
|
n -= 1
|
||
|
|
m += 1
|
||
|
|
else:
|
||
|
|
if m == 0:
|
||
|
|
return r*Heaviside(x - point)
|
||
|
|
else:
|
||
|
|
return r*DiracDelta(x,m-1)
|
||
|
|
# In some very weak sense, x=0 is still a singularity,
|
||
|
|
# but we hope will not be of any practical consequence.
|
||
|
|
return S.Zero
|
||
|
|
return None
|