Help with numerical differentiation

Hi! I have problems when I'm trying to do numerical differentiation.

I have an array with 9 data points that represents f(x) for 9 different x. Note that these are actual measurements with error and there is no "real" function f. 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?

Here is the interpolated curve f(x): http://i.imgur.com/j3AZl.png

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

[Joakim]

Thank you for your reply! I got help from John D'Errico and his amazing SLM toolbox:

http://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling

It did the trick and I got the result I was expecting!

Answer 2

[Jan]

What do you expect as 2nd derivative of noisy data? Any level of smoothness will be purely injected by the computational method. Therefore the results will be too artificial to be reliable.

Do you have any information about the physical object behind the measurement? Then a good strategy is to fit the model parameters to the measurement data. E.g. if the values are the position of a pendulum, you adjust the string length until the trajetory matchs your data in a best way. Then the 2nd derivative will be easy to calculate, smooth and reliabale. But imagine, the values are stock prices - then a 2nd derivative would be even meaningless.

Without considering the physical model, smooth 2nd derivatives are as magic as predicting the value of a die based on any combination of the previous values.

Answer 3

[Teja Muppirala]

If you have the Curve Fitting Toolbox, there are some functions in there that may be of use, like FIT and DIFFERENTIATE:

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]

xvec = linspace(x(1),x(end),1001);

F = fit(x(:),f(:), 'smoothingspline'); 

figure;

subplot(211);
plot(x,f,'ro',xvec,F(xvec),'b');
title('Interpolated function'); 
grid 

subplot(212);
[~,deriv2] = differentiate(F,xvec);
plot(xvec,deriv2);
title('Second Derivative Estimate')
grid

[Joakim]

This did actually work pretty well. At least better than the methods I came up with my self! However, I found out that the SLM tool did the trick! I really recommend it!

http://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling