Coverage for /usr/lib/python3/dist-packages/sympy/integrals/deltafunctions.py: 10%

1from sympy.core.mul import Mul 

2from sympy.core.singleton import S 

3from sympy.core.sorting import default_sort_key 

4from sympy.functions import DiracDelta, Heaviside 

5from .integrals import Integral, integrate 

6 

7 

8def change_mul(node, x): 

9 """change_mul(node, x) 

10 

11 Rearranges the operands of a product, bringing to front any simple 

12 DiracDelta expression. 

13 

14 Explanation 

15 =========== 

16 

17 If no simple DiracDelta expression was found, then all the DiracDelta 

18 expressions are simplified (using DiracDelta.expand(diracdelta=True, wrt=x)). 

19 

20 Return: (dirac, new node) 

21 Where: 

22 o dirac is either a simple DiracDelta expression or None (if no simple 

23 expression was found); 

24 o new node is either a simplified DiracDelta expressions or None (if it 

25 could not be simplified). 

26 

27 Examples 

28 ======== 

29 

30 >>> from sympy import DiracDelta, cos 

31 >>> from sympy.integrals.deltafunctions import change_mul 

32 >>> from sympy.abc import x, y 

33 >>> change_mul(x*y*DiracDelta(x)*cos(x), x) 

34 (DiracDelta(x), x*y*cos(x)) 

35 >>> change_mul(x*y*DiracDelta(x**2 - 1)*cos(x), x) 

36 (None, x*y*cos(x)*DiracDelta(x - 1)/2 + x*y*cos(x)*DiracDelta(x + 1)/2) 

37 >>> change_mul(x*y*DiracDelta(cos(x))*cos(x), x) 

38 (None, None) 

39 

40 See Also 

41 ======== 

42 

43 sympy.functions.special.delta_functions.DiracDelta 

44 deltaintegrate 

45 """ 

46 

47 new_args = [] 

48 dirac = None 

49 

50 #Sorting is needed so that we consistently collapse the same delta; 

51 #However, we must preserve the ordering of non-commutative terms 

52 c, nc = node.args_cnc() 

53 sorted_args = sorted(c, key=default_sort_key) 

54 sorted_args.extend(nc) 

55 

56 for arg in sorted_args: 

57 if arg.is_Pow and isinstance(arg.base, DiracDelta): 

58 new_args.append(arg.func(arg.base, arg.exp - 1)) 

59 arg = arg.base 

60 if dirac is None and (isinstance(arg, DiracDelta) and arg.is_simple(x)): 

61 dirac = arg 

62 else: 

63 new_args.append(arg) 

64 if not dirac: # there was no simple dirac 

65 new_args = [] 

66 for arg in sorted_args: 

67 if isinstance(arg, DiracDelta): 

68 new_args.append(arg.expand(diracdelta=True, wrt=x)) 

69 elif arg.is_Pow and isinstance(arg.base, DiracDelta): 

70 new_args.append(arg.func(arg.base.expand(diracdelta=True, wrt=x), arg.exp)) 

71 else: 

72 new_args.append(arg) 

73 if new_args != sorted_args: 

74 nnode = Mul(*new_args).expand() 

75 else: # if the node didn't change there is nothing to do 

76 nnode = None 

77 return (None, nnode) 

78 return (dirac, Mul(*new_args)) 

79 

80 

81def deltaintegrate(f, x): 

82 """ 

83 deltaintegrate(f, x) 

84 

85 Explanation 

86 =========== 

87 

88 The idea for integration is the following: 

89 

90 - If we are dealing with a DiracDelta expression, i.e. DiracDelta(g(x)), 

91 we try to simplify it. 

92 

93 If we could simplify it, then we integrate the resulting expression. 

94 We already know we can integrate a simplified expression, because only 

95 simple DiracDelta expressions are involved. 

96 

97 If we couldn't simplify it, there are two cases: 

98 

99 1) The expression is a simple expression: we return the integral, 

100 taking care if we are dealing with a Derivative or with a proper 

101 DiracDelta. 

102 

103 2) The expression is not simple (i.e. DiracDelta(cos(x))): we can do 

104 nothing at all. 

105 

106 - If the node is a multiplication node having a DiracDelta term: 

107 

108 First we expand it. 

109 

110 If the expansion did work, then we try to integrate the expansion. 

111 

112 If not, we try to extract a simple DiracDelta term, then we have two 

113 cases: 

114 

115 1) We have a simple DiracDelta term, so we return the integral. 

116 

117 2) We didn't have a simple term, but we do have an expression with 

118 simplified DiracDelta terms, so we integrate this expression. 

119 

120 Examples 

121 ======== 

122 

123 >>> from sympy.abc import x, y, z 

124 >>> from sympy.integrals.deltafunctions import deltaintegrate 

125 >>> from sympy import sin, cos, DiracDelta 

126 >>> deltaintegrate(x*sin(x)*cos(x)*DiracDelta(x - 1), x) 

127 sin(1)*cos(1)*Heaviside(x - 1) 

128 >>> deltaintegrate(y**2*DiracDelta(x - z)*DiracDelta(y - z), y) 

129 z**2*DiracDelta(x - z)*Heaviside(y - z) 

130 

131 See Also 

132 ======== 

133 

134 sympy.functions.special.delta_functions.DiracDelta 

135 sympy.integrals.integrals.Integral 

136 """ 

137 if not f.has(DiracDelta): 

138 return None 

139 

140 # g(x) = DiracDelta(h(x)) 

141 if f.func == DiracDelta: 

142 h = f.expand(diracdelta=True, wrt=x) 

143 if h == f: # can't simplify the expression 

144 #FIXME: the second term tells whether is DeltaDirac or Derivative 

145 #For integrating derivatives of DiracDelta we need the chain rule 

146 if f.is_simple(x): 

147 if (len(f.args) <= 1 or f.args[1] == 0): 

148 return Heaviside(f.args[0]) 

149 else: 

150 return (DiracDelta(f.args[0], f.args[1] - 1) / 

151 f.args[0].as_poly().LC()) 

152 else: # let's try to integrate the simplified expression 

153 fh = integrate(h, x) 

154 return fh 

155 elif f.is_Mul or f.is_Pow: # g(x) = a*b*c*f(DiracDelta(h(x)))*d*e 

156 g = f.expand() 

157 if f != g: # the expansion worked 

158 fh = integrate(g, x) 

159 if fh is not None and not isinstance(fh, Integral): 

160 return fh 

161 else: 

162 # no expansion performed, try to extract a simple DiracDelta term 

163 deltaterm, rest_mult = change_mul(f, x) 

164 

165 if not deltaterm: 

166 if rest_mult: 

167 fh = integrate(rest_mult, x) 

168 return fh 

169 else: 

170 from sympy.solvers import solve 

171 deltaterm = deltaterm.expand(diracdelta=True, wrt=x) 

172 if deltaterm.is_Mul: # Take out any extracted factors 

173 deltaterm, rest_mult_2 = change_mul(deltaterm, x) 

174 rest_mult = rest_mult*rest_mult_2 

175 point = solve(deltaterm.args[0], x)[0] 

176 

177 # Return the largest hyperreal term left after 

178 # repeated integration by parts. For example, 

179 # 

180 # integrate(y*DiracDelta(x, 1),x) == y*DiracDelta(x,0), not 0 

181 # 

182 # This is so Integral(y*DiracDelta(x).diff(x),x).doit() 

183 # will return y*DiracDelta(x) instead of 0 or DiracDelta(x), 

184 # both of which are correct everywhere the value is defined 

185 # but give wrong answers for nested integration. 

186 n = (0 if len(deltaterm.args)==1 else deltaterm.args[1]) 

187 m = 0 

188 while n >= 0: 

189 r = S.NegativeOne**n*rest_mult.diff(x, n).subs(x, point) 

190 if r.is_zero: 

191 n -= 1 

192 m += 1 

193 else: 

194 if m == 0: 

195 return r*Heaviside(x - point) 

196 else: 

197 return r*DiracDelta(x,m-1) 

198 # In some very weak sense, x=0 is still a singularity, 

199 # but we hope will not be of any practical consequence. 

200 return S.Zero 

201 return None