34
$\begingroup$

Update: (1) By V11, not sure of the exact version, the derivative IntegerPart' has been given a symbolic definition. (2) The numeric derivative computed has changed from an order 8 approximation (see @acl's answer) to an order 16 one. I'm not sure how best to update this Q&A. -- @MichaelE2

Just found the following while debugging a problem.

Mathematica is calculating the derivative of IntegerPart[x] in some odd way:

Plot[{IntegerPart@u, D[IntegerPart[x], x] /.x -> u}, {u, 0, 3},PlotRange -> Full]

enter image description here

Do you know what is the violet curve in the plot above, and why Mathematica thinks it is the derivative of IntegerPart[x]?

Edit

Or consider for example the following: 

f[x_] := Sin[x] + IntegerPart[x];
Plot[{f[x], Derivative[1][f][x]}, {x, 0, Pi/2}]

enter image description here

Or similarly 

Plot[{FractionalPart[x], Derivative[1][FractionalPart][x]}, {x, 0, 3}]

enter image description here

What is happening here?

Improve this question
$\endgroup$
11
  • 1
    $\begingroup$ Strictly speaking IntegerPart[x] is not a differentiable function. $\endgroup$
    – nilo de roock
    Commented Jun 13, 2011 at 8:19
  • 4
    $\begingroup$ @ndroock1 It is differentiable in the space of distributions (functions like Dirac delta-function). The result roughly should be proportional to the Dirac grid - a discrete sum of Dirac delta-functions Sum[DiracDelta[x-n],{n,-Infinity,Infinity}] - DiracDelta[x] (no delta-function at zero). What is being seen on the plot looks like a discrete analog of this, along the lines Sin[k x]/(Pi x), when k->Infinity - this is one well-known limiting representation of a delta-function. $\endgroup$
    – Leonid Shifrin
    Commented Jun 13, 2011 at 9:08
  • $\begingroup$ @Leonid Indeed. Just try Integrate[Derivative[1][IntegerPart][x], {x, .9, 1.1}]. But that doesn't seem to explain the plot $\endgroup$
    – Dr. belisarius
    Commented Jun 13, 2011 at 11:52
  • $\begingroup$ @belisarius True. This seems to reflect an imperfection on the side of symbolic treatment of things related to differentiation of piecewise functions and the like. I am pretty sure that work is underway in WRI to address these issues. $\endgroup$
    – Leonid Shifrin
    Commented Jun 13, 2011 at 13:12
  • $\begingroup$ @belisarius, I was not aware that Derivative is able to numerically approximate the derivatives of completely arbitrary user-defined functions, but now I see that this is mentioned in the docs. The method it uses is not mentioned. What I don't understand though is why this gives a different result than yours: Plot[D[Quotient[x, 1], x] // Evaluate, {x, 0, 3}, PlotRange -> Full] $\endgroup$
    – Szabolcs
    Commented Jun 13, 2011 at 14:23

2 Answers 2

Reset to default
40
$\begingroup$

I've completely overhauled my answer. I believe this now answers the questions posed (why mma thinks the violet line is the derivative of IntegerPart'[x]).

Let's first look at ND, simply because its internals are easier to access and we may obtain some insight. Try:

Needs["NumericalCalculus`"]
nd[x_, opts___] := ND[IntegerPart[u], u, x, opts]
Manipulate[
 Plot[nd[x, Scale \[Rule] s, Terms \[Rule] n], {x, 1.5, 2.5}, 
  PlotRange -> Full],
  {{s, 1}, .01, 4},
  {{n, 5}, 1, 20, 1}]

which allows to vary the two parameters of ND, namely the Scale and number of Terms in the method it uses. Typical results look like this:

enter image description here

Now, the results of ND are different from those of N[D[...]], but allow us to guess what is happening.

In particular, ND is easy enough to "reverse-engineer"; all we have to do is ask nicely: 

ND[f[u], u, x, Terms -> 2, Scale -> 1] 
(*
 ->  4. (- f[0. + x] + f[0.5 + x]) - 1. (- f[0. + x] + f[1. + x])
*)

