Repository: sympy/sympy
Pull request: Modified Heaviside(t) to have a standard value at Heaviside(0)=1/2
#### References to other Issues or PRs

related to https://github.com/sympy/sympy/pull/20411
fixes #21392

#### Brief description of what is fixed or changed

Modified Heaviside(t) to have a standard value at Heaviside(0)=1/2 with all tests adapted.

The test_issue_15923() are removed, they are not necessary anymore.

See discussion in https://github.com/sympy/sympy/pull/20411 for more details.

#### Other comments

#### Release Notes

<!-- BEGIN RELEASE NOTES -->

* functions
  * Heaviside(0) now returns 1/2 instead of being undefined.

<!-- END RELEASE NOTES -->


Problem statement:
get_points() generates error when used with Heaviside() function
The following code used to work with older versions of sympy (e.g. 1.3) but generates an error with 1.7.1

```
import sympy as sym

t = sym.symbols('t')
p = sym.plot(sym.Heaviside(t), (t, -5, 5), show=False)
points = p[0].get_points()
```

results in

`AttributeError: 'ImmutableDenseNDimArray' object has no attribute 'as_coefficient'`

Hints from issue discussion:

Commit message: Post-review changes:
* Heaviside(0,nan) now evaluates and returns nan
* Tests tidied up
* Documentation improved

Patch to review:
diff --git a/sympy/functions/elementary/complexes.py b/sympy/functions/elementary/complexes.py
index f4cdfc9f9791..3c3e8d21b3cf 100644
--- a/sympy/functions/elementary/complexes.py
+++ b/sympy/functions/elementary/complexes.py
@@ -442,7 +442,7 @@ def _eval_rewrite_as_Piecewise(self, arg, **kwargs):
     def _eval_rewrite_as_Heaviside(self, arg, **kwargs):
         from sympy.functions.special.delta_functions import Heaviside
         if arg.is_extended_real:
-            return Heaviside(arg, H0=S(1)/2) * 2 - 1
+            return Heaviside(arg) * 2 - 1
 
     def _eval_rewrite_as_Abs(self, arg, **kwargs):
         return Piecewise((0, Eq(arg, 0)), (arg / Abs(arg), True))
diff --git a/sympy/functions/elementary/tests/test_trigonometric.py b/sympy/functions/elementary/tests/test_trigonometric.py
index 65b068349b09..61b63a319899 100644
--- a/sympy/functions/elementary/tests/test_trigonometric.py
+++ b/sympy/functions/elementary/tests/test_trigonometric.py
@@ -1049,7 +1049,7 @@ def test_atan2():
     assert atan2(0, u) == pi
 
     assert atan2(y, oo) ==  0
-    assert atan2(y, -oo)==  2*pi*Heaviside(re(y)) - pi
+    assert atan2(y, -oo)==  2*pi*Heaviside(re(y), S.Half) - pi
 
     assert atan2(y, x).rewrite(log) == -I*log((x + I*y)/sqrt(x**2 + y**2))
     assert atan2(0, 0) is S.NaN
diff --git a/sympy/functions/special/delta_functions.py b/sympy/functions/special/delta_functions.py
index e5315e181d51..480f76fe525f 100644
--- a/sympy/functions/special/delta_functions.py
+++ b/sympy/functions/special/delta_functions.py
@@ -404,8 +404,11 @@ def _sage_(self):
 
 
 class Heaviside(Function):
-    r"""
-    Heaviside Piecewise function.
+    r"""Heaviside step function.
+
+    .. versionchanged:: 1.9
+        ``Heaviside(0)`` now returns 1/2; previously, the default was that
+        the value at 0 was undefined.
 
     Explanation
     ===========
@@ -413,35 +416,34 @@ class Heaviside(Function):
     Heaviside function has the following properties:
 
     1) $\frac{d}{d x} \theta(x) = \delta(x)$
-    2) $\theta(x) = \begin{cases} 0 & \text{for}\: x < 0 \\ \text{undefined} &
+    2) $\theta(x) = \begin{cases} 0 & \text{for}\: x < 0 \\ \frac{1}{2} &
        \text{for}\: x = 0 \\1 & \text{for}\: x > 0 \end{cases}$
     3) $\frac{d}{d x} \max(x, 0) = \theta(x)$
 
     Heaviside(x) is printed as $\theta(x)$ with the SymPy LaTeX printer.
 
