# Derivatives of piecewise functions of functions

### Background

Derivatives of piecewise functions in Mathematica are computed according to special rules. According the Piecewise documentation (see Possible Issues),

Derivatives are computed piece-by-piece, unless the function is univariate in a real variable.

This distinction is important for piecewise functions that have a pointwise definition (such as a discontinuity). In this case, if Mathematica can determine that the variable is real, it finds the derivative at the discontinuity by looking at the derivatives on each side, rather that at the function definition at the discontinuity. This is somewhat discussed here.

In general, this behavior is intelligent and useful. Consider the function

$$f(b, p) = \sqrt{\frac{b - 1}{e^{(b - 1)p} - 1}}$$

which is, in fact, continuous, but evaluates to Indeterminate at $b = 1$. Its derivatives with respect to the first argument will have the same problem. Now, If we define $f$ via Piecewise, we can put the appropriate limit in by hand.

f[b_, p_] := Piecewise[
{{Sqrt[(b - 1)/(Exp[(b - 1) p] - 1)], b < 1 || b > 1},
{Sqrt[1/p], b == 1}}
]


To get derivatives, we can differentiate as normal:

D[f[b, p], b] // Simplify


$$\begin{cases} -\frac{1}{4 \sqrt{\frac{1}{p}}} & b = 1 \\[6pt] \frac{e^{(b-1) p} (1 + p - b p) - 1}{2 \sqrt{\frac{b-1}{e^{(b-1) p} - 1}} \left(e^{(b-1) p}- 1 \right)^2} & \text{True} \end{cases}$$

Even though the function value at $b = 1$ does not depend on $b$, the derivative is correctly found because Mathematica looks at the surrounding region. Multiple derivatives are fine, too.

### The Problem

However, if the variable defining the piecewise regions is a function of the variable of differentiation, the differentiation always proceeds piece-by-piece.

D[f[g[b], p], b] // Simplify


$$\begin{cases} 0 & g(b) = 1 \\[6pt] \frac{e^{(g(b)-1) p} (1 + p - g(b) p) - 1}{2 \sqrt{\frac{g(b)-1}{e^{(g(b)-1) p} - 1}} \left(e^{(g(b)-1) p}- 1 \right)^2} g^\prime (b) & \text{True} \end{cases}$$

The derivative at the discontinuity vanishes. Note that D[f[g, p], b, NonConstants -> {g}] and Dt[f[g, p], b, Constants -> {b}] also exhibit this behavior. Neither does it matter that g is unspecified (g = Sin will do the same thing).

