Numerical differentiation using 9 data points

Question

I have problems when I'm trying to do numerical differentiation in Matlab. But my question might me more about numerical analysis than about Matlab.

I have an array with 9 data points that represents f(x) for 9 different x. I need to find f''(x) numerically. The values I have for x and f(x) are

x = [2271.38, 2555.30, 2697.26, 2768.24, 2839.22, 2910.20, 2981.18, 3123.14, 3407.06]

f(x) = [577.4063, 311.3341, 193.0833, 141.3048, 95.1501, 58.8130 32.4931, 6.9511, 0.1481]

and I can interpolate to get a smooth curve. I use spline interpolation but is some other interpolation preferable when you are going to differentiate?

I've tried different methods:

Just simple forward,backward and central difference quotients

A wavelet based method: http://www.mathworks.com/matlabcentral/fileexchange/13948-numerical-differentiation-based-on-wavelet-transforms

and the derivest suite: http://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

Non of this have worked satisfactory. The second derivative is very unstable with respect to the step length and the adaptive method in the derivest suite works terribly bad. Maybe I'm just using it in the wrong way!

Any help is appreciated!

Thanks in advance

Answer 1

I imagine it was you who asked a similar question on MATLAB Central the other day. You did not post your data there, so I could not really give a good answer then.

Estimation of the second derivative is a difficult thing to do. It is an ill-posed problem. Differentiation is itself a noise amplifier, so estimating the second derivative is "twice" as bad. It is simply not an easy thing to do, certainly not well.

Using this set of points, I chose to estimate a spline model using my SLM toolbox.

x = [2271.38, 2555.30, 2697.26, 2768.24, 2839.22, 2910.20, 2981.18, 3123.14, 3407.06];
f = [577.4063, 311.3341, 193.0833, 141.3048, 95.1501, 58.8130 32.4931, 6.9511, 0.1481];

First of all, plot the data. What can I learn from that plot? Are there any inferences I might choose to make?

The simple plot tells me, along with your comments, that I expect this function to be a monotonic decreasing function. It appears to be asymptotically linear at each end, like a hyperbolic segment, with positive curvature over the entire domain.

So now I would use this information to build a model for your data, using my SLM toolbox.

slm = slmengine(x,f,'plot','on','decreasing','on','knots',20, ...
    'concaveup','on','endconditions','natural');

slmengine is designed to take information from you, in the form of prescriptions for the shape of the curve. What you will find is that by providing such information, it greatly regularizes the shape of the result, to conform to your knowledge of the process. Here I was just making a few guesses about the curve shape from your comments.

In the above call, I instructed SLM to:

• make a plot of the result when done
• create a monotone decreasing function of x
• use 20 equally spaced knots
• force the curve to have an everywhere positive second derivative
• set the second derivatives at the end to be zero

The plot as generated is a gui itself, allowing you to plot the function and data, but also to plot derivatives of the result. The vertical green lines are the knot locations.

Here we see that the curve fit is a reasonable approximation to what you are looking for.

How about the second derivative plot? Of course, SLM is a piecewise cubic tool. Therefore the second derivatives are only piecewise linear. Is this a problem? Will you ask me to provide a tool for higher order splines? Too bad, but no, I won't. Those higher order derivatives are far too poorly estimated to ask for a highly smooth result. In fact, I would be quite happy with this prediction. Note that the glitches in the second derivative were consistent. If I used more knots or fewer, they were still there. This is a nice way to know if the shape you see is a feature of the curve, or merely an artifact of knot placement.

See that the constraints I placed on the shape of the curve resulted in a fit that was quite reasonable, despite the fact that I used far more knots than I had data points. SLM had no trouble in the estimation.

If I wish to try for a smoother estimate of the second derivative, simply use more knots. SLM is relatively fast. So with 50 knots, we get a very similar result for the second derivative curve.

You can find SLM (here) on MATLAB Central. It does require the optimization toolbox.