# Numerical differentiation methods

Is it possible to write code in Mathematica that implements various differentiation methods (like forward, central, extrapolated, etc.)?

• Long story - short, yeah, you can, you having troubles implementing it/them ? Sep 7, 2013 at 21:50
• See this, that, or the other for various approaches. Sep 9, 2013 at 16:01
• See this as well. May 5, 2015 at 7:09

In my experience almost all finite difference formulas can be implemented very efficiently using ListCorrelate. Let's look at how to implement a few difference formulas in 1D with periodic boundary conditions on a uniform mesh with spacing h:

Central Difference

CDx[v_List, h_]:=ListCorrelate[ {-1, 0, 1}/(2 h), v, 2];


Backward Difference

BDx[v_List, h_] := ListCorrelate[{-1, 1}/(h), v, 2];


Forward Difference

FDx[v_List, h_] := ListCorrelate[{-1, 1}/(h), v];


We can check that they are performing the correct behavior via a simple symbolic test

v = Table[x[i], {i, 1, 5}];
CDx[v, h]
BDx[v, h]
FDx[v, h]


FiniteDifferenceDerivative Function

With that said, there is also a very nice function hidden inside of the NDSolve context called FiniteDifferenceDerivative. Search for "tutorial/NDSolveMethodOfLines" inside of the Documentation Center and look for the section "Spatial Derivative Approximations". This function allows for differentiation on a nonuniform mesh. The user must supply a grid and values. For example the 2D Lapacian operator is implemented via:

lap[vals_, {xgrid_, ygrid_}] :=
NDSolveFiniteDifferenceDerivative[{2, 0}, {xgrid, ygrid}, vals] +
NDSolveFiniteDifferenceDerivative[{0, 2}, {xgrid, ygrid}, vals]


Let's test it out:

n = 16;
h = 1.0/n;
xgrid = Table[x, {x, 0, 1 - h, h}];
ygrid = Table[y, {y, 0, 1 - h, h}];
vals = Table[Sin[2 Pi x]*Sin[2 Pi*y], {x, xgrid}, {y, ygrid}];
valsDiff = lap[vals, {xgrid, ygrid}];
(* Plot to see solution *)
ListPlot3D[vals]
ListPlot3D[valsDiff]

• You probably want to remove the Head Real from the definition of CDx otherwise your symbolic test CDx[v,h] will not evaluate. Sep 7, 2013 at 23:05
• @RunnKine Good catch. I've made the edit. Sep 8, 2013 at 1:37

To find the finite-difference formula using m grid points, a polynomial of order m-1 is differentiated. This method is recommended over using Taylor series expansion when the number of grid points becomes large. The Fornberg formula (implemented as function FDFormula in the above link, down the page) can also be used to generate the finite-difference formula efficiently for larger numbers of grid points.

I have a demo that uses this method in the pipeline for the demonstration web site, but if you like you can try it from here

A report on this using Mathematica is here. Examples of using polynomial method:

## Examples

order = 2;
nPoints = 3;
points = makePoints[nPoints, "centered"];
mkFormula[order, points, x, h, f]


order = 1;
nPoints = 2;
points = makePoints[nPoints, "forward"];
mkFormula[order, points, x, h, f]


order = 1;
nPoints = 2;
points = makePoints[nPoints, "backward"];
mkFormula[order, points, x, h, f]


order = 1;
nPoints = 10;
points = makePoints[nPoints, "backward"];
mkFormula[order, points, x, h, f]


## Code

