Numerical differentiation methods

Question

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

Answer 1

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_}] :=
NDSolve`FiniteDifferenceDerivative[{2, 0}, {xgrid, ygrid}, vals] +
NDSolve`FiniteDifferenceDerivative[{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]

Answer 2

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

Answer 3

Short answer: yes.

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.

Answer 4

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)

gradient[l_,dt_]:=     Flatten[{(l[[2]]-l[[1]])/dt,(l[[3;;]]- 
                   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.

Answer 5

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)

Answer 6

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]