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$
  • $\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$

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.