-    Regarding to the value at 0, Mathematica defines $\theta(0)=1$, but Maple
-    uses $\theta(0) = \text{undefined}$. Different application areas may have
-    specific conventions. For example, in control theory, it is common practice
-    to assume $\theta(0) = 0$ to match the Laplace transform of a DiracDelta
-    distribution.
+    The value at 0 is set differently in different fields. SymPy uses 1/2,
+    which is a convention from electronics and signal processing, and is
+    consistent with solving improper integrals by Fourier transform and
+    convolution.
 
-    To specify the value of Heaviside at ``x=0``, a second argument can be
-    given. Omit this 2nd argument or pass ``None`` to recover the default
-    behavior.
+    To specify a different value of Heaviside at ``x=0``, a second argument
+    can be given. Using ``Heaviside(x, nan)`` gives an expression that will
+    evaluate to nan for x=0.
 
     Examples
     ========
 
-    >>> from sympy import Heaviside, S
+    >>> from sympy import Heaviside, nan
     >>> from sympy.abc import x
     >>> Heaviside(9)
     1
     >>> Heaviside(-9)
     0
     >>> Heaviside(0)
-    Heaviside(0)
-    >>> Heaviside(0, S.Half)
     1/2
+    >>> Heaviside(0, nan)
+    nan
     >>> (Heaviside(x) + 1).replace(Heaviside(x), Heaviside(x, 1))
     Heaviside(x, 1) + 1
 
@@ -487,21 +489,17 @@ def fdiff(self, argindex=1):
 
         """
         if argindex == 1:
-            # property number 1
             return DiracDelta(self.args[0])
         else:
             raise ArgumentIndexError(self, argindex)
 
-    def __new__(cls, arg, H0=None, **options):
+    def __new__(cls, arg, H0=S.Half, **options):
         if isinstance(H0, Heaviside) and len(H0.args) == 1:
-            H0 = None
-
-        if H0 is None:
-            return super(cls, cls).__new__(cls, arg, **options)
+            H0 = S.Half
         return super(cls, cls).__new__(cls, arg, H0, **options)
 
     @classmethod
-    def eval(cls, arg, H0=None):
+    def eval(cls, arg, H0=S.Half):
         """
         Returns a simplified form or a value of Heaviside depending on the
         argument passed by the Heaviside object.
@@ -528,7 +526,7 @@ class is about to be instantiated and it returns either some simplified
         1
 
         >>> Heaviside(0)
-        Heaviside(0)
+        1/2
 
         >>> Heaviside(0, 1)
         1
@@ -539,7 +537,7 @@ class is about to be instantiated and it returns either some simplified
         >>> Heaviside(S.NaN)
         nan
 
-        >>> Heaviside(x).eval(100)
+        >>> Heaviside(x).eval(42)
         1
 
         >>> Heaviside(x - 100).subs(x, 5)
@@ -551,8 +549,9 @@ class is about to be instantiated and it returns either some simplified
         Parameters
         ==========
 
-        arg : argument passed by HeaviSide object
-        HO : boolean flag for HeaviSide Object is set to True
+        arg : argument passed by Heaviside object
+
+        H0 : value of Heaviside(0)
 
         """
         H0 = sympify(H0)
@@ -575,58 +574,61 @@ def _eval_rewrite_as_Piecewise(self, arg, H0=None, **kwargs):
         Examples
         ========
 
-        >>> from sympy import Heaviside, Piecewise, Symbol
+        >>> from sympy import Heaviside, Piecewise, Symbol, nan
         >>> x = Symbol('x')
 
         >>> Heaviside(x).rewrite(Piecewise)
-        Piecewise((0, x < 0), (Heaviside(0), Eq(x, 0)), (1, x > 0))
+        Piecewise((0, x < 0), (1/2, Eq(x, 0)), (1, x > 0))
+
+        >>> Heaviside(x,nan).rewrite(Piecewise)
+        Piecewise((0, x < 0), (nan, Eq(x, 0)), (1, x > 0))
 
         >>> Heaviside(x - 5).rewrite(Piecewise)
-        Piecewise((0, x - 5 < 0), (Heaviside(0), Eq(x - 5, 0)), (1, x - 5 > 0))
+        Piecewise((0, x - 5 < 0), (1/2, Eq(x - 5, 0)), (1, x - 5 > 0))
 
         >>> Heaviside(x**2 - 1).rewrite(Piecewise)
-        Piecewise((0, x**2 - 1 < 0), (Heaviside(0), Eq(x**2 - 1, 0)), (1, x**2 - 1 > 0))
+        Piecewise((0, x**2 - 1 < 0), (1/2, Eq(x**2 - 1, 0)), (1, x**2 - 1 > 0))
 
         """
