3
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ Why you don't switch in rational numbers? $\endgroup$
    – Dimitris
    May 28 '15 at 15:34
  • $\begingroup$ I can do that, but eventually I'll need to make it work with real numbers. $\endgroup$
    – psachdeva
    May 28 '15 at 15:39
  • $\begingroup$ 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}]; $\endgroup$
    – Histograms
    May 28 '15 at 15:42
1
$\begingroup$

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]

Mathematica graphics

$\endgroup$
1
  • $\begingroup$ Dr.belisarius (@ does not work.): This is not a workaround since the continuity of the derivative is used in deriv[fun_, pp_] := Limit[D[fun[p], p], p -> pp], but this continuity is not established. $\endgroup$
    – user64494
    Mar 11 at 17:38
2
$\begingroup$

You might also find the NumericalCalculuspackage helpful. In particular the function ND

Needs["NumericalCalculus`"]

ND[edens[p], p, p0]
(* 155.59 *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.