# Symbolic derivatives are being calculated numerically

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]


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}]


Or similarly

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


What is happening here?

-

## migrated from stackoverflow.comJul 28 '13 at 6:31

This question came from our site for professional and enthusiast programmers.

Strictly speaking IntegerPart[x] is not a differentiable function. – ndroock1 Jun 13 '11 at 8:19
@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. – Leonid Shifrin Jun 13 '11 at 9:08
@Leonid Indeed. Just try Integrate[Derivative[1][IntegerPart][x], {x, .9, 1.1}]. But that doesn't seem to explain the plot – Dr. belisarius Jun 13 '11 at 11:52
@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. – Leonid Shifrin Jun 13 '11 at 13:12
@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] – Szabolcs Jun 13 '11 at 14:23

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:

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

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:

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

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

-
@acl Very nice answer. Hope you don't mind I added a new chapter (and please check if you agree on my update!) – Dr. belisarius Jun 14 '11 at 17:24
@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! – acl Jun 14 '11 at 17:56
@acl Yep. It's right. See update – Dr. belisarius Jun 14 '11 at 17:59
@belisarius wow, it didn't occur to me at all to check the results by looking them up! brilliant! – acl Jun 14 '11 at 18:04
@belisarius actually it feels almost exactly like doing physics but without the paper- and grant-writing, political manoeuvring etc, of course... hmmm... – acl Jun 14 '11 at 19:01

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[%]


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

-
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. – Szabolcs Jun 13 '11 at 14:39
@Szabolcs Not that I can think of off the top of my head. – Brett Champion Jun 13 '11 at 15:18
@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. – acl Jun 14 '11 at 12:07