The Second Derivative.

This demonstration shows some examples of computing numerical second derivatives by the formula:

f''(x) = ((f(x+h) - 2 f(x) + f(x-h))/(h^2))

We first investigate the stability of this formula by trying it out on sin(x) for various values of h.

f[x_] := Sin[x];
fp[x_] = D[f[xx],xx] /. xx->x;
Approx[x_,h_] := (f[x+h] - 2 f[x] + f[x-h])/(h^2);


Manipulate[
h := 10^(-n);
DerivData := Table[{i/10.0,Approx[i/10.0,h]},{i,0,31.40}];
Show[Plot[-f[xx],{xx,0,Pi}],ListLinePlot[DerivData,Filling->Axis,PlotRange->Full],
ListPlot[DerivData,PlotStyle->{PointSize[Large]}]],{n,5,18,Appearance->\!\(\*
StyleBox["\"",
FontSize->16]\)Labeled"}]

We now focus on computing a particular value of the derivative, at x = 1, numerically so that we can see the error as a function of h. Here we see the number of correct digits in the derivative with h = (10^(-n)) as a function of n.

NewData :=Table[{n,Abs[Log[10,Abs[Approx[1,10.0^(-n)] + Sin[1]]]]},{n,1,20}]
Show[ListLinePlot[NewData,Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}],
ListPlot[NewData,PlotStyle->{PointSize[Large]}]]

There are a couple of things worth noting about this plot. We still reach about  8 correct digits in the answer, even though we are doing the computation in double precision. However, we start to lose digits much more quickly than we did during the computation of the first derivative.

Here's the actual data.

N[-Sin[1],8]
Column[Table[{n,Approx[1,10.0^(-n)]},{n,1,20}]]

-0.84147098480789650665250232163029899962`8.

({{{1,-0.8407699926874289`}}, {{2,-0.8414639725740081`}}, {{3,-0.8414709147253063`}}, {{4,-0.841470981782777`}}, {{5,-0.8414713370541448`}}, {{6,-0.8415490526658687`}}, {{7,-0.8437694987151191`}}, {{8,-2.2204460492503126`}}, {{9,0.`}}, {{10,0.`}}, {{11,0.`}}, {{12,0.`}}, {{13,0.`}}, {{14,0.`}}, {{15,0.`}}, {{16,0.`}}, {{17,0.`}}, {{18,0.`}}, {{19,0.`}}, {{20,0.`}}})