Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
Are you aware of ND? –  VLC Nov 11 '12 at 17:43
Have a look at "NumericalCalculus`". –  b.gatessucks Nov 11 '12 at 17:43

4 Answers 4

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]


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[l, pad, "Extrapolated"], {1}],



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

share|improve this answer
Something like N[D[f[x], x] /. x -> 5] also works. There is a difference between ND and D, see here –  acl Nov 11 '12 at 22:58
@acl, yes - and I had upvoted your linked answer back then, too! It's very exhaustive. –  Jens Nov 11 '12 at 23:00
Lettuce knot forget scenario 3, of Automatic Differentiation. –  Daniel Lichtblau Nov 12 '12 at 15:56
@DanielLichtblau If I look at this Wolfram demo, it seems to me that the technique is better classified as a symbolic differentiation trick. But maybe that's a matter of taste. My criterion would be that a numerical differentiation should be able to deal with functions that are calculated by numerical root finding, or using FixedPoint etc. Automatic differentiation requires that the function be calculable using elementary operations. –  Jens Nov 12 '12 at 19:20

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

plot of exact and approximate derivatives

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.

share|improve this answer
For the discrete case, see this. There's also the possibility of using a Savitzky-Golay filter for differentiating... –  J. M. is back. Nov 13 '12 at 22:10

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.

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.

share|improve this answer
I think the difference between this and ND is that ND takes the derivative from one side while N[D[]] takes it symmetrically. –  acl Nov 12 '12 at 20:32
@acl Sorry about this, it seems you already warned me once before ... –  Szabolcs Jun 16 '14 at 14:48

There is another possibility, and that is to differentiate an InterpolatingFunction.

Here is an example of a rather poorly behaving function whose derivative we can recover by this technique. The transfer function of a Butterworth filter looks like

h[s_] := 1 / (1 + 2s + 2s^2 + s^3)

for s complex. The Arg of this function, which, for purely imaginary frequencies, is the phase of the filter, suffers a branch cut at an inconvenient place:

Plot[Arg[h[I s]], {s, 0, 6}]

enter image description here

The symbolic derivative of this thing, which is group delay for the filter, bogusly produces Complex values, even if we didn't have the problem of the inconvenient branch cut:

Table[D[Arg[h[I s]], s] /. s -> t, {t, 0, 2}]
{ -2 I Arg', 
  (-(1/2) + 2 I) Arg'[-(1/2) - I/2], 
  (824/4225 - (118 I)/4225) Arg'[-(7/65) + (4 I)/65] }

We can get a better-behaved function by adding 2 Pi and then modding it out:

Plot[Mod[2 Pi + Arg[h[I s]], 2 Pi], {s, 0, 6}, PlotRange -> All]

enter image description here

But now we have a "dimple" where we had the branch cut, and, sure enough, ND from NumericalCalculus has problems:

<< NumericalCalculus`
Plot[ND[Mod[2 Pi + Arg[h[I w]], 2 Pi], w, s], {s, 0, 6}, PlotRange -> Full]

enter image description here

The DerivativeFilter method above works great to fix this, but here is another way. First, sample the original Arg or phase function:

ListPlot[args$ = Table[{s, Mod[2 Pi + Arg[h[I s]], 2 Pi]}, {s, .05, 6, .05}]]

enter image description here

Build an interpolation function from it:

iargs$ = Interpolation[args$]
InterpolatingFunction[{{0.05, 6.}}, <>]

and differentiate that (notice carefully the tick mark in the following expression, a convenient shorthand for derivatives of functions of single variables):

Plot[iargs$'[s], {s, 0, 6}]

enter image description here

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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