How can I differentiate Numerically?

Question

Mathematica has two ways to integrate: Integrate and NIntegrate.
But what about D? D and Derivative are for symbolic differentiation.

How can I differentiate numerically?

Answer 1

As others mentioned, there's the ND function from the NumericalCalculus` package.

It's a bit less widely known that Derivative is also able to approximate derivatives purely numerically.

Let's create a numerical black box function:

f[x_?NumericQ] := x^2

_?NumericQ makes sure that the innards of f are inaccessible to Derivative, so f'[x] returns unevaluated. But plugging in an explicit numerical value,

In[]:= f'[3.]
Out[]= 6.

Caveat: In my experience, this method doesn't always give terribly reliable answers (generally worse than ND, at least for the particular functions I tried it for).

---

Update: Please also see a discussion by @acl concluding that the difference between Derivative and ND is that ND[...] takes the derivative from one side while N[D[...]] takes it symmetrically.

Answer 2

There are two rather different scenarios for numerical derivatives:

1. Differentiating a continuous function that's only defined numerically
2. Approximating the derivative of a list of data that could itself be generated numerically

For scenario 1, here is an example function and its derivative:

f[x_?NumericQ] := BesselJ[1, x]

Needs["NumericalCalculus`"]

Plot[ND[f[x], x, y], {y, 0, 10}]

For scenario 2, here I discretize the above function to get a table of values in some interval:

l = Table[f[x], {x, 0, 10, .1}];

Now the derivative can be taken using DerivativeFilter. This doesn't require the "NumericalCalculus" package:

dl = With[{pad = 10},
   ArrayPad[
    DerivativeFilter[
     ArrayPad[l, pad, "Extrapolated"], {1}],
    -pad]
   ];

ListPlot[dl]

Here, I took the list l and padded it (optional) with pad extra entries at the start and end, before taking the numerical derivative using DerivativeFilter. The purpose of the padding is that it allows me to extrapolate the data points so that the derivative at the interval boundaries will look smooth.

After doing the derivative, I remove the padding by using a negative argument in ArrayPad[..., -pad].

Answer 3

Had ND[] not been implemented in the NumericalCalculus` package, one could implement any number of numerical differentiation methods in Mathematica. I gave a number of warnings on the use of, as well as methods for, numerical differentiation in this MO answer; in particular, the last two methods I described lend themselves to somewhat compact implementations:

(* Cauchy method *)
ND1[f_, {x_, x0_}, opts___] :=
    Chop[NIntegrate[Exp[-I t] Function[x, f][x0 + Exp[I t]], {t, -Pi, Pi}, 
                    AccuracyGoal -> Infinity, opts, Method -> "DoubleExponential"]/(2 Pi)]

(* Lanczos method *)
ND2[f_, {x_, x0_}, opts___] := Module[{prec = Precision[x0], h, pr2}, 
    If[prec === Infinity, prec = MachinePrecision]; 
    h = # (Abs[x0] + #) &[(10^-(prec/2))];
    If[prec === MachinePrecision, pr2 = $MachinePrecision, pr2 = prec];
    N[3 NIntegrate[t Function[x, f][t + x0], {t, -h, h}, opts, 
                   WorkingPrecision -> pr2]/(2 h^3), prec]]

Both approaches are easily generalized; the Lanczos method can be generalized to arbitrary integer-order derivatives, and the Cauchy method can be generalized to arbitrary complex-order derivatives. I won't be discussing these generalizations further in this answer, tho.

A demonstration:

Plot[{-BesselJ[1, x], ND1[BesselJ[0, t], {t, x}]}, {x, 0, 10}, Axes -> None,
     Frame -> True, PlotStyle -> {Directive[Gray, Thick, Dashed], Blue}]

It's not too hard to implement the Richardsonian method I alluded to in my MO answer, but the routine is somewhat longer; I'll edit this answer to include it if there's interest.