Outstanding reply, John, thank you so much.
"John D'Errico" <woodchips@rochester.rr.com> wrote in message <hmjdgc$47h$1@fred.mathworks.com>...
> "Evan Ruzanski" <ruzanski.02@engr.colostate.edu> wrote in message <hmi7pe$2or$1@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 higherorder 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, SavitskyGolay 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 SavitskyGolay 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 illpolsed 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