mkFormula[order_Integer /; order >= 1,
grid_List,(*list of points*)
x_,(*symbol to use,has no value*)
h_,(*symbol to use for spacing,has no value*)
f_(*symbol to use for f[
x] has no value*)] /; (order <= Length[grid] - 1) :=
Module[{points, e, z, p},
points = {x + # h, f[x + # h]} & /@ grid;
p = InterpolatingPolynomial[points, z];
e = D[p, {z, order}];
e = e /. z -> x;
Simplify[Together[Numerator[e]]/Denominator[e]]];

makePoints[nPoints_Integer /; nPoints >= 2,
finiteDifferenceType_String] := Module[{d = nPoints - 1},
Which[finiteDifferenceType == "centered", Table[n, {n, -d/2, d/2}],
finiteDifferenceType == "forward", Table[n, {n, 0, nPoints - 1}],
finiteDifferenceType == "backward", Table[n, {n, -nPoints + 1, 0}]
]
]


Slightly longer answer: this is the simplest way that I could think of, several years ago. I hope the code still works on newer versions of mma.

forwardDiff[f_, h_:1] := Drop[RotateLeft[f] - f, -1]/h
centralDiff[f_, h_:1] :=
Rest[Drop[RotateLeft[f] - RotateRight[f], -1]]/(2h)
backwardDiff[f_, h_:1] := Drop[f - RotateRight[f], 1]/h


You can test it.

f[x_] := 3Sin[x]Exp[-.5x];
data = Table[f[x], {x, 0, 7, .2}];
ListPlot[data, Frame -> True]

d1 = forwardDiff[data, 0.2];
ListPlot[d1, Frame -> True, PlotRange -> All]


And so on with the other schemes. Edit: please note that with every differentiation you lose a point in your data. These schemes are not indicated for implementation of multiple differentiations (I am not aware if they can be useful at all in the real world, perhaps this form can be used in the microcontroller's code of an embedded system...) More efficient methods compute an approximation of the (nth) derivative using Taylor's polynomial, the (easily computable) derivatives of a suitable interpolating polynomial or, if you really want to show off :-) a Cauchy integral.

• Alternative: forwardDiff[f_, h_:1] := (f[[2;;]]-f[[;;-2]])/h Sep 7, 2013 at 22:09

You can implement them as :

ForwardDiff[l_,dt_]:=  (l[[2;;]]-l[[;;-2]])/dt

BackwardDiff[l_,dt_]:= (l[[;;-2]]-l[[2;;]])/-dt

CentralDiff[l_,dt_]:=  (l[[3;;]]-l[[;;-3]])/(2 dt)

l[[1;;-3]])/(2 dt),(l[[-1]]-l[[-2]])/dt}]


Let us take some sample condition:

position = Table[t^2 + 2 t, {t, 0, 1, 0.1}] (*Length = 11*)
{0., 0.21, 0.44, 0.69, 0.96, 1.25, 1.56, 1.89, 2.24, 2.61, 3.}

velocity = ForwardDiff[position, 0.1] (*Length = 10*)
{2.1, 2.3, 2.5, 2.7, 2.9, 3.1, 3.3, 3.5, 3.7, 3.9}

velocity = BackwardDiff[position, 0.1] (*Length = 10*)
{2.1, 2.3, 2.5, 2.7, 2.9, 3.1, 3.3, 3.5, 3.7, 3.9}

velocity = CentralDiff[position, 0.1] (*Length = 9*)
{2.2, 2.4, 2.6, 2.8, 3., 3.2, 3.4, 3.6, 3.8}


I have added one more function gradient, which will keep the length equal to 11.

velocity = gradient[position, 0.1] (*Length = 11*)
{2.1, 2.2, 2.4, 2.6, 2.8, 3., 3.2, 3.4, 3.6, 3.8, 3.9}


Actual Velocity Vector:

velocity= Table[2 t + 2, {t, 0, 1, 0.1}] (*Length = 11*)
{2., 2.2, 2.4, 2.6, 2.8, 3., 3.2, 3.4, 3.6, 3.8, 4.}


I hope it will work.

• Hi! Welcome to Mathematica.SE! You might find interesting that in Span (;;), when the range starts/stops at the beginning/end of the list resp., the 1/Length[l] may be omitted, as in Niki's comment. You can see the equivalence of three methods here: Block[{l = Range@4}, { {l[[2 ;; Length[l]]], Rest[l], l[[2 ;;]]}, {l[[1 ;; Length[l] - 1]], Most[l], l[[;; -2]]} }] Jun 6, 2021 at 16:16
• Thanks @MichaelE2 for pointing it out. Now I have edited it. I have also seen the NIKI's implementation, I just wanted vectorize it without using any other function. Jun 7, 2021 at 10:13

