Piecewise functions in Mathematica

Adam Strzebonski

Wolfram Research, Inc.

Abstract

Mathematica 5.1 contains extensive support for piecewise functions throughout the system. An arbitrary piecewise function (with a finite number of pieces) can be represented using the new Piecewise construct. The function PiecewiseExpand allows to transform an arbitrary composition of piecewise functions into a single piecewise function. These tools, as well as Mathematica's condition solving capabilities, allow for efficient handling of piecewise functions in many contexts, including equation and inequality solving, optimization, computation of limits, derivatives, series, indefinite and definite integrals, solving of differential equations, and graphics.

Motivation

• piecewise in teaching

• piecewise in numerics (e.g. interpolation, quadrature, ode solvers, pde solvers, ...)

• piecewise in graphics (e.g. lines, polygons, splines, nurbs, ...)

• piecewise in electrical engineering (e.g. nonlinear circuits, step responses, ...)

• piecewise in control engineering (e.g. step responses, hybrid systems, fuzzy control, ... )

• piecewise in biomedical engineering (e.g. neuron firing, gene production, ...)

• piecewise in symbolic computation (e.g. equation and inequality solving, optimization, limits, derivatives, series, indefinite and definite integrals, differential equations, ...)

Piecewise Definition and Evaluation

This defines a piecewise function.

pw =  
      Piecewise[{{Sin[x]/x,  x < 0},  {1, x == 0}}, -x^2/100  + 
          1]

It can be evaluated as any mathematical function.

pw /. {{x -> -5}, {x -> 0}, {x -> 5}}

Plot[ pw,  {x,  -10, 10}]

Note that evaluating pw at zero does not produce error messages, because a value expression in a Piecewise function is evaluated only if the corresponding condition does not evaluate to False, and no previous condition evaluates to True.

Piecewise Expansion

The following property is fundamental for symbolic computation with piecewise functions and a key to deep support for piecewise functions in Mathematica.

An expression is in piecewise normal form if it is a single piecewise function, with conditions and values free of piecewise functions.

Property: An arbitrary expression containing piecewise functions can be expanded out into an equivalent expression in piecewise normal form.

Proof: Let E be an expression containing piecewise functions. Let C be the set of all conditions of all piecewise functions in E, with each condition involving piecewise functions replaced with a finite collection of non-piecewise conditions obtained by replacing the piecewise functions with all possible combinations of their value expressions. Now we can construct a single piecewise function F equivalent to E by taking its conditions to be the set of all possible combinations of truth values of elements of C. The corresponding value expressions can be obtained from E by noting that whenever all conditions in C have fixed truth values, each of the piecewise functions in E simplifies to one of its value expressions.

By piecewise expansion we will mean transforming an expression containing piecewise functions into an equivalent expression in piecewise normal form. The proof above gives an effective method of computing piecewise expansion. However, this method is rarely practical. Usually Mathematica computes piecewise expansion by a depth-first expansion of subexpressions based on the following three transformation rules, combined with simplification of intermediate results.

Distributive Law 1: Value Denesting.

Distributive Law 2: Condition Denesting.

Distributive Law 3: Function Composition.

PiecewiseExpand

This computes piecewise expansion of a nested piecewise function.

(Subscript[pw, 1]) = 
      Piecewise[{{Piecewise[{{1, x ≥ 0}}, 2], 
            Piecewise[{{x, x ≤ 1}}, x/2]^2 ≥ 1/2}}, 
        3]

PiecewiseExpand[(Subscript[pw, 1])]

Complexity of Piecewise Expansion

A piecewise normal form of F[(Subscript[pw, 1]), (Subscript[pw, 2]),  ... , (Subscript[pw, k])], where each (Subscript[pw, i]) has (Subscript[n, i]) piecewise cases will in the worst case have in (Subscript[n, 1]) (Subscript[n, 2]) ⋯ (Subscript[n, k]) piecewise cases.

Here is an example displaying the worst case complexity.

spw[k_] := 
      Piecewise[{{2^(k - 1),  (Subscript[c, k])}}]
PiecewiseExpand[spw[1] + spw[2] + spw[3]]

Table[
      Length[PiecewiseExpand[Sum[spw[i], {i, n}]][[1]]] + 1, {n, 
        10}]

Could we use univariate linear inequality conditions and still get an exponential number of cases?

f[i_Integer, n_Integer, x_] := 
      Piecewise[{{(2^(i - 1)), 
            Or @@ Table[
                j <= x < j + (2^(i - 1)), {j, (2^(i - 1)), (2^n), 
                  (2^i)}]}}]
