How does D work with numerical arguments?

Question

If I define

f[x_?NumericQ] := (Sow[x]; SawtoothWave[x])
D[f[x], x] /. x -> 1.1
(* 5.29432 *)

numeric differentiation fails, despite the function being smooth and well-defined in a region around 1.1.

We can see the values of x at which f is evaluated

Reap[D[f[x], x] /. x -> 1.1]
(* {5.29432, {{1.1, 1.15263, 1.20526, 1.25789, 1.31053, 1.36316, 1.41579,
    1.46842, 1.52105, 1.04737, 0.994737, 0.942105, 0.889474, 0.836842,
    0.784211, 0.731579, 0.678947}}} *)

Mathematica provides ND for numerical differentiation, which behaves sensibly

Needs["NumericalCalculus`"]
Reap[ND[f[x], x, 1.1, Scale -> 0.01]]
(* {1., {{1.1, 1.11, 1.105, 1.1025, 1.10125, 1.10063, 1.10031, 1.10016}}} *)

The documentation for D does not provide obvious clues on what is to be expected with numeric arguments.

The sensible advice is not to use D with numeric functions. However, it may not always be obvious to the user that D will be applied to an expression.

Related question Derivative of mod gives unacceptable results

D[f[x],x] gets turned into Derivative[1][f][x] (f'), so this is really a question of how Derivative behaves for numerical functions. You can investigate what's happening with: f[x_?NumericQ] := (Echo[{x, SawtoothWave[x]}]; SawtoothWave[x]); D[f[x], x] /. x -> 1.1. It seems to naively sample some points with distance O(1), and then generate a result. — QuantumDot, Aug 29 '16 at 11:26
Somewhere on this site someone discusses the numerical method used by Derivative, which, as I recall, is slightly different from ND. It could be discussed as a side-issue in a Q&A about something else, but for whatever reason, I can't locate it... — Michael E2, Aug 29 '16 at 13:05
I think this is the one I was thinking of: (29329) --The question is somewhat different, but @acl's answer shows that the eighth-order central difference is being used to numerically compute the derivative. That seems to be true here, too. — Michael E2, Aug 30 '16 at 00:58

mikado · Answer 1 · 2016-08-29T13:16:30.420

Here is a partial explanation

f[x_?NumericQ] := (Sow[x]; SawtoothWave[x])

Collecting 'x' and 'f[x]' and defining an interpolation function of the maximum supported order

data = {#, f[#]} & /@ Flatten[Last[Reap[D[f[x], x] /. x -> 1.1]]]
(* {{1.1, 0.1}, {1.15263, 0.152632}, {1.20526, 0.205263}, {1.25789, 
  0.257895}, {1.31053, 0.310526}, {1.36316, 0.363158}, {1.41579, 
  0.415789}, {1.46842, 0.468421}, {1.52105, 0.521053}, {1.04737, 
  0.0473684}, {0.994737, 0.994737}, {0.942105, 0.942105}, {0.889474, 
  0.889474}, {0.836842, 0.836842}, {0.784211, 0.784211}, {0.731579, 
  0.731579}, {0.678947, 0.678947}} *)

intdata = Interpolation[data, InterpolationOrder -> 16];

The interpolating function passes through the points, but the order is much too high to give a useful result here.

    Show[{ListPlot[data], Plot[intdata[t], {t, 0.8, 1.4}]}, PlotRange -> {-1, 2}]

We have a match to the value returned by D.

{intdata'[1.1], D[f[x], x] /. x -> 1.1}
(* {5.29432, 5.29432} *)

This doesn't explain the sample values chosen, those these are obviously equally spaced.

EDIT

It seems that the sample values have a spacing of 1/19, independent of where the derivative is evaluated.

How does D work with numerical arguments?

1 Answers1