# Derivative of Continuous and Differentiable Piecewise function is indeterminate

I'm attempting to analytically stitch together two linear functions with two quadratic functions. The function I'm considering, then, is a piecewise function consisting of 4 parts. Below are the values I'm considering for the scenario as well as the function itself.

k1 = 1.41; k2 = 1.45; eps1 = 2*^8; eps2 = 3.5*^8; p0 = 1.04*^8; delP = 0.01*^8;
kmax = (1/2)(k1 + k2) + (k2 p0 + eps2 - k1 p0 - eps1)/delP;
edens[p_] := Piecewise[{{k1 p + eps1, p <= p0 - delP},
{(1/2)(kmax - k1)/delP * (p-p0)^2 + kmax(p-p0) + (k1 p0 + eps1) + (1/2)delP * (kmax - k1), p0-delP <= p <= p0},
{(1/2)(k2 - kmax)/delP * (p-p0)^2 + kmax(p-p0) + (k2 p0 + eps2) + (1/2)delP * (k2 - kmax), p0 <= p <= p0+delP},
{k2 p + eps2, p >= p0 + delP}}];


This function is continuous and differentiable at the boundary points, which you can see by plotting it:

Plot[edens[p], {p, p0 - 2*delP, p0 + 2*delP}, Exclusions -> None] However, if I try and evaluate the derivative of this function at the boundary points (p0, p0-delP, and p0+delP), I get Indeterminate even though the value is the same on both sides.

I have a feeling this is because I'm using approximate real numbers. Is there anyway to force Mathematica to say that the derivative at a boundary point is a specific value?

• Why you don't switch in rational numbers? – Dimitris May 28 '15 at 15:34
• I can do that, but eventually I'll need to make it work with real numbers. – psachdeva May 28 '15 at 15:39
• if you do this before you define the function it should work. {k1,k2,eps1,eps2,p0,delP,kmax}=Rationalize[{k1,k2,eps1,eps2,p0,delP, kmax}]; – Histograms May 28 '15 at 15:42

You may define:

deriv[fun_, pp_] := Limit[D[fun[p], p], p -> pp]
deriv[edens, p0]
(* 155.59 *)


In fact, this way you can also obtain a symbolic form for the derivative

Clear[k1, k2, eps1, eps2, p0, delP, kmax];
deriv[edens, p] // InputForm
(* Piecewise[{{k1, delP + p - p0 <= 0},
{kmax + ((-k1 + kmax)*(p - p0))/delP, delP + p - p0 >= 0 && p - p0 <= 0},
{kmax + ((k2 - kmax)*(p - p0))/delP, p - p0 >= 0 && delP - p + p0 >= 0},
{k2, delP - p + p0 <= 0}},
0] *)


To get a plot you could use:

Plot[deriv[edens, p] /. p -> pp, {pp, p0 - 2*delP, p0 + 2*delP},  Exclusions -> None] You might also find the NumericalCalculuspackage helpful. In particular the function ND

Needs["NumericalCalculus"]

ND[edens[p], p, p0]
(* 155.59 *)
`