Table[
      Length[PiecewiseExpand[Sum[f[i, n, x], {i, n}]][[1]]] + 1, {n, 
        10}]

Here are graphs of f[1, n, x], ..., f[n, n, x], for n=4. The graph of their sum is given in red.

Show[
      Append[Table[
          Plot[Evaluate@f[i, 4, x], {x, 0, 16}, PlotRange -> {0, 16}, 
            PlotStyle -> Hue[i/4 - 1/8], 
            DisplayFunction -> Identity], {i, 
            4}], 
Plot[
          Evaluate@Sum[f[i, 4, x], {i, 4}], {x, 0, 16}, 
          DisplayFunction -> Identity, PlotStyle -> Red]], 
      DisplayFunction -> $DisplayFunction]

The sum indeed has 16 different cases.

PiecewiseExpand[Sum[f[i, 4, x], {i, 4}]]

But the complexity here is exponential only if we measure it based solely on the number of cases. If we look at the number of inequalities used to specify the conditions in f[1, n, x], we see that the input to PiecewiseExpand already has an exponential size.

f[1, 4, x]

In fact, if the conditions in all piecewise functions involved are univariate polynomial equations and inequalities in the same variable, which is quite often the case in practice, we have a much more optimistic estimate of the complexity of piecewise expansion. Namely the worst case complexity of piecewise expansion is linear in the size of input measured by the number of equations and inequalities it contains.

Suppose the piecewise conditions in (Subscript[pw, 1]), (Subscript[pw, 2]),  ... , (Subscript[pw, k]) are univariate polynomial equations and inequalities in a real variable x. Let n be the number of real roots of all polynomials present in the conditions. Then any expression F[(Subscript[pw, 1]), (Subscript[pw, 2]),  ... , (Subscript[pw, k])] has a piecewise normal form with at most 2 n+1 cases.

Simplification of Piecewise Functions

Piecewise Removal

Piecewise[{{e1,  True},  {e2,  d2},  {e3,  d3}},  e]

Dead Case Removal (True)

Piecewise[{{e1,  d1},  {e2,  d2},  {e3,  True},  {e4,  
          d4},  {e5,  d5}}, e]

Dead Case Removal (False)

Piecewise[{{e1,  d1},  {e2,  d2},  {e3,  False},  {e4,  
          d4},  {e5,  d5}}, e]

Dead Case Removal (Implication)

Piecewise[ {{e1,  d1},  {e2,  d2},  {e3,  d2 && d3},  {e4,  
          d4}}, e]

Piecewise[ {{e1,  d1},  {e2,  d2 || d3},  {e3, d3},  {e4,  d4}}, 
      e]

Case Merging

Piecewise[ {{e1,  d1},  {e23,  d2},  {e23,  d3},  {e4,  d4}}, 
      e]

Reordering Transformation

Piecewise[{{e[l, a, r, g, e], d1}, {e2, d2}, {e[l, a, r, g, e], d3}}, 
        e] // Simplify

Orderless Transformation

This transformation is important in many applications. The resulting piecewise function has all cases disjoint and therefore the cases can be arbitrarily reordered.

PiecewiseExpand[ 
      Piecewise[ {{e1,  d1},  {e2,  d2},  {e3,  d3}}],  
      Method ->  {"OrderlessConditions"  ->  
            True}]

Mathematica Applications of Piecewise Functions

• special piecewise functions and simplification

• equations, inequalities and quantifiers

• optimization

• derivatives, limits and series

• indefinite and definite integration

• differential equations

• visualization

Finite Special Piecewise Functions

There are already many special piecewise functions in Mathematica.

They are all special cases of Piecewise, and can be transformed to Piecewise using PiecewiseExpand.

PiecewiseExpand[Max[x, Min[y, z]]]

In this case the expansion will only happen if you assume your variables to be real, since for complex arguments Abs is not a piecewise function.

PiecewiseExpand[Max[Abs[x],  Abs[y]],  Reals]

Complicated combinations of special piecewise functions may often have a simple piecewise normal form representation.

