Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Accurate numerical differentiation in MATLAB?

Subject: Accurate numerical differentiation in MATLAB?

From: Evan Ruzanski

Date: 2 Mar, 2010 05:30:22

Message: 1 of 6

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!

Subject: Accurate numerical differentiation in MATLAB?

From: Matt Fig

Date: 2 Mar, 2010 05:41:03

Message: 2 of 6

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

Subject: Accurate numerical differentiation in MATLAB?

From: Mark Shore

Date: 2 Mar, 2010 12:31:06

Message: 3 of 6

"Matt Fig" <spamanon@yahoo.com> wrote in message <hmi8df$bgl$1@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.

Subject: Accurate numerical differentiation in MATLAB?

From: Mark Shore

Date: 2 Mar, 2010 13:12:20

Message: 4 of 6


> 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

Subject: Accurate numerical differentiation in MATLAB?

From: John D'Errico

Date: 2 Mar, 2010 16:14:04

Message: 5 of 6

"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 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

Subject: Accurate numerical differentiation in MATLAB?

From: Evan Ruzanski

Date: 3 Mar, 2010 09:02:06

Message: 6 of 6

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 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

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us