This is an old question, but it seems the function ND hasn't been mentioned yet. If you're unsure of what formulas it's using exactly, you can evaluate it with symbolic functions as well to get the formula:

Needs["NumericalCalculus"]
ND[f[x], {x, 1}, x0, Terms -> 2, Scale -> h, WorkingPrecision -> Infinity] // Simplify


-((3 f[x0] - 4 f[h/2 + x0] + f[h + x0])/h)

In the spirit of sharing ancient code I wrote on this topic that is still working and relevant today, here is some code in a notebook I called "difference differentials". You should like this, I certainly got a kick out of it when I did it. It was a 4 beer late night curiosity event that was never shared until now. This was from the early undergrad days, a private curiosity, and I was basically "hacking" the Mathematica with few skills, so don't pick on me, I actually have grown by light years since this.

My original description:

"This code uses central differences to approximate the derivatives of a function. Change lo, hi, del, and f[x_] as you please, then go to Evaluation->Evaluate Notebook for a new result."

 lo = -12;
hi = 12;
del = 1;
f[x_] := (x + 16) (x + 5) (x) (x - 11) (x - 14)/10000

plot = Plot[f[x], {x, lo, hi}, AxesLabel -> {"x", "f[x]"},
PlotStyle -> Directive[Opacity[.75], Thick, Blue]];
plottable = Table[{x, (f[x])}, {x, lo, hi, del}];
plot1prime =
Plot[f'[x], {x, lo, hi}, AxesLabel -> {"x", "f[x]"},
PlotStyle -> Directive[Opacity[.75], Thick, Purple]];
plot1primetable = Table[{x, (f'[x])}, {x, lo, hi, del}];
plot2prime =
Plot[f''[x], {x, lo, hi}, AxesLabel -> {"x", "f[x]"},
PlotStyle -> Directive[Opacity[.75], Thick, Brown]];
plot2primetable = Table[{x, (f''[x])}, {x, lo, hi, del}];
dif1plottable =
Table[{x, (f[x + del] - f[x - del])/(2*del)}, {x, lo, hi,
del}];(*centered difference*)
p1plot = ListPlot[dif1plottable, AxesLabel -> {"t", "P(t)"},
PlotStyle -> Directive[PointSize[Medium], Red]];
dif2plottable =
Table[{x, (f[x + del] - 2*f[x] + f[x - del])/(del^2)}, {x, lo, hi,
del}];(*centered second difference*)
p2plot = ListPlot[dif2plottable, AxesLabel -> {"t", "P(t)"},
PlotStyle -> Directive[PointSize[Medium], Black]];
Show[plot, plot1prime, plot2prime, p1plot, p2plot, AxesOrigin -> {0, 0},
PlotRange -> Automatic]


I am entertained by the result of that output even today, but there was alot I did not know about writing good code at that time that could be improved.

Here is the table output (WOW there was so much I did not know back then!):

 table = Table[{i}, {i, lo, hi, del}];
Table[AppendTo[table[[i]],
plottable[[i]][[2]]], {i, ((hi - lo)/del) + 1}];
Table[AppendTo[table[[i]],
plot1primetable[[i]][[2]]], {i, ((hi - lo)/del) + 1}];
Table[AppendTo[table[[i]],
dif1plottable[[i]][[2]]], {i, ((hi - lo)/del) + 1}];
Table[AppendTo[table[[i]],
plot2primetable[[i]][[2]]], {i, ((hi - lo)/del) + 1}];
Table[AppendTo[table[[i]],
dif2plottable[[i]][[2]]], {i, ((hi - lo)/del) + 1}];
TableForm[N[table],
TableHeadings -> {None, {"x", "f[x]", "f'[x]", "1st differences",
"f''[x]", "2nd differences"}}, TableAlignments -> Center]
`