(Subscript[pw, 2]) = 
      Nest[Max[(# + 2  UnitStep[x] - 7),  #^2 + 1]/(#^2 - 
                UnitStep[x - 1] + 3) &, UnitStep[x], 3]

PiecewiseExpand[(Subscript[pw, 2])]

Infinite Special Piecewise Functions

Another important group of special piecewise functions consist of functions with potentially infinitely many piecewise cases. These can be transformed to Piecewise if the restrictions on their argument make the number of possible cases finite.

Here the restriction on x makes the number of possible cases finite.

PiecewiseExpand[Mod[x^2, Ceiling[x]], 1 < x < 3]

In this case Mathematica can compute finite bounds on x and hence peform the expansion.

PiecewiseExpand[0 ≤ x  -  Floor[x]^2 ≤ 3,  
      Reals]

Piecewise Equations and Inequalities

Mathematica can solve systems of equations and inequalities involving piecewise functions, and can solve quantifier elimination problems for such systems. This can be done provided the piecewise conditions are themselves equations and inequalities.

This solves two inequalities involving special piecewise functions.

Reduce[1  ≤ x  +  Ceiling[x + Floor[x]] < 10,  
      x]

This solves an equation involving a piecewise function.

Reduce[
      (x^2)  +  
          Piecewise[{{x, 0 < x < 1}, {(x^3), 
                1 ≤ x < 2}}]  ==  5, x]

This finds a particular solution to a system of inequalities and equations.

FindInstance[
      2 ≤ Max[(x^2), (y^2)] ≤ 3 && 
        y == Floor[x] + 1, {x, y}, Reals]

This finds the conditions on a such that resulting inequalities have a solution, i.e. performs quantfier elimination.

Resolve[(Subscript[∃ , x]) Piecewise[{{x + 7, (x^2) < a}}, (x^3) - 2 a]^2 + 
            x ≤ 2, Reals]

This tells for which a a square with side a fits in the unit circle.

\(Resolve[∀ \_{x, y}\((x\^2 + y\^2)\)\ Boole[
              Max[Abs[x], Abs[y]] < a\/2] < 1, Reals]\)

Piecewise Optimization

Mathematica solve optimization problems involving piecewise functions in the objective function as well as in constraints. This can be done provided the piecewise conditions are equations and inequalities.

This minimizes a piecewise function over an interval.

Minimize[{3 + (((x - 2))^UnitStep)[ 1 -  Abs[x - 3/2]] ,  
        0 ≤ x < 4},  {x}]

This maximizes a nested piecewise function over piecewise defined inequality constraints.

Maximize[{Exp[ Floor[x]^Sign[x]],  
        0 < x + Ceiling[x]^3 < 12},  {x}]

Piecewise Derivatives

This computes the derivative of a piecewise function.

D[Piecewise[{{Sin[x]/x,  x < 0},  {1, x == 0}}, -x^2/
            100  + 1], x]

For univariate piecewise functions over the reals Mathematica checks that the derivatives at case boundaries indeed exist.

D[Piecewise[{{Sin[x]/x,  x < 0},  {0, x == 0}}, -x^2/
            100  + 1], x]

In this case the piecewise function is expanded and the derivative for each segment is given.

D[Max[x,  (x^2)],  x]

Piecewise Limit

This computes the limit at z->I of the characteristic function of the unit circle. By default the limit is taken in the direction of the negative real axis.

Limit[Piecewise[{{1,  Abs[z] ≤ 1}}], 
      z -> I]

This computes the limit in the direction of the positive imaginary axis.

Limit[Piecewise[{{1,  Abs[z] ≤ 1}}], z -> I,  
      Direction -> I]

For problems involving parameters the answer may be piecewise.

Limit[Max[a x^2 + 1, x], x -> 2]

Piecewise Series

This function is analytic at x==1.

Series[
      Exp[Piecewise[ {{(E^x), x < 0}, {Sin[x], 
              x ≥ 0}} ]],  {x,  1,  2}]

In this case the series around infinity is well defined.

Series[Min[x, (x^2)],  {x,  Infinity,  3}]

For problems involving parameters the answer may be piecewise.

Series[(E^Max)[3 (x^2) + 1, (x^3) + 3 x],  {x,  a,  
        2}]

Piecewise Indefinite Integrals

Here is an integral of a piecewise function.

int = 
      Integrate[
        Piecewise[{{(E^x), (x^2) ≤ 1}}, Sin[x]],  
        x]

Integration constants are chosen to make the result continuous (if Mathematica can find the case boundaries).

Plot[int, {x, -2, 2}]

This function is not holomorphic, and hence it cannot be integrated over the complex numbers.

Integrate[Abs[2 - Abs[x]],  x]

However, over the reals it is piecewise linear and can be integrated.

Assuming[x ∈ Reals, 
      int2 = Integrate[Abs[2 - Abs[x]],  x]]

The integrand is continuous, hence the integral is smooth ((C^1)).

Plot[int2, {x, -4, 4}]

Piecewise Definite Integrals over Regions

When integrating over a region defined by a logical combination of inequalities it is often convenient to use the new piecewise function Boole, where Boole acts as a characteristic function or indicator function, being 1 when its argument is True and 0 otherwise. In this case it is important to understand that the complexity depends not only on the function to be integrated, but also on the complexity of the region over which you are integrating.

This computes the area of a unit disk.

Integrate[
      Boole[(x^2) + (y^2) < 
          1],  {x,  -∞,  ∞},  {y,  -∞,  ∞}]

NIntegrate[
      Boole[(x^2) + (y^2) < 
          1],  {x,  -∞,  ∞},  {y,  -∞,  ∞}]

This integrates a polynomial over a region defined as the intersection of a ball and a double cone.

Integrate[(2  (x^2) + 3  (y^2)) Boole[
            (x^2) + (y^2) + (z^2) ≤ 1 ∧ 
              3  (x^2) + 3  (y^2) ≤ 
                (z^2)],  {x, -∞, ∞} , {y, -∞ ,  ∞},  {z,  -∞,  ∞}] // Timing

NIntegrate[(2  (x^2) + 3  (y^2)) Boole[
          (x^2) + (y^2) + (z^2) ≤ 1 ∧ 
            3  (x^2) + 3  (y^2) ≤ 
              (z^2)],  {x, -∞, ∞} , {y, -∞ ,  ∞},  {z,  -∞,  ∞}]

When there are parameters in the region description Mathematica will attempt to compute all possible cases. In this case different parameter values can result in three fundamentally different regions.

Integrate[
      Boole[a (x^2) + (y^2) < 1],  {x,  0,  1},  {y,  0, 1}]

The approach above works in full generality for piecewise algebraic inequalities, however in some cases it will also work for transcendental inequalities.

Integrate[(E^x) Boole[
          x + 2 > (E^x)], {x, -∞, ∞}]

In this case there is an infinite number of intervals defined by the inequlities.

Integrate[
      Boole[1/2 < Sin[x] < 2/3]/
        (x^2), {x, -∞, ∞}]

General Piecewise Definite Integrals

For most purposes using the Boole construct is sufficient. However, Mathematica can work with arbitrary piecewise functions, as long as it can compute the integrals over the regions defined by the piecewise conditions.

Here is an integral of a function with 319 discontinuities in the integration interval.

Integrate[Piecewise[{{(x^2) , Sin[(x^3)] > 0}}], {x, 0, 10}]

Here is a plot of the integrand.

Plot[Piecewise[{{(x^2) , Sin[(x^3)] > 0}}], {x,  0,  10},  
      PlotPoints  ->  1000]

This integrates a nested piecewise function over a rectangle.

Integrate[
        Max[Abs[ (x^2) - y], Abs[x - (y^2)]] 
          Boole[Abs[x] + Abs[y] < 
              1], {x,  -∞,  ∞},  {y,  -∞,  ∞} ] // Timing

One can integrate many other types of special piecewise functions.

Integrate[
      Exp[Floor[(x^2) + y]] 
        Boole[Abs[x] + Abs[y] < 
            1], {x,  -∞,  ∞},  {y,  -∞,  ∞} ]

In this case we have nested piecewise functions both in the integrand and in the region description.

Integrate[
      FractionalPart[(x^2) + 2  Floor[x]] 
        Boole[1 < x + Floor[x] < 
            5], {x,  -∞,  ∞}]

Here is an example involving six different piecewise functions.

Integrate[
      Sin[(x^2) UnitStep[x] Log[Ceiling[Abs[x - 2  y]]]] 
        Boole[Abs[x] + Abs[y] < 
            1], {x,  -∞,  ∞},  {y,  -∞,  ∞} ]

Piecewise Differential Equations

Mathematica can solve piecewise defined differential equations. In the current implementation the piecewise conditions should depend only on the independent variable.

This solves a differential equation that is piecewise defined.

DSolve[ y'[t] == 
          Piecewise[{{y[t], t ≤ 2}}, -y[t]/2]  &&  
        y[0] == 1,  y,  t]

In this case the equation is formulated using special piecewise functions.

DSolve[y'[t]  +  Max[t, (t^2)] y[t]  ==  
          Boole[t ≥ 0]  &&  y[0] == 1,  y[t],  
      t ]

Example: Piecewise Polynomial Interpolation

Implementation

\(\(Options[
          PiecewiseInterpolatingPolynomial] = {InterpolationOrder -> 
            3};\)
\)
PiecewiseInterpolatingPolynomial[data_List, x_, opts___?OptionQ] := 
      Module[{d, pdata}, 
          d = InterpolationOrder /. {opts} /. 
              Options[PiecewiseInterpolatingPolynomial]; 
pdata = 
            Partition[data, d + 1, 
              d, {1, 1}, {}]; 
Piecewise[
            Transpose[{Map[InterpolatingPolynomial[#, x] &, pdata], 
                Map[#[[1, 1]] ≤ x < #[[-1, 1]] &, 
                  pdata]}]]] /; ArrayDepth[data] == 2

Examples

Runge's example of failure of interpolating polynomial to converge when degrees increase.

RungeInterpolate[ n_,  x_ ]  :=  
      InterpolatingPolynomial[ 
        Table[ {x,  (1/(1 + 25  (x^2)))},  {x,  -1. ,  1. ,  
            2. /(n - 1)}],  x ]

Table[ 
        Plot[ Evaluate @ {(1/(1 + 25 (x^2))),  
              RungeInterpolate[n,  x]},  {x,  -1,  1},  
          PlotStyle ->  {{Red},  {Blue}},  
          PlotRange  -> {-1, 1}],  {n,  5,  20,  
          5}];

RungeInterpolate[10, x]

Here is Runge's example, using piecewise polynomial interpolation instead. As can be seen we get very fast convergence across the whole interval.

RungePiecewiseInterpolate[ n_,  x_ ]  :=  
      PiecewiseInterpolatingPolynomial[ 
        Table[ {x,  (1/(1 + 25  (x^2)))},  {x,  -1. ,  1. ,  
            2. /(n - 1)}],  x ]

Table[ 
        Plot[ Evaluate @ {(1/(1 + 25 (x^2))),  
              RungePiecewiseInterpolate[n,  x]},  {x,  -1,  1},  
          PlotStyle ->  {{Red},  {Blue}},  
          PlotRange  ->  {{-1, 1}, {-0.2, 1.2}}],  {n,  5,  
          20,  5}];

RungePiecewiseInterpolate[10, x]

Example: Lebesgue integration

This is an illustration of step function approximations used to define the Lebesgue intergal,
and Lebesgue's convergence theorems. We define two sequences of step functions pointwise convergent to the given function f(x) on the interval [a, b]. We assume that 0≤f(x)≤M.

\(LowerApprox[f_, x_, n_, {a_, b_}, M_] := 
      Piecewise[
        Table[{k\ M/2^n, 
            Reduce[k\ M/2^n ≤ f < \((k + 1)\) M/2^n && 
                a ≤ x ≤ b, x]}, {k, 
            2^n}]]
\)
UpperApprox[f_, x_, n_, {a_, b_}, M_] := 
      Piecewise[
        Table[{k M/2^n, 
            Reduce[(k - 1) M/2^n < f ≤ k M/2^n && 
                a ≤ x ≤ b, x]}, {k, 2^n}]]

The functions of the first sequence are increasing and <= f, so their integrals tend to the integral of f from below by Lebesgue's monotonic convergence theorem. The functions of the second sequence are decreasing , <= M, and > f, so their integrals tend to the integral of f from above by Lebesgue's bounded convergence theorem.

f = 2 (Sqrt[1 - (x^2)]);

Table[
        Plot[Evaluate@{LowerApprox[f, x, i, {-1, 1}, 2], 
              UpperApprox[f, x, i, {-1, 1}, 2], f}, {x, -1, 1}, 
          PlotRange -> {{-1, 1}, {0, 2}}, 
          PlotStyle -> {{Blue}, {Red}, {Green}}], {i, 3, 
          6}];

Let's see how close are the integrals of the step functions to the value of the integral of f.

measure[s_Or] := Plus @@ (measure /@ List @@ s)
measure[s_Inequality] := s[[5]] - s[[1]]
measure[s_Equal] := 0
StepIntegrate[f_Piecewise] := 
      Plus @@ (#[[1]]*measure[#[[2]]] & /@ 
            f[[1]])

(Subscript[l, 6]) = 
      StepIntegrate[LowerApprox[f, x, 6, {-1, 1}, 2]] // 
        Expand

(Subscript[u, 6]) = 
      StepIntegrate[UpperApprox[f, x, 6, {-1, 1}, 2]] // 
        Expand

N[{(Subscript[l, 6]), (Subscript[u, 6])}]

Summary

• piecewise used everywhere

• piecewise normal form

• piecewise symbolic,numeric and graphic solvers

• easy to do prototyping,research and teaching

• many more possible applications