-        if H0 is None:
-            return Piecewise((0, arg < 0), (Heaviside(0), Eq(arg, 0)), (1, arg > 0))
         if H0 == 0:
             return Piecewise((0, arg <= 0), (1, arg > 0))
         if H0 == 1:
             return Piecewise((0, arg < 0), (1, arg >= 0))
         return Piecewise((0, arg < 0), (H0, Eq(arg, 0)), (1, arg > 0))
 
-    def _eval_rewrite_as_sign(self, arg, H0=None, **kwargs):
+    def _eval_rewrite_as_sign(self, arg, H0=S.Half, **kwargs):
         """
         Represents the Heaviside function in the form of sign function.
 
         Explanation
         ===========
 
-        The value of the second argument of Heaviside must specify Heaviside(0)
-        = 1/2 for rewritting as sign to be strictly equivalent. For easier
-        usage, we also allow this rewriting when Heaviside(0) is undefined.
+        The value of Heaviside(0) must be 1/2 for rewritting as sign to be
+        strictly equivalent. For easier usage, we also allow this rewriting
+        when Heaviside(0) is undefined.
 
         Examples
         ========
 
-        >>> from sympy import Heaviside, Symbol, sign, S
+        >>> from sympy import Heaviside, Symbol, sign, nan
         >>> x = Symbol('x', real=True)
+        >>> y = Symbol('y')
 
-        >>> Heaviside(x, H0=S.Half).rewrite(sign)
+        >>> Heaviside(x).rewrite(sign)
         sign(x)/2 + 1/2
 
         >>> Heaviside(x, 0).rewrite(sign)
         Piecewise((sign(x)/2 + 1/2, Ne(x, 0)), (0, True))
 
-        >>> Heaviside(x - 2, H0=S.Half).rewrite(sign)
+        >>> Heaviside(x, nan).rewrite(sign)
+        Piecewise((sign(x)/2 + 1/2, Ne(x, 0)), (nan, True))
+
+        >>> Heaviside(x - 2).rewrite(sign)
         sign(x - 2)/2 + 1/2
 
-        >>> Heaviside(x**2 - 2*x + 1, H0=S.Half).rewrite(sign)
+        >>> Heaviside(x**2 - 2*x + 1).rewrite(sign)
         sign(x**2 - 2*x + 1)/2 + 1/2
 
-        >>> y = Symbol('y')
-
         >>> Heaviside(y).rewrite(sign)
         Heaviside(y)
 
@@ -648,7 +650,7 @@ def _eval_rewrite_as_sign(self, arg, H0=None, **kwargs):
                 (pw1, True))
             return pw2
 
-    def _eval_rewrite_as_SingularityFunction(self, args, **kwargs):
+    def _eval_rewrite_as_SingularityFunction(self, args, H0=S.Half, **kwargs):
         """
         Returns the Heaviside expression written in the form of Singularity
         Functions.
diff --git a/sympy/functions/special/tests/test_delta_functions.py b/sympy/functions/special/tests/test_delta_functions.py
index 72ba93e2c96e..69722eac1bc8 100644
--- a/sympy/functions/special/tests/test_delta_functions.py
+++ b/sympy/functions/special/tests/test_delta_functions.py
@@ -6,6 +6,8 @@
 
 from sympy.testing.pytest import raises, warns_deprecated_sympy
 
+from sympy.core.expr import unchanged
+
 from sympy.core.function import ArgumentIndexError
 
 
@@ -80,19 +82,13 @@ def test_DiracDelta():
 
 
 def test_heaviside():
-    assert Heaviside(0).func == Heaviside
     assert Heaviside(-5) == 0
     assert Heaviside(1) == 1
-    assert Heaviside(nan) is nan
+    assert Heaviside(0) == S.Half
 
     assert Heaviside(0, x) == x
