Take the 2-minute tour ×
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
1  
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

3 Answers 3

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
3  
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 at 14:48

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

ND

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]

filter

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
1  
@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. Nov 13 '12 at 22:10

Your Answer

 
discard

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.