# Numerical derivative from data points

I want to take the second derivative of a set of data points (at the point $x = 0$). Let's assume that this set looks like that

{{0, 1}, {5/999, 0.999925}, {5/333, 0.9997}, {25/999, 0.999325}, {35/
999, 0.998801}, {5/111, 0.998127}, {55/999, 0.997305}, {65/999,
0.996335}, {25/333, 0.995217}, {85/999, 0.993952}, {95/999,
0.992542}, {35/333, 0.990986}, {115/999, 0.989287}, {125/999,
0.987444}, {5/37, 0.98546}, {145/999, 0.983336}, {155/999,
0.981073}, {55/333, 0.978672}, {175/999, 0.976136}, {5/27,
0.973465}, {65/333, 0.970662}, {205/999, 0.967728}, {215/999,
0.964665}, {25/111, 0.961475}, {235/999, 0.95816}, {245/999,
0.954722}, {85/333, 0.951163}, {265/999, 0.947486}, {275/999,
0.943692}, {95/333, 0.939783}, {295/999, 0.935763}, {305/999,
0.931632}, {35/111, 0.927395}, {325/999, 0.923052}, {335/999,
0.918607}, {115/333, 0.914062}, {355/999, 0.909419}, {365/999,
0.904682}, {125/333, 0.899851}, {385/999, 0.894931}, {395/999,
0.889923}}


One possibility would to take a function and fitting it to the data set. After that I could take the second derivative at $x=0$. But there the question arises, which function one should take.

Is there any Mathematica routine which can do this "Derivation of numerical data" by its own?

Edit: Thanks for your answers. My problem is probably a bit more advanced, because I'm going to add Gaussian noise to the data. This results in the new data set:

  lst = {{0, 1}, {5/999, 0.9997470723678498}, {5/333,
0.9998479841238097}, {25/999, 0.9953393855010988}, {35/999,
0.9981898503272589}, {5/111, 0.9974638849540485}, {55/999,
0.9971895950455869}, {65/999, 0.9951395169325857}, {25/333,
0.988116112128079}, {85/999, 0.9925548979371955}, {95/999,
0.9852737364775073}, {35/333, 1.000135408043998}, {115/999,
0.9945341860655752}, {125/999, 0.9886336384450369}, {5/37,
0.9931098663739206}, {145/999, 0.9824257991680876}, {155/999,
0.9859582499945159}, {55/333, 0.9850262632503664}, {175/999,
0.9670026311764188}, {5/27, 0.9516452584645868}, {65/333,
0.9492761241829594}}


By using this set, "Interpolation" and a fit with a polynomial fit of degree 4 fails. Do you have any idea how to handle this?

• Are you aware of Interpolation? Possible duplicate of mathematica.stackexchange.com/q/14023/131 and/or mathematica.stackexchange.com/q/10986/131 (amongst others). – Yves Klett Apr 21 '16 at 13:39
• Theres NDSolveFiniteDifferenceDerivative. See for example mathematica.stackexchange.com/q/31832 or The Numerical Method of Lines for more. – Michael E2 Apr 21 '16 at 14:11
• To address your edits: I think that the problem is unclear formulated. It can have several answers. To prefer one of them one needs to know some more information that you probably keep in head, but do not write. (1) Where did the data come from? (2) Why do you add a noise to it? (3) What is expected from the function in question? – Alexei Boulbitch Apr 22 '16 at 7:32
• Continuation: You should understand that taking the same data and few times adding random numbers to it (depending upon the noise amplitude) you may obtain for the derivatives results very much different from one another. Just because numerically the first derivative is (y2-y1)/(x2-x1) and you add to y2 and y1 some different random numbers. The derivative takes the form: (y2+dy2-y1-dy1)/(x2-x1). Then everything depends upon the ratios dy2/y2 and dy1/y1. It is still worse with the second derivative. – Alexei Boulbitch Apr 22 '16 at 7:37
• Thank you. I try to answer your questions. (1) The data is coming from a calculation, where only discrete time steps are taken to simulate a situation in experiment. (so the data without noise is correct in the limit of the theoretical model) (2) To simulate the situation in experiment. The noise are not just random numbers, but Gaussian noise. (3) As said before, I have to take the second derivative of that function at the time t=0. Without noise, the poly(4) (where this shall denote now a polynomial of the degree 4) fit leads to a very good result in comparison to the exact result. – QuantumMechanics Apr 22 '16 at 15:00

