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

How can I differentiate a function numerically?

  • 2
    $\begingroup$ Are you aware of ND? $\endgroup$
    – VLC
    Commented Nov 11, 2012 at 17:43
  • 1
    $\begingroup$ Have a look at "NumericalCalculus`". $\endgroup$ Commented Nov 11, 2012 at 17:43
  • $\begingroup$ For continuously defined functions e.g. f[x_] := Sin[x]/x + E^(-x); Plot[f'[x], {x, -1, 9}] $\endgroup$
    – Narasimham
    Commented Mar 10, 2016 at 10:43

5 Answers 5


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

  • 1
    $\begingroup$ Something like N[D[f[x], x] /. x -> 5] also works. There is a difference between ND and D, see here $\endgroup$
    – acl
    Commented Nov 11, 2012 at 22:58
  • 1
    $\begingroup$ @acl, yes - and I had upvoted your linked answer back then, too! It's very exhaustive. $\endgroup$
    – Jens
    Commented Nov 11, 2012 at 23:00
  • 2
    $\begingroup$ Lettuce knot forget scenario 3, of Automatic Differentiation. $\endgroup$ Commented Nov 12, 2012 at 15:56
  • $\begingroup$ @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. $\endgroup$
    – Jens
    Commented Nov 12, 2012 at 19:20

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.

  • 3
    $\begingroup$ I think the difference between this and ND is that ND takes the derivative from one side while N[D[]] takes it symmetrically. $\endgroup$
    – acl
    Commented Nov 12, 2012 at 20:32
  • $\begingroup$ @acl Sorry about this, it seems you already warned me once before ... $\endgroup$
    – Szabolcs
    Commented Jun 16, 2014 at 14:48

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.

  • $\begingroup$ For the discrete case, see this. There's also the possibility of using a Savitzky-Golay filter for differentiating... $\endgroup$ Commented Nov 13, 2012 at 22:10
  • $\begingroup$ can we also use acegen? $\endgroup$
    Commented Mar 19, 2019 at 1:10

I present here another approach to numerically differentiating a function $f(x)$. While the methods in my other answer approximate $f^\prime(x)$ at a single specified value, the following method constructs an approximate function $p(x)$ such that $p(x)\approx f(x)$ within a given interval $[a,b]$, and then proceeds to directly generate $p^\prime(x)$.

The key here is to replace the function with its Chebyshev series, and then differentiate that series. This is similar to what is done in the Chebfun package for MATLAB, and previously presented in a number of answers by Michael E2 and yours truly.

To recall, this is how one might derive the coefficients of the Chebyshev series of a function over a given interval (as previously shown here):

f[x_] := BesselJ[0, x];
{a, b} = {0, 19}; (* interval of approximation *)

n = 128; (* arbitrarily chosen integer *)
prec = 20; (* precision *)
cnodes = Rescale[N[Cos[π Range[0, n]/n], prec], {-1, 1}, {a, b}];
fc = f /@ cnodes;
cc = Sqrt[2/n] FourierDCT[fc, 1];
cc[[{1, -1}]] /= 2;

Then, one can either take the dot product of the Chebyshev coefficients cc with an appropriate list of suitably scaled Chebyshev polynomials (of the first kind), or use a more efficient evaluation method such as Clenshaw's algorithm, which was presented in these two answers. For simplicity, I will use the dot product method here:

Plot[{f[x], cc.ChebyshevT[Range[0, n], Rescale[x, {a, b}, {-1, 1}]]},
     {x, a, b}, PlotStyle -> {AbsoluteThickness[4], AbsoluteThickness[1]}]

approximating a Bessel function

Now, one might be tempted to use D[] on the constructed Chebyshev series, but the problem with this approach is that ChebyshevT[] automatically converts to its power basis form for non-numeric arguments, thus negating any stability gained from remaining in the Chebyshev basis. Instead, one should directly construct the Chebyshev coefficients of the derivative.

The way to do that is to use a Chebyshev spectral differentiation matrix. What follows is a routine (previously posted here) to generate it.

chebmat[n_Integer?Positive, prec_: MachinePrecision] := 
 Module[{fac, xm},
        xm = ConstantArray[N[-Cos[π Range[0, n]/n], prec], n + 1];
        fac = Flatten[{2, PadRight[{}, n - 1, {-1, 1}], 2 - 4 Mod[n, 2]}];
        xm = Outer[Times, fac, 1/fac]/(xm - Transpose[xm] + IdentityMatrix[n + 1]);
        xm - DiagonalMatrix[Total[xm, {2}]]]

Now to generate the Chebyshev coefficients of the first derivative, do this:

dc = chebmat[n, prec].fc;
cd = FourierDCT[dc, 1] Sqrt[2/n] (2/(b - a));
cd[[{1, -1}]] /= 2;


Plot[{f'[x], cd.ChebyshevT[Range[0, n], Rescale[x, {a, b}, {-1, 1}]]},
     {x, a, b}, PlotStyle -> {AbsoluteThickness[4], AbsoluteThickness[1]}]

approximating the derivative, Chebyshev-style

This approach is easily extensible to arbitrary order derivatives, just by multiplying fc with an appropriate power of chebmat[n]. For instance, here's how to generate the Chebyshev coefficients of the third derivative:

k = 3; (* derivative order *)
dck = MatrixPower[chebmat[n, prec], k, fc];
cdk = FourierDCT[dck, 1] Sqrt[2/n] (2/(b - a))^k;
cdk[[{1, -1}]] /= 2;
  • $\begingroup$ I think you've got a typo: dc should read dc = chebmat[n, prec].(f /@ cnodes). Might oughtta check dck, too. (I guess I didn't check this when I upvoted...Probably too busy.) $\endgroup$
    – Michael E2
    Commented Mar 9, 2016 at 19:09
  • $\begingroup$ @Michael, whoops. I must've mixed up lines during pasting. Should be fixed now. $\endgroup$ Commented Mar 9, 2016 at 23:22

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


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.