Accurate numerical differentiation in MATLAB?
in [Matlab]
|
From: Evan Ruzanski on 3 Mar 2010 04:02 Outstanding reply, John, thank you so much. "John D'Errico" <woodchips(a)rochester.rr.com> wrote in message <hmjdgc$47h$1(a)fred.mathworks.com>... > "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 |