First of all there is somewhere in Mma a package for numerical calculation of derivatives, but I did not manage to find a reference.

To offer a way to calculate the derivative. You could use the interpolation function. Here is your list:

lst = {{0, 1}, {5/999, 0.999925}, {5/333, 0.9997}, {25/999,
0.999325}, {35/999, 0.998801}, {5/111, 0.998127}, {55/999,
0.997305}, {65/999, 0.996335}, {25/333, 0.995217}, {85/999,
0.993952}, {95/999, 0.992542}, {35/333, 0.990986}, {115/999,
0.989287}, {125/999, 0.987444}, {5/37, 0.98546}, {145/999,
0.983336}, {155/999, 0.981073}, {55/333, 0.978672}, {175/999,
0.976136}, {5/27, 0.973465}, {65/333, 0.970662}, {205/999,
0.967728}, {215/999, 0.964665}, {25/111, 0.961475}, {235/999,
0.95816}, {245/999, 0.954722}, {85/333, 0.951163}, {265/999,
0.947486}, {275/999, 0.943692}, {95/333, 0.939783}, {295/999,
0.935763}, {305/999, 0.931632}, {35/111, 0.927395}, {325/999,
0.923052}, {335/999, 0.918607}, {115/333, 0.914062}, {355/999,
0.909419}, {365/999, 0.904682}, {125/333, 0.899851}, {385/999,
0.894931}, {395/999, 0.889923}};


This is the interpolation function and its plot:

 f = Interpolation[lst, InterpolationOrder -> 3]
Plot[f[x], {x, 0, 395/999}]


looking as follows:

Now here is its second derivative and the corresponding plot:

g[x_] = D[f[x], {x, 2}];
Plot[g[x], {x, 0, 395/999}]


I do not like this solution because of the fast variation in the vicinity of zero. It is most probably related to the discreteness of the data, but you might be satisfied with this solution.

If, in contrast, your solution is expected to be smooth together with the second derivative, I would still make the fitting. It is not difficult to guess the model:

model = 1 - a*x^2 - b*x^4;
ff = FindFit[lst, model, {a, b}, x]
Show[{
ListPlot[lst],
Plot[model /. ff, {x, 0, 395/999},
PlotStyle -> Directive[Red, Thickness[0.003]]]
}]


(* {a -> 0.786669, b -> -0.543841} *)

Then

 der2 = D[model /. ff, {x, 2}];
Plot[der2, {x, 0, 395/999}]


yielding

Have fun!

• Thank you! I've edited my question to a more realistic but also advanced case. – QuantumMechanics Apr 21 '16 at 14:43

In version 10.2 or later you can use the EXPERIMENTAL function FindFormula

lst = {{0, 1}, {5/999, 0.999925}, {5/333, 0.9997}, {25/999,
0.999325}, {35/999, 0.998801}, {5/111, 0.998127}, {55/999,
0.997305}, {65/999, 0.996335}, {25/333, 0.995217}, {85/999,
0.993952}, {95/999, 0.992542}, {35/333, 0.990986}, {115/999,
0.989287}, {125/999, 0.987444}, {5/37, 0.98546}, {145/999,
0.983336}, {155/999, 0.981073}, {55/333, 0.978672}, {175/999,
0.976136}, {5/27, 0.973465}, {65/333, 0.970662}, {205/999,
0.967728}, {215/999, 0.964665}, {25/111, 0.961475}, {235/999,
0.95816}, {245/999, 0.954722}, {85/333, 0.951163}, {265/999,
0.947486}, {275/999, 0.943692}, {95/333, 0.939783}, {295/999,
0.935763}, {305/999, 0.931632}, {35/111, 0.927395}, {325/999,
0.923052}, {335/999, 0.918607}, {115/333, 0.914062}, {355/999,
0.909419}, {365/999, 0.904682}, {125/333, 0.899851}, {385/999,
0.894931}, {395/999, 0.889923}};

f[x_] = FindFormula[lst][x]

(*  0.999983 - 0.0071517 x - 0.756732 x^2 + 0.0662479 x^3 + 0.285862 x^4  *)

{xmin, xmax} = MinMax[lst[[All, 1]]];

Plot[f''[x], {x, xmin, max}]
`

• Thank you! I've edited my question to a more realistic but also advanced case. – QuantumMechanics Apr 21 '16 at 14:43