-    assert Heaviside(0, nan) is nan
-    assert Heaviside(x, None) == Heaviside(x)
-    assert Heaviside(0, None) == Heaviside(0)
-
-    # we do not want None and Heaviside(0) in the args:
-    assert Heaviside(x, H0=None).args == (x,)
-    assert Heaviside(x, H0=Heaviside(0)).args == (x,)
+    assert unchanged(Heaviside,x, nan)
+    assert Heaviside(0, nan) == nan
 
     assert adjoint(Heaviside(x)) == Heaviside(x)
     assert adjoint(Heaviside(x - y)) == Heaviside(x - y)
@@ -122,6 +118,8 @@ def test_rewrite():
         Piecewise((0, x <= 0), (1, x > 0)))
     assert Heaviside(x, 1).rewrite(Piecewise) == (
         Piecewise((0, x < 0), (1, x >= 0)))
+    assert Heaviside(x, nan).rewrite(Piecewise) == (
+        Piecewise((0, x < 0), (nan, Eq(x, 0)), (1, x > 0)))
 
     assert Heaviside(x).rewrite(sign) == \
         Heaviside(x, H0=Heaviside(0)).rewrite(sign) == \
@@ -153,13 +151,4 @@ def test_rewrite():
     assert Heaviside(x).rewrite(SingularityFunction) == SingularityFunction(x, 0, 0)
     assert 5*x*y*Heaviside(y + 1).rewrite(SingularityFunction) == 5*x*y*SingularityFunction(y, -1, 0)
     assert ((x - 3)**3*Heaviside(x - 3)).rewrite(SingularityFunction) == (x - 3)**3*SingularityFunction(x, 3, 0)
-    assert Heaviside(0).rewrite(SingularityFunction) == SingularityFunction(0, 0, 0)
-
-def test_issue_15923():
-    x = Symbol('x', real=True)
-    assert Heaviside(x).rewrite(Piecewise, H0=0) == (
-        Piecewise((0, x <= 0), (1, True)))
-    assert Heaviside(x).rewrite(Piecewise, H0=1) == (
-        Piecewise((0, x < 0), (1, True)))
-    assert Heaviside(x).rewrite(Piecewise, H0=S.Half) == (
-        Piecewise((0, x < 0), (S.Half, Eq(x, 0)), (1, x > 0)))
+    assert Heaviside(0).rewrite(SingularityFunction) == S.Half
diff --git a/sympy/integrals/transforms.py b/sympy/integrals/transforms.py
index 739450dd9391..4f5218052bab 100644
--- a/sympy/integrals/transforms.py
+++ b/sympy/integrals/transforms.py
@@ -1234,17 +1234,17 @@ def pw_simp(*args):
 
     u = Dummy('u')
 
-    def simp_heaviside(arg):
+    def simp_heaviside(arg, H0=S.Half):
         a = arg.subs(exp(-t), u)
         if a.has(t):
-            return Heaviside(arg)
+            return Heaviside(arg,H0)
         rel = _solve_inequality(a > 0, u)
         if rel.lts == u:
             k = log(rel.gts)
-            return Heaviside(t + k)
+            return Heaviside(t + k,H0)
         else:
             k = log(rel.lts)
-            return Heaviside(-(t + k))
+            return Heaviside(-(t + k),H0)
     f = f.replace(Heaviside, simp_heaviside)
 
     def simp_exp(arg):
diff --git a/sympy/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py
index 42fae775d346..eaede3c00c44 100644
--- a/sympy/printing/pretty/tests/test_pretty.py
+++ b/sympy/printing/pretty/tests/test_pretty.py
@@ -7224,8 +7224,8 @@ def test_pretty_misc_functions():
     assert upretty(fresnelc(x)) == 'C(x)'
     assert pretty(fresnels(x)) == 'S(x)'
     assert upretty(fresnels(x)) == 'S(x)'
-    assert pretty(Heaviside(x)) == 'Heaviside(x)'
-    assert upretty(Heaviside(x)) == 'θ(x)'
+    assert pretty(Heaviside(x)) == 'Heaviside(x, 1/2)'
+    assert upretty(Heaviside(x)) == 'θ(x, 1/2)'
     assert pretty(Heaviside(x, y)) == 'Heaviside(x, y)'
     assert upretty(Heaviside(x, y)) == 'θ(x, y)'
     assert pretty(dirichlet_eta(x)) == 'dirichlet_eta(x)'
