Accurate numerical differentiation in MATLAB?
in [Matlab]
|
From: Evan Ruzanski on 2 Mar 2010 00:30 Hello, I'm trying to compute the derivative of a time series. Previous research has shown that simple differences [i.e., using y = diff(x)] produce large errors for large values of derivatives for this particular quantity and that computing the derivative using higher-order polynomials from piecewise Lagrangian interpolation polynomials is about 2 orders of magnitude more accurate. Does anyone know of a MATLAB function/package to perform this (or similar) type of accurate numerical differentiation??? Many thanks!
From: Matt Fig on 2 Mar 2010 00:41 I am not sure if this is all you will need, but John D'Errico does top-notch work. You might want to check out some of his tools, especially this one: http://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling
From: Mark Shore on 2 Mar 2010 07:31 "Matt Fig" <spamanon(a)yahoo.com> wrote in message <hmi8df$bgl$1(a)fred.mathworks.com>... > I am not sure if this is all you will need, but John D'Errico does top-notch work. You might want to check out some of his tools, especially this one: > > http://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling You might also want to look at http://www.mathworks.com/matlabcentral/fileexchange/13948-numerical-differentiation-based-on-wavelet-transforms.
From: Mark Shore on 2 Mar 2010 08:12 > You might also want to look at http://www.mathworks.com/matlabcentral/fileexchange/13948-numerical-differentiation-based-on-wavelet-transforms. I broke the link by putting a period at the end. http://www.mathworks.com/matlabcentral/fileexchange/13948-numerical-differentiation-based-on-wavelet-transforms
From: John D'Errico on 2 Mar 2010 11:14
"Evan Ruzanski" <ruzanski.02(a)engr.colostate.edu> wrote in message <hmi7pe$2or$1(a)fred.mathworks.com>... > Hello, > > I'm trying to compute the derivative of a time series. Previous research has shown that simple differences [i.e., using y = diff(x)] produce large errors for large values of derivatives for this particular quantity and that computing the derivative using higher-order polynomials from piecewise Lagrangian interpolation polynomials is about 2 orders of magnitude more accurate. > > Does anyone know of a MATLAB function/package to perform this (or similar) type of accurate numerical differentiation??? > Numerical differentiation from a data series is a common problem. While you can use high order Lagrange interpolating polynomials, I never recommend such a thing. Yes, I showed how to do interpolating polynomials, in one of my blogs, but I will never recommend anything over about 3rd or 5th order polynomials. Even then, a Lagrange polynomial has poor behavior in terms of continuity of the derivatives. Worse, those polynomials can do seriously nasty things between the data points if the polynomial has a high enough order. One solution is a cubic interpolating spline. These are easily fit for data series, since they involve a simply solved linear system of equations. (Usually tridiagonal.) Differentiate the cubic spline, and you have a derivative estimator. Since the cubic spline will be a C2 function, the derivative estimates will generally be decent. Of course I can always find a way to trip up any such scheme. For example, compute the derivatives of a simple series... t = 0:20; y = double(t >= 10); This function has zero derivative everywhere, except for a derivative singularity between 9 and 10 in t. If you have the splines toolbox, then you can easily plot these things. f = spline(t,y); fp = fnder(f); fnplt(f) fnplt(fp) If you don't have the splines toolbox, then it is still possible to differentiate a cubic spline as a function. Just do this: fp = f; fp.coefs = f.coefs*[3 0 0;0 2 0;0 0 1;0 0 0]; fp.order = 3; And we can plot these curves as easily, by evaluating the curves at a set of points. that = linspace(t(1),t(end),1000); yhat = ppval(that,f); yphat = ppval(that,fp); plot(t,y,'o',that,yhat,'r-') Importantly, see that the spline interpolant shows ringing behavior near the derivative singularity. This is logical of course, since a cubic spline interpolant simply cannot well approximate a function with a singularity in it. Lagrange interpolating polynomials would do as poorly of course, especially high order polynomials of that form. We can do better for this class of underlying function, by using a spline form that represents curves with such singularities. Use the pchip function to build the spline. f = pchip(t,y); See that this curve is a far better approximation to the original step function than was the simple cubic spline. There are other methods one might employ for these problems. For example, Savitsky-Golay methods are nice to build derivative estimates of a times series. They are fast and efficient for long series of data, as they can be computed using the function filter. There are several such Savitsky-Golay variants on the file exchange. I've got one up there too, called movingslope. http://www.mathworks.com/matlabcentral/fileexchange/16997 Note that these methods are essentially fast versions of a differentiated Lagrange interpolating polynomial. In the case of movingslope, it allows you to smooth the data, in effect no longer a pure Lagrange interpolating polynomial. Thus, for linear polynomials, movingslope gives these estimates: fp_s = movingslope(y) fp_s = 0 0 0 0 0 0 0 0 0 0.5 0.5 0 0 0 0 0 0 0 0 0 But for higher order approximants, we still see ringing: fp_s = movingslope(y,5,3) fp_s = 0 0 0 0 0 0 0 0 -0.08333 0.58333 0.58333 -0.08333 0 0 0 0 0 0 0 0 Instead, try this method on a smooth function, so that ringing is not a problem. y = sin(t/3); Here the derivative function is known to be simple: (1/3)*cos(t/3). And movingslope does a decent job of estimation. fp_s = movingslope(y,5,4) fp_s = Columns 1 through 11 0.33265 0.31515 0.26186 0.18003 0.078381 -0.031895 -0.13866 -0.23016 -0.29632 -0.32986 -0.32709 Columns 12 through 21 -0.28831 -0.21779 -0.1233 -0.015229 0.094516 0.19386 0.27185 0.31993 0.33311 0.30833 Suppose we have noisy data though? The problem with noise and differentiation is that differentiation is a noise amplifier. It is a variant of ill-polsed problem. Any noise in the data will be greatly amplified when we differentiate. Here a lower order approximation will be best, as well, some smoothing will be useful. We might consider using the SLM tools as someone suggested, but they require the user to specify the knots (breaks) of the spline. This may be acceptable, but you will need to employ slightly more thought in the process. The point is, if you understand your data as well as have an understanding of the methods used to find the derivative estimates, you can find the best tool for the job. John |