It appears that this occurs because with a compound variable Piecewise no longer knows that the variable is real (I'm not sure why), so the derivative is computed piece-by-piece, and the derivative of the constant (in $b$) piece at $b = 1$ goes to zero. One can check this by adding some $b$-dependence and verifying that it behaves as expected for piece-by-piece differentiation:

h[b_, p_] := Piecewise[
{{Sqrt[(b - 1)/(Exp[(b - 1) p] - 1)], b < 1 || b > 1},
{b Sqrt[1/p], b == 1}}
]

D[h[g[b], p], b] // Simplify


$$\begin{cases} \sqrt{\frac{1}{p}} \, g^\prime (b) & g(b) = 1 \\[6pt] \frac{e^{(g(b)-1) p} (1 + p - g(b) p) - 1}{2 \sqrt{\frac{g(b)-1}{e^{(g(b)-1) p} - 1}} \left(e^{(g(b)-1) p}- 1 \right)^2} g^\prime (b) & \text{True} \end{cases}$$

Of course, this is only an issue for piecewise functions that have a pointwise definition that doesn't carry the full variable dependence necessary to compute derivatives.

Is this behavior intended? How do I work around it?

Ideally, one should be able to take as many derivatives of $f(g(b), p)$ with respect to $b$ (providing $g$ is well-behaved) as are possible of $f(b, p)$. For the purposes of this question, treat $g$ as an unspecified, but well-behaved function -- our lack of knowledge about its behavior should be parametrized by the appearance of $g^\prime(b)$ and higher derivatives in the output. Essentially, Mathematica is unable to work out the chain rule for piecewise functions with a pointwise definition, since it differentiates them piece-by-piece, and I want to know the best way to deal with this.

I ran into this problem because I need to be able to take derivatives of a piecewise function of a dependent variable with respect to an independent variable so that I can use equations containing that function in NDSolve. Practically, it is easy enough to get one or two correct derivatives by defining the function value at the tricky point as some number of terms in its series expansion to give it explicit variable-dependence, but I am interested in a more general solution.

Any suggestions are welcome.

• How do you know if a general function of g[b] is differentiable at all? I mean, when we do not know the definition we can't guarantee the derivative exists at that point for the point-wise function. If it exists, and is equal from left and right, i guess we can define a function to use the definition of a derivative using limits, or, do a numerical differentiation? – MathX Apr 2 '16 at 3:27
• @MathX While you wrote this, I answered with essentially the same general thought... – Jens Apr 2 '16 at 3:28
• @Jens yeah I was editing when I saw an answer was posted haha you beat me to it. This question actually made me curious to see if Mathematica has a built-in Frechet derivative function which apparently it doesn't. – MathX Apr 2 '16 at 3:33
• @MathX In this case, I know that the function g[b] is differentiable. – Virgil Apr 2 '16 at 4:01
• @MathX also, take a look at DifferenceQuotient. – Virgil Apr 2 '16 at 4:03

Obviously you can get derivatives only if the function actually has a derivative, and this is true for the example in the question. Your function is even continuous at the point in question. This means singularity is removable, and therefore we should not have to use Piecewise for that single point at all. Here is how you can avoid its use when you suspect that the derivative exists everywhere:

f[x_, p_] :=
Module[{b}, Limit[Sqrt[(b - 1)/(Exp[(b - 1) p] - 1)], b -> x]]

D[f[b, p], b] // Simplify

(*
==> (-1 +
E^((-1 + b) p) (1 + p - b p))/(2 Sqrt[(-1 + b)/(-1 +
E^((-1 + b) p))] (-1 + E^((-1 + b) p))^2)
*)

D[f[g[b], p], b] // Simplify

(*
==> -(((1 - E^(p (-1 + g[b])) (1 + p) +
E^(p (-1 + g[b])) p g[b]) Derivative[g][b])/(
2 (-1 + E^(p (-1 + g[b])))^2 Sqrt[(-1 + g[b])/(-1 + E^(
p (-1 + g[b])))]))
*)


The trick is simply to evaluate the function as a Limit. With that, the result contains no case distinctions.

Edit

Having the derivatives as above, you may then of course encounter the same removable singularity and therefore would need to evaluate these results using Limit as well. I.e., to do the substitution mentioned in the comment, you would have to do

Limit[D[f[g[b], p], b], g[b] -> 1]

(* ==> -(Derivative[g][b]/(4 Sqrt[1/p])) *)


Edit 2

The comments made it clear that the a formulation in terms of Piecewise is indeed preferred. So here is an alternative approach that relies on the idea (mentioned in the comment) that f itself may only be needed for numerical arguments:

ClearAll[f];
fOrig[x_, p_] :=
D[Piecewise[{{Sqrt[(x - 1)/(Exp[(x - 1) p] - 1)],
x < 1 || x > 1}, {Sqrt[1/p], x == 1}}], x]

Derivative[1, 0][f][b_, p_] :=
Module[{x}, Simplify[fOrig[x, p]] /. x -> b]

f[b_?NumericQ, p_] := fOrig[b, p]

D[f[b, p], b]


$$\begin{cases} -\frac{1}{4 \sqrt{\frac{1}{p}}} & b=1 \\ \frac{e^{(b-1) p} (-b p+p+1)-1}{2 \sqrt{\frac{b-1}{e^{(b-1) p}-1}} \left(e^{(b-1) p}-1\right)^2} & \text{True} \end{cases}$$

D[f[g[b], p], b]


$$g'(b) \begin{cases} -\frac{1}{4 \sqrt{\frac{1}{p}}} & g(b)=1 \\ \frac{e^{p (g(b)-1)} (-g(b) p+p+1)-1}{2 \left(-1+e^{p (g(b)-1)}\right)^2 \sqrt{\frac{g(b)-1}{-1+e^{p (g(b)-1)}}}} & \text{True} \\ \end{cases}$$

This gives the correct results because I defined the derivative directly to use the first differentiation that worked properly (without inner function), and then manually added the chain rule.

• The difficulty is that I want the derivatives, which are also continuous, to give the appropriate limit at the point in question. However with this Limit-based definition, (D[f[g[b], p], b]) /. g[b] -> 1 returns Indeterminate. – Virgil Apr 2 '16 at 3:59
• @Virgil The substitution g[b] -> 1 isn't going to work anyway because it leaves g' untouched. But I think I understand what you mean. You want to keep g unspecified until the end, not give it a definition at the beginning, right? Then I guess the limit approach wouldn't help. I think it only works if the inner function is already defined so that the singular point in the outer function can actually be identified. – Jens Apr 2 '16 at 4:45
• @Virgil I added the "workaround" for the case you mentioned, showing the extent to which Limit can do what you ask. However, it's probably not as "automated" as you want. However, not knowing anything about g, I don't see a better alternative right now. – Jens Apr 2 '16 at 4:57
• Correct. In the end, g is going to be a dependent variable in NDSolve, so it needs to be unspecified. I also need the equations it occurs in (which include the piecewise function and derivatives thereof) to be well-defined (as in not Indeterminate) everywhere. I played around with Limit myself for a while, but abandoned it for Piecewise, since it cannot pass along the limit to derivatives. – Virgil Apr 2 '16 at 5:06
• I would like to automate it. Perhaps a new myLimit-type function that has HoldFirst, passes D inside, and evaluates only for numerical arguments would work – Virgil Apr 2 '16 at 5:10

I have accepted Jens' method, since it works and with a minimum of inconvenience. However, it is more useful when generalized to multiple derivatives, so I will show what I did to accomplish that below. I also want to present two other work-arounds that I came up with. I do not like either as well as Jens', but they may be useful in cases where one does not want to use a numeric function.

## Jens' method for general derivatives

To correctly obtain more complicated derivatives, we simply need to make a general Derivative definition for $f$ that is numeric in the appropriate argument:

fOrig[b_, p_] := Piecewise[
{{Sqrt[(b - 1)/(Exp[(b - 1) p] - 1)], b < 1 || b > 1},
{Sqrt[1/p], b == 1}}
]

f[b_?NumericQ, p_] := fOrig[b, p]

Derivative[n__][f][b_?NumericQ, p_] := Derivative[n][fOrig][b, p]


This is the approach I have decided to take in my own work.

## A special derivative

It one want's to avoid numeric definitions as in the above, one option is to define a special derivative that uses Inactivate and Activate to hold particular functions inactive during the course of derivation, restoring them after any chain and/or product rules have been worked out. The effect is the same as that achieved by making the argument numeric.

SetAttributes[myD, HoldFirst];
Options@myD = Options@D;
myD[expr_, patt_, x__, opts : OptionsPattern[]] := Activate[
D[Inactivate[expr, patt], x, opts],
patt]


Here, the argument patt takes a pattern, and any parts of expr matching that pattern are inactivated before D is applied and reactivated after. It can be applied as so:

myD[fOrig[g[b], p], fOrig, b]


$$g^\prime (b) \, \begin{cases} - \frac{1}{4\sqrt{\frac{1}{p}}} & g(b) = 1 \\[6pt] \frac{e^{(g(b)-1) p} (1 + p - g(b) p) - 1}{2 \sqrt{\frac{g(b)-1}{e^{(g(b)-1) p} - 1}} \left(e^{(g(b)-1) p}- 1 \right)^2} & \text{True} \end{cases}$$

Higher-order derivatives work as well.

I am not as fond of this method as it requires the use of a special new derivative, which I find cumbersome.

## A special limit

One other thing one can do is to work with the piece-by-piece differentiation of Piecewise and modify the limiting value itself, so that when differentiated it gives the proper result. I came up with this:

SeriesLimit[f_, {x_, x0_?NumericQ}] := SeriesCoefficient[f, {x, x0, 0}]
Derivative[0, {0, n_}][seriesLimit][f_, {x_, x0_}] := seriesLimit[D[f, {x, n}], {x, x0}]
Derivative[n_, {0, 0}][seriesLimit][f_, {x_, x0_}] := seriesLimit[1, {x, x0}]
Times[a___, seriesLimit[f_, {x_, x0_}], b___] ^:= seriesLimit[a f b, {x_, x0_}]


seriesLimit is a wrapper that holds its arguments until the argument x0 is numeric, at which point it takes the first term in the series expansion of $f$ around $x_0$, that is, $f(x_0)$. Under differentiation, it passes derivatives with respect to x0 inside to f as derivatives with respect to x, and collects terms generated by derivatives with respect to other arguments back inside. This behavior is probably not robust, but it is sufficient here.

D[seriesLimit[g[x, p], {x, b}], b]

seriesLimit[Derivative[1, 0][g][x, p], {x, y}]

D[seriesLimit[g[x, p], {x, h[b]}], b, b]

seriesLimit[h''[b] Derivative[1, 0][g][x, p], {x, h[b]}] + seriesLimit[h'[b]^2 Derivative[2, 0][g][x, p], {x, h[b]}]

seriesLimit[Sqrt[(x - 1)/(Exp[(x - 1) p] - 1)], {x, 1}]

Sqrt[1/p]


It can be used in Piecewise to alleviate the difficulty at hand:

fNew[b_, p_] := Piecewise[
{{Sqrt[(b - 1)/(Exp[(b - 1) p] - 1)], b < 1 || b > 1},
{seriesLimit[Sqrt[(x - 1)/(Exp[(x - 1) p] - 1)], {x, b}], b == 1}}
]

D[fNew[g[b], p], b]


$$\begin{cases} g^\prime (b) \, \frac{e^{(g(b)-1) p} (1 + p - g(b) p) - 1}{2 \sqrt{\frac{g(b)-1}{e^{(g(b)-1) p} - 1}} \left(e^{(g(b)-1) p}- 1 \right)^2} & g(b) > 1 \| \, g(b) < 1 \\[6pt] \text{seriesLimit} \Big[ g^\prime (b) \, \frac{e^{(x-1) p} (1 + p - x p) - 1}{2 \sqrt{\frac{x-1}{e^{(x-1) p} - 1}} \left(e^{(x-1) p}- 1 \right)^2}, \{x, g(b)\} \Big] & \text{True} \end{cases}$$

In the case that the piecewise variable is not a function (i.e. $g(b) \rightarrow b$), under differentiation the derivative at the discontinuity is calculated by the limiting procedure outlined in the question, and the seriesLimit drops out.