and playing with the Terms and Scale gives the game up: it is a finite difference method and the points at which the function to be differentiated is evaluated are always the same (for the same Terms and Scale). Notice how the derivative being calculated is always the derivative from the right.

Now, let us check (numerically) which points are evaluated when IntegerPart is differentiated a la belisarius. Define 

ClearAll[mip];
mip[x_?NumericQ] := (Sow[x]; IntegerPart[x])
ClearAll[getoffsets]
getoffsets[x_] := Reap[mip'[x]][[2, 1]] - x

then getoffsets[x] returns a list containing the offsets of the points used for evaluating the derivative from x. Thus, 

getoffsets[.9]
(*
->{0.,0.0526316, 0.105263, 0.157895, 0.210526, 0.263158, 0.315789, 0.368421, 0.421053,
     -0.0526316,-0.105263,-0.157895,-0.210526,-0.263158,-0.315789,-0.368421,-0.421053}
*)

In fact, these offsets are independent of the point where the derivative is evaluated: 

getoffsets[#] & /@ {.1, .5, .9, 1.1, 1.5, 143.} // Differences // Chop
(*
-> {{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},  
    {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
    {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}}
*)

(ie, the offsets from x of the points used for calculating the derivative at x are the same for all the x given in the list above).

Armed with this information, we can now plot the function IntegerPart'[x] together with vertical lines located at the points that differ from 1 (the singular point of IntegerPart'[x]) by any of the offsets. That is, at the vertical lines, one of the points used to estimate the derivate by mma crosses 1, where the value of IntegerPart[x] jumps:

 offsets = getoffsets[1.1];
 Plot[IntegerPart'[x],{x, .5, 1.5},
  Epilog -> (Line[{{#, 0}, {#, 10}}] & /@(1 - offsets)),
  PlotRange -> Full]

this gives

out

Evidently, the jumps in the plot occur whenever one of the points used to evaluate the derivative crosses 1, where IntegerPart[x] changes from 0 to 1. So it is completely natural that the plot has discontinuities there, given the method used to obtain the derivative.

To clarify: The numerical derivative is, effectively a sum of the form $\sum_i a_i \, f(x_i)$, with $f$ the function of which we are evaluating the derivative and $x_i$ the points given by x plus the $i$th offset (see list above).
Whenever x lies on one of the vertical lines, $f(x_i)$ for one of the possible $i$ jumps from 0 to 1, hence the value of the sum jumps by $a_i$. This is indeed what we are seeing.

With a bit more work one would probably be able to work out the coefficients being used in this finite difference scheme, but I think this already answers the question. If not, let me know what is missing.

Finally, consider what happens with a step function if we allow/do not allow the plotter to look at the symbolic form of the function:

ClearAll[mop];
mop[x_?NumericQ] := (Sow[x]; HeavisideTheta[x])
Plot[mop'[x], {x, -.5, .5}, PlotRange -> Full]
Plot[HeavisideTheta'[x], {x, -.5, .5}, PlotRange -> Full]

In the first plot, we prevent the plotter from seeing the symbolic expression; in the second, we do not. Here is what is produced:

enter image description here

So mma recognizes what is going on in the second case and plots the right thing. In the first it can't, and we get a behaviour similar to what belisarius saw.

It appears that compensating for such discontinuous functions is done "per case", is not done by looking at the numerical behaviour but by looking at the symbolic form, and IntegerPart isn't one of the covered cases.

Edit (by belisarius)

Elaborating a little over this great answer, one can find the finite differences coefficients that Mma uses for calculating the derivatives.

Assuming a centered differences method: 

ClearAll[mip];
mip[x_?NumericQ] := (Sow[x]; IntegerPart@x)
ClearAll[getoffsets]
getoffsets[x_] := Reap[mip'[x]][[2, 1]] - x;
k = getoffsets[1.];
Rationalize[
 Table[c[r], {r, Length@k}] k[[2]] /. ToRules@
   Chop@
    Reduce[
     And @@
      Join[
       Table[N[mip'[i]] == Sum[c[r] mip[i + k[[r]]], {r, Length@k}], 
                                                   {i, 1, 2, 1/Length@k}],
       Table[c[p] == -c[p + (Length@k - 1)/2], {p, 2, (Length@k + 1)/2}]],
     Table[c[r], {r, Length@k}]], 10^-11]
(* ->
{0, 8/9, -(14/45), 56/495, -(7/198), 56/6435, -(2/1287), 8/45045, -(1/ 102960),
  -(8/9), 14/45, -(56/495), 7/198, -(56/6435), 2/1287, -(8/45045),  1/102960}
*)

I am too lazy to compute the actual coefficients for the eighth order difference to check if Mma is using a standard algorithm, so I used Google to search for them

enter image description here So, everything is in place: @acl's answer is right, and Mathematica is using the standard eighth order centered approximation for calculating derivatives :)

Improve this answer
$\endgroup$
7
  • $\begingroup$ @acl Very nice answer. Hope you don't mind I added a new chapter (and please check if you agree on my update!) $\endgroup$
    – Dr. belisarius
    Commented Jun 14, 2011 at 17:24
  • $\begingroup$ @belisarius it should be... I did a similar calculation yesterday but didn't have time to either clean it up or make sure it's OK. so go ahead! $\endgroup$
    – acl
    Commented Jun 14, 2011 at 17:56
  • $\begingroup$ @acl Yep. It's right. See update $\endgroup$
    – Dr. belisarius
    Commented Jun 14, 2011 at 17:59
  • 3
    $\begingroup$ @belisarius actually it feels almost exactly like doing physics but without the paper- and grant-writing, political manoeuvring etc, of course... hmmm... $\endgroup$
    – acl
    Commented Jun 14, 2011 at 19:01
  • 1
    $\begingroup$ For reference, the method shown in this paper can be used to calculate the eighth order centered derivative approximation coefficients: With[{m = 1, s = 8, n = 16}, CoefficientList[Normal[Series[x^s Log[x]^m, {x, 1, n}]], x]] $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Feb 26, 2016 at 3:25
7
$\begingroup$

The violet curve is IntegerPart'[x], of course. Here's the behavior, outside of Plot:

In[37]:= f = Function[x, Evaluate[D[IntegerPart[x], x]]];

In[42]:= f /@ Range[0.5, 1.5, 0.025]

Out[42]= {0., 0., 0., 0., -0.000184538, -0.000184538, 0.00318987, \
0.00318987, -0.0263362, -0.0263362, 0.13901, 0.13901, -0.532708, \
-0.532708, 1.61679, 1.61679, -4.29432, -4.29432, 12.5946, 12.5946, \
12.5946, 12.5946, 12.5946, -4.29432, -4.29432, 1.61679, 1.61679, \
-0.532708, -0.532708, 0.13901, 0.13901, -0.0263362, -0.0263362, \
0.00318987, 0.00318987, -0.000184538, -0.000184538, -2.75067*10^-15, \
-2.75067*10^-15, -2.75067*10^-15, -2.75067*10^-15}

In[43]:= ListLinePlot[%]

enter image description here

Beyond that, I'm not sure how the derivative is being computed. Some sort of finite difference, perhaps.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ so it looks like Mathematica is able to approximate the numerical derivative of any function. Do you know if the function that does this is accessible separately? E.g. N[Integrate[...]] does the same as NIntegrate[], but if we use NIntegrate directly, we have control over the parameters of the numerical methods used. Is there any way to control how a numerical derivative is calculated when using N[Derivative[1][f][x]]? There's ND[], but it's part of a package and gives different (and much better) results. $\endgroup$
    – Szabolcs
    Commented Jun 13, 2011 at 14:39
  • $\begingroup$ @Szabolcs Not that I can think of off the top of my head. $\endgroup$
    – Brett Champion
    Commented Jun 13, 2011 at 15:18
  • $\begingroup$ @Szabolcs ND might give better results, but it always takes the derivative from the right (try ND[rf[u], u, x, Terms \[Rule] 2, Scale \[Rule] s] to see this--rf is to be undefined here). N[D[...]] seems to be more careful, but it seems it is using a different (finite difference) method. $\endgroup$
    – acl
    Commented Jun 14, 2011 at 12:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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