diff --git a/sympy/printing/tests/test_mathml.py b/sympy/printing/tests/test_mathml.py
index ffd59826093d..a2404718a911 100644
--- a/sympy/printing/tests/test_mathml.py
+++ b/sympy/printing/tests/test_mathml.py
@@ -1530,7 +1530,7 @@ def test_print_different_functions():
     assert mpp.doprint(primeomega(x)) == '<mrow><mi>&#x3A9;</mi><mfenced><mi>x</mi></mfenced></mrow>'
     assert mpp.doprint(fresnels(x)) == '<mrow><mi>S</mi><mfenced><mi>x</mi></mfenced></mrow>'
     assert mpp.doprint(fresnelc(x)) ==  '<mrow><mi>C</mi><mfenced><mi>x</mi></mfenced></mrow>'
-    assert mpp.doprint(Heaviside(x)) == '<mrow><mi>&#x398;</mi><mfenced><mi>x</mi></mfenced></mrow>'
+    assert mpp.doprint(Heaviside(x)) == '<mrow><mi>&#x398;</mi><mfenced><mi>x</mi><mfrac><mn>1</mn><mn>2</mn></mfrac></mfenced></mrow>'
 
 
 def test_mathml_builtins():
diff --git a/sympy/printing/tests/test_octave.py b/sympy/printing/tests/test_octave.py
index bdd0dca02e12..b39116d4606d 100644
--- a/sympy/printing/tests/test_octave.py
+++ b/sympy/printing/tests/test_octave.py
@@ -90,7 +90,7 @@ def test_Function_change_name():
     assert mcode(RisingFactorial(x, y)) == "pochhammer(x, y)"
     assert mcode(DiracDelta(x)) == "dirac(x)"
     assert mcode(DiracDelta(x, 3)) == "dirac(3, x)"
-    assert mcode(Heaviside(x)) == "heaviside(x)"
+    assert mcode(Heaviside(x)) == "heaviside(x, 1/2)"
     assert mcode(Heaviside(x, y)) == "heaviside(x, y)"
     assert mcode(binomial(x, y)) == "bincoeff(x, y)"
     assert mcode(Mod(x, y)) == "mod(x, y)"


Code review comments:
diff_hunk: @@ -7224,8 +7224,8 @@ def test_pretty_misc_functions():
     assert upretty(fresnelc(x)) == 'C(x)'
     assert pretty(fresnels(x)) == 'S(x)'
     assert upretty(fresnels(x)) == 'S(x)'
-    assert pretty(Heaviside(x)) == 'Heaviside(x)'
-    assert upretty(Heaviside(x)) == 'θ(x)'
+    assert pretty(Heaviside(x)) == 'Heaviside(x, 1/2)'
+    assert upretty(Heaviside(x)) == 'θ(x, 1/2)'
path: sympy/printing/pretty/tests/test_pretty.py
text: @user1:
I wonder if the extra arg should be omitted in the printing for the default value of `1/2`.

@author:
That would be a very good idea, but I do not know where to do it. Would I need to write a `def _print_Function` in the `class Heaviside`?

@user1:
It's the other way round. You add a `_print_Heaviside` method to the pretty printer:
https://github.com/sympy/sympy/blob/92e7cb0b27bb810d67013a2cdbbd1f8e7f7aba45/sympy/printing/pretty/pretty.py#L1574-L1576

@author:
OK, I did it ;)  Push comes as soon as I have run all checks I think I need to run locally.

@author:
I wondered whether we should also change the normal printer, but I think we should just leave it as it is, such that people using this will see immediately that `Heaviside` has changed. Is that OK?

@user1:
What do you mean by the normal printer? SymPy has so many printers...

@author:
The printer that activates for in-line documentation and is tested through `bin/doctest`. I don't actually know which one that is ;)

P.S. It is the same that does `print(Heaviside(x))`

Since all went through now, I suggest to take what we have now, proceed with DiracDelta, and later return to Heaviside if necessary.

diff_hunk: @@ -404,44 +404,46 @@ def _sage_(self):
 
 
 class Heaviside(Function):
-    r"""
-    Heaviside Piecewise function.
+    r"""Heaviside step function.
+
+    .. versionchanged:: 1.9
+        ``Heaviside(0)`` now returns 1/2; previously, the default was that
+        the value at 0 was undefined.
path: sympy/functions/special/delta_functions.py
text: @user1:
We don't have anything about this in the style guide but I wouldn't put the versionchanged at the top of the docstring. Probably it should go at the bottom.

@author:
Will be moved ;)