1from sympy.core import S, sympify, diff 2from sympy.core.decorators import deprecated 3from sympy.core.function import Function, ArgumentIndexError 4from sympy.core.logic import fuzzy_not 5from sympy.core.relational import Eq, Ne 6from sympy.functions.elementary.complexes import im, sign 7from sympy.functions.elementary.piecewise import Piecewise 8from sympy.polys.polyerrors import PolynomialError 9from sympy.utilities import filldedent 10 11 12############################################################################### 13################################ DELTA FUNCTION ############################### 14############################################################################### 15 16 17class DiracDelta(Function): 18 r""" 19 The DiracDelta function and its derivatives. 20 21 Explanation 22 =========== 23 24 DiracDelta is not an ordinary function. It can be rigorously defined either 25 as a distribution or as a measure. 26 27 DiracDelta only makes sense in definite integrals, and in particular, 28 integrals of the form ``Integral(f(x)*DiracDelta(x - x0), (x, a, b))``, 29 where it equals ``f(x0)`` if ``a <= x0 <= b`` and ``0`` otherwise. Formally, 30 DiracDelta acts in some ways like a function that is ``0`` everywhere except 31 at ``0``, but in many ways it also does not. It can often be useful to treat 32 DiracDelta in formal ways, building up and manipulating expressions with 33 delta functions (which may eventually be integrated), but care must be taken 34 to not treat it as a real function. SymPy's ``oo`` is similar. It only 35 truly makes sense formally in certain contexts (such as integration limits), 36 but SymPy allows its use everywhere, and it tries to be consistent with 37 operations on it (like ``1/oo``), but it is easy to get into trouble and get 38 wrong results if ``oo`` is treated too much like a number. Similarly, if 39 DiracDelta is treated too much like a function, it is easy to get wrong or 40 nonsensical results. 41 42 DiracDelta function has the following properties: 43 44 1) $\frac{d}{d x} \theta(x) = \delta(x)$ 45 2) $\int_{-\infty}^\infty \delta(x - a)f(x)\, dx = f(a)$ and $\int_{a- 46 \epsilon}^{a+\epsilon} \delta(x - a)f(x)\, dx = f(a)$ 47 3) $\delta(x) = 0$ for all $x \neq 0$ 48 4) $\delta(g(x)) = \sum_i \frac{\delta(x - x_i)}{\|g'(x_i)\|}$ where $x_i$ 49 are the roots of $g$ 50 5) $\delta(-x) = \delta(x)$ 51 52 Derivatives of ``k``-th order of DiracDelta have the following properties: 53 54 6) $\delta(x, k) = 0$ for all $x \neq 0$ 55 7) $\delta(-x, k) = -\delta(x, k)$ for odd $k$ 56 8) $\delta(-x, k) = \delta(x, k)$ for even $k$ 57 58 Examples 59 ======== 60 61 >>> from sympy import DiracDelta, diff, pi 62 >>> from sympy.abc import x, y 63 64 >>> DiracDelta(x) 65 DiracDelta(x) 66 >>> DiracDelta(1) 67 0 68 >>> DiracDelta(-1) 69 0 70 >>> DiracDelta(pi) 71 0 72 >>> DiracDelta(x - 4).subs(x, 4) 73 DiracDelta(0) 74 >>> diff(DiracDelta(x)) 75 DiracDelta(x, 1) 76 >>> diff(DiracDelta(x - 1),x,2) 77 DiracDelta(x - 1, 2) 78 >>> diff(DiracDelta(x**2 - 1),x,2) 79 2*(2*x**2*DiracDelta(x**2 - 1, 2) + DiracDelta(x**2 - 1, 1)) 80 >>> DiracDelta(3*x).is_simple(x) 81 True 82 >>> DiracDelta(x**2).is_simple(x) 83 False 84 >>> DiracDelta((x**2 - 1)*y).expand(diracdelta=True, wrt=x) 85 DiracDelta(x - 1)/(2*Abs(y)) + DiracDelta(x + 1)/(2*Abs(y)) 86 87 See Also 88 ======== 89 90 Heaviside 91 sympy.simplify.simplify.simplify, is_simple 92 sympy.functions.special.tensor_functions.KroneckerDelta 93 94 References 95 ========== 96 97 .. [1] http://mathworld.wolfram.com/DeltaFunction.html 98 99 """ 100 101 is_real = True 102 103 def fdiff(self, argindex=1): 104 """ 105 Returns the first derivative of a DiracDelta Function. 106 107 Explanation 108 =========== 109 110 The difference between ``diff()`` and ``fdiff()`` is: ``diff()`` is the 111 user-level function and ``fdiff()`` is an object method. ``fdiff()`` is 112 a convenience method available in the ``Function`` class. It returns 113 the derivative of the function without considering the chain rule. 114 ``diff(function, x)`` calls ``Function._eval_derivative`` which in turn 115 calls ``fdiff()`` internally to compute the derivative of the function. 116 117 Examples 118 ======== 119 120 >>> from sympy import DiracDelta, diff 121 >>> from sympy.abc import x 122 123 >>> DiracDelta(x).fdiff() 124 DiracDelta(x, 1) 125 126 >>> DiracDelta(x, 1).fdiff() 127 DiracDelta(x, 2) 128 129 >>> DiracDelta(x**2 - 1).fdiff() 130 DiracDelta(x**2 - 1, 1) 131 132 >>> diff(DiracDelta(x, 1)).fdiff() 133 DiracDelta(x, 3) 134 135 Parameters 136 ========== 137 138 argindex : integer 139 degree of derivative 140 141 """ 142 if argindex == 1: 143 #I didn't know if there is a better way to handle default arguments 144 k = 0 145 if len(self.args) > 1: 146 k = self.args[1] 147 return self.func(self.args[0], k + 1) 148 else: 149 raise ArgumentIndexError(self, argindex) 150 151 @classmethod 152 def eval(cls, arg, k=0): 153 """ 154 Returns a simplified form or a value of DiracDelta depending on the 155 argument passed by the DiracDelta object. 156 157 Explanation 158 =========== 159 160 The ``eval()`` method is automatically called when the ``DiracDelta`` 161 class is about to be instantiated and it returns either some simplified 162 instance or the unevaluated instance depending on the argument passed. 163 In other words, ``eval()`` method is not needed to be called explicitly, 164 it is being called and evaluated once the object is called. 165 166 Examples 167 ======== 168 169 >>> from sympy import DiracDelta, S 170 >>> from sympy.abc import x 171 172 >>> DiracDelta(x) 173 DiracDelta(x) 174 175 >>> DiracDelta(-x, 1) 176 -DiracDelta(x, 1) 177 178 >>> DiracDelta(1) 179 0 180 181 >>> DiracDelta(5, 1) 182 0 183 184 >>> DiracDelta(0) 185 DiracDelta(0) 186 187 >>> DiracDelta(-1) 188 0 189 190 >>> DiracDelta(S.NaN) 191 nan 192 193 >>> DiracDelta(x).eval(1) 194 0 195 196 >>> DiracDelta(x - 100).subs(x, 5) 197 0 198 199 >>> DiracDelta(x - 100).subs(x, 100) 200 DiracDelta(0) 201 202 Parameters 203 ========== 204 205 k : integer 206 order of derivative 207 208 arg : argument passed to DiracDelta 209 210 """ 211 k = sympify(k) 212 if not k.is_Integer or k.is_negative: 213 raise ValueError("Error: the second argument of DiracDelta must be \ 214 a non-negative integer, %s given instead." % (k,)) 215 arg = sympify(arg) 216 if arg is S.NaN: 217 return S.NaN 218 if arg.is_nonzero: 219 return S.Zero 220 if fuzzy_not(im(arg).is_zero): 221 raise ValueError(filldedent(''' 222 Function defined only for Real Values. 223 Complex part: %s found in %s .''' % ( 224 repr(im(arg)), repr(arg)))) 225 c, nc = arg.args_cnc() 226 if c and c[0] is S.NegativeOne: 227 # keep this fast and simple instead of using 228 # could_extract_minus_sign 229 if k.is_odd: 230 return -cls(-arg, k) 231 elif k.is_even: 232 return cls(-arg, k) if k else cls(-arg) 233 234 @deprecated(useinstead="expand(diracdelta=True, wrt=x)", issue=12859, deprecated_since_version="1.1") 235 def simplify(self, x, **kwargs): 236 return self.expand(diracdelta=True, wrt=x) 237 238 def _eval_expand_diracdelta(self, **hints): 239 """ 240 Compute a simplified representation of the function using 241 property number 4. Pass ``wrt`` as a hint to expand the expression 242 with respect to a particular variable. 243 244 Explanation 245 =========== 246 247 ``wrt`` is: 248 249 - a variable with respect to which a DiracDelta expression will 250 get expanded. 251 252 Examples 253 ======== 254 255 >>> from sympy import DiracDelta 256 >>> from sympy.abc import x, y 257 258 >>> DiracDelta(x*y).expand(diracdelta=True, wrt=x) 259 DiracDelta(x)/Abs(y) 260 >>> DiracDelta(x*y).expand(diracdelta=True, wrt=y) 261 DiracDelta(y)/Abs(x) 262 263 >>> DiracDelta(x**2 + x - 2).expand(diracdelta=True, wrt=x) 264 DiracDelta(x - 1)/3 + DiracDelta(x + 2)/3 265 266 See Also 267 ======== 268 269 is_simple, Diracdelta 270 271 """ 272 from sympy.polys.polyroots import roots 273 274 wrt = hints.get('wrt', None) 275 if wrt is None: 276 free = self.free_symbols 277 if len(free) == 1: 278 wrt = free.pop() 279 else: 280 raise TypeError(filldedent(''' 281 When there is more than 1 free symbol or variable in the expression, 282 the 'wrt' keyword is required as a hint to expand when using the 283 DiracDelta hint.''')) 284 285 if not self.args[0].has(wrt) or (len(self.args) > 1 and self.args[1] != 0 ): 286 return self 287 try: 288 argroots = roots(self.args[0], wrt) 289 result = 0 290 valid = True 291 darg = abs(diff(self.args[0], wrt)) 292 for r, m in argroots.items(): 293 if r.is_real is not False and m == 1: 294 result += self.func(wrt - r)/darg.subs(wrt, r) 295 else: 296 # don't handle non-real and if m != 1 then 297 # a polynomial will have a zero in the derivative (darg) 298 # at r 299 valid = False 300 break 301 if valid: 302 return result 303 except PolynomialError: 304 pass 305 return self 306 307 def is_simple(self, x): 308 """ 309 Tells whether the argument(args[0]) of DiracDelta is a linear 310 expression in *x*. 311 312 Examples 313 ======== 314 315 >>> from sympy import DiracDelta, cos 316 >>> from sympy.abc import x, y 317 318 >>> DiracDelta(x*y).is_simple(x) 319 True 320 >>> DiracDelta(x*y).is_simple(y) 321 True 322 323 >>> DiracDelta(x**2 + x - 2).is_simple(x) 324 False 325 326 >>> DiracDelta(cos(x)).is_simple(x) 327 False 328 329 Parameters 330 ========== 331 332 x : can be a symbol 333 334 See Also 335 ======== 336 337 sympy.simplify.simplify.simplify, DiracDelta 338 339 """ 340 p = self.args[0].as_poly(x) 341 if p: 342 return p.degree() == 1 343 return False 344 345 def _eval_rewrite_as_Piecewise(self, *args, **kwargs): 346 """ 347 Represents DiracDelta in a piecewise form. 348 349 Examples 350 ======== 351 352 >>> from sympy import DiracDelta, Piecewise, Symbol 353 >>> x = Symbol('x') 354 355 >>> DiracDelta(x).rewrite(Piecewise) 356 Piecewise((DiracDelta(0), Eq(x, 0)), (0, True)) 357 358 >>> DiracDelta(x - 5).rewrite(Piecewise) 359 Piecewise((DiracDelta(0), Eq(x - 5, 0)), (0, True)) 360 361 >>> DiracDelta(x**2 - 5).rewrite(Piecewise) 362 Piecewise((DiracDelta(0), Eq(x**2 - 5, 0)), (0, True)) 363 364 >>> DiracDelta(x - 5, 4).rewrite(Piecewise) 365 DiracDelta(x - 5, 4) 366 367 """ 368 if len(args) == 1: 369 return Piecewise((DiracDelta(0), Eq(args[0], 0)), (0, True)) 370 371 def _eval_rewrite_as_SingularityFunction(self, *args, **kwargs): 372 """ 373 Returns the DiracDelta expression written in the form of Singularity 374 Functions. 375 376 """ 377 from sympy.solvers import solve 378 from sympy.functions import SingularityFunction 379 if self == DiracDelta(0): 380 return SingularityFunction(0, 0, -1) 381 if self == DiracDelta(0, 1): 382 return SingularityFunction(0, 0, -2) 383 free = self.free_symbols 384 if len(free) == 1: 385 x = (free.pop()) 386 if len(args) == 1: 387 return SingularityFunction(x, solve(args[0], x)[0], -1) 388 return SingularityFunction(x, solve(args[0], x)[0], -args[1] - 1) 389 else: 390 # I don't know how to handle the case for DiracDelta expressions 391 # having arguments with more than one variable. 392 raise TypeError(filldedent(''' 393 rewrite(SingularityFunction) doesn't support 394 arguments with more that 1 variable.''')) 395 396 397############################################################################### 398############################## HEAVISIDE FUNCTION ############################# 399############################################################################### 400 401 402class Heaviside(Function): 403 r""" 404 Heaviside step function. 405 406 Explanation 407 =========== 408 409 The Heaviside step function has the following properties: 410 411 1) $\frac{d}{d x} \theta(x) = \delta(x)$ 412 2) $\theta(x) = \begin{cases} 0 & \text{for}\: x < 0 \\ \frac{1}{2} & 413 \text{for}\: x = 0 \\1 & \text{for}\: x > 0 \end{cases}$ 414 3) $\frac{d}{d x} \max(x, 0) = \theta(x)$ 415 416 Heaviside(x) is printed as $\theta(x)$ with the SymPy LaTeX printer. 417 418 The value at 0 is set differently in different fields. SymPy uses 1/2, 419 which is a convention from electronics and signal processing, and is 420 consistent with solving improper integrals by Fourier transform and 421 convolution. 422 423 To specify a different value of Heaviside at ``x=0``, a second argument 424 can be given. Using ``Heaviside(x, nan)`` gives an expression that will 425 evaluate to nan for x=0. 426 427 .. versionchanged:: 1.9 ``Heaviside(0)`` now returns 1/2 (before: undefined) 428 429 Examples 430 ======== 431 432 >>> from sympy import Heaviside, nan 433 >>> from sympy.abc import x 434 >>> Heaviside(9) 435 1 436 >>> Heaviside(-9) 437 0 438 >>> Heaviside(0) 439 1/2 440 >>> Heaviside(0, nan) 441 nan 442 >>> (Heaviside(x) + 1).replace(Heaviside(x), Heaviside(x, 1)) 443 Heaviside(x, 1) + 1 444 445 See Also 446 ======== 447 448 DiracDelta 449 450 References 451 ========== 452 453 .. [1] http://mathworld.wolfram.com/HeavisideStepFunction.html 454 .. [2] http://dlmf.nist.gov/1.16#iv 455 456 """ 457 458 is_real = True 459 460 def fdiff(self, argindex=1): 461 """ 462 Returns the first derivative of a Heaviside Function. 463 464 Examples 465 ======== 466 467 >>> from sympy import Heaviside, diff 468 >>> from sympy.abc import x 469 470 >>> Heaviside(x).fdiff() 471 DiracDelta(x) 472 473 >>> Heaviside(x**2 - 1).fdiff() 474 DiracDelta(x**2 - 1) 475 476 >>> diff(Heaviside(x)).fdiff() 477 DiracDelta(x, 1) 478 479 Parameters 480 ========== 481 482 argindex : integer 483 order of derivative 484 485 """ 486 if argindex == 1: 487 return DiracDelta(self.args[0]) 488 else: 489 raise ArgumentIndexError(self, argindex) 490 491 def __new__(cls, arg, H0=S.Half, **options): 492 if isinstance(H0, Heaviside) and len(H0.args) == 1: 493 H0 = S.Half 494 return super(cls, cls).__new__(cls, arg, H0, **options) 495 496 @property 497 def pargs(self): 498 """Args without default S.Half""" 499 args = self.args 500 if args[1] is S.Half: 501 args = args[:1] 502 return args 503 504 @classmethod 505 def eval(cls, arg, H0=S.Half): 506 """ 507 Returns a simplified form or a value of Heaviside depending on the 508 argument passed by the Heaviside object. 509 510 Explanation 511 =========== 512 513 The ``eval()`` method is automatically called when the ``Heaviside`` 514 class is about to be instantiated and it returns either some simplified 515 instance or the unevaluated instance depending on the argument passed. 516 In other words, ``eval()`` method is not needed to be called explicitly, 517 it is being called and evaluated once the object is called. 518 519 Examples 520 ======== 521 522 >>> from sympy import Heaviside, S 523 >>> from sympy.abc import x 524 525 >>> Heaviside(x) 526 Heaviside(x) 527 528 >>> Heaviside(19) 529 1 530 531 >>> Heaviside(0) 532 1/2 533 534 >>> Heaviside(0, 1) 535 1 536 537 >>> Heaviside(-5) 538 0 539 540 >>> Heaviside(S.NaN) 541 nan 542 543 >>> Heaviside(x).eval(42) 544 1 545 546 >>> Heaviside(x - 100).subs(x, 5) 547 0 548 549 >>> Heaviside(x - 100).subs(x, 105) 550 1 551 552 Parameters 553 ========== 554 555 arg : argument passed by Heaviside object 556 557 H0 : value of Heaviside(0) 558 559 """ 560 H0 = sympify(H0) 561 arg = sympify(arg) 562 if arg.is_extended_negative: 563 return S.Zero 564 elif arg.is_extended_positive: 565 return S.One 566 elif arg.is_zero: 567 return H0 568 elif arg is S.NaN: 569 return S.NaN 570 elif fuzzy_not(im(arg).is_zero): 571 raise ValueError("Function defined only for Real Values. Complex part: %s found in %s ." % (repr(im(arg)), repr(arg)) ) 572 573 def _eval_rewrite_as_Piecewise(self, arg, H0=None, **kwargs): 574 """ 575 Represents Heaviside in a Piecewise form. 576 577 Examples 578 ======== 579 580 >>> from sympy import Heaviside, Piecewise, Symbol, nan 581 >>> x = Symbol('x') 582 583 >>> Heaviside(x).rewrite(Piecewise) 584 Piecewise((0, x < 0), (1/2, Eq(x, 0)), (1, x > 0)) 585 586 >>> Heaviside(x,nan).rewrite(Piecewise) 587 Piecewise((0, x < 0), (nan, Eq(x, 0)), (1, x > 0)) 588 589 >>> Heaviside(x - 5).rewrite(Piecewise) 590 Piecewise((0, x - 5 < 0), (1/2, Eq(x - 5, 0)), (1, x - 5 > 0)) 591 592 >>> Heaviside(x**2 - 1).rewrite(Piecewise) 593 Piecewise((0, x**2 - 1 < 0), (1/2, Eq(x**2 - 1, 0)), (1, x**2 - 1 > 0)) 594 595 """ 596 if H0 == 0: 597 return Piecewise((0, arg <= 0), (1, arg > 0)) 598 if H0 == 1: 599 return Piecewise((0, arg < 0), (1, arg >= 0)) 600 return Piecewise((0, arg < 0), (H0, Eq(arg, 0)), (1, arg > 0)) 601 602 def _eval_rewrite_as_sign(self, arg, H0=S.Half, **kwargs): 603 """ 604 Represents the Heaviside function in the form of sign function. 605 606 Explanation 607 =========== 608 609 The value of Heaviside(0) must be 1/2 for rewritting as sign to be 610 strictly equivalent. For easier usage, we also allow this rewriting 611 when Heaviside(0) is undefined. 612 613 Examples 614 ======== 615 616 >>> from sympy import Heaviside, Symbol, sign, nan 617 >>> x = Symbol('x', real=True) 618 >>> y = Symbol('y') 619 620 >>> Heaviside(x).rewrite(sign) 621 sign(x)/2 + 1/2 622 623 >>> Heaviside(x, 0).rewrite(sign) 624 Piecewise((sign(x)/2 + 1/2, Ne(x, 0)), (0, True)) 625 626 >>> Heaviside(x, nan).rewrite(sign) 627 Piecewise((sign(x)/2 + 1/2, Ne(x, 0)), (nan, True)) 628 629 >>> Heaviside(x - 2).rewrite(sign) 630 sign(x - 2)/2 + 1/2 631 632 >>> Heaviside(x**2 - 2*x + 1).rewrite(sign) 633 sign(x**2 - 2*x + 1)/2 + 1/2 634 635 >>> Heaviside(y).rewrite(sign) 636 Heaviside(y) 637 638 >>> Heaviside(y**2 - 2*y + 1).rewrite(sign) 639 Heaviside(y**2 - 2*y + 1) 640 641 See Also 642 ======== 643 644 sign 645 646 """ 647 if arg.is_extended_real: 648 pw1 = Piecewise( 649 ((sign(arg) + 1)/2, Ne(arg, 0)), 650 (Heaviside(0, H0=H0), True)) 651 pw2 = Piecewise( 652 ((sign(arg) + 1)/2, Eq(Heaviside(0, H0=H0), S(1)/2)), 653 (pw1, True)) 654 return pw2 655 656 def _eval_rewrite_as_SingularityFunction(self, args, H0=S.Half, **kwargs): 657 """ 658 Returns the Heaviside expression written in the form of Singularity 659 Functions. 660 661 """ 662 from sympy.solvers import solve 663 from sympy.functions import SingularityFunction 664 if self == Heaviside(0): 665 return SingularityFunction(0, 0, 0) 666 free = self.free_symbols 667 if len(free) == 1: 668 x = (free.pop()) 669 return SingularityFunction(x, solve(args, x)[0], 0) 670 # TODO 671 # ((x - 5)**3*Heaviside(x - 5)).rewrite(SingularityFunction) should output 672 # SingularityFunction(x, 5, 0) instead of (x - 5)**3*SingularityFunction(x, 5, 0) 673 else: 674 # I don't know how to handle the case for Heaviside expressions 675 # having arguments with more than one variable. 676 raise TypeError(filldedent(''' 677 rewrite(SingularityFunction) doesn't 678 support arguments with more that 1 variable.''')) 679