Central Differences

September 26, 2008

The most common way of computing numerical derivative of a function f(x) at any point x^* is to approximate f(x) by some polynomial P_m(x) in the neighborhood of x^* . It is expected that if selected neighborhood of x^* is sufficiently small then P_m(x) approximates f(x) near x^* well and we can assume that $f’(x^*)\approx P’_m(x^*)$ .

Let’s consider this approach in details (or go directly to the table of formulas).

At first, we sample f(x) at the ( is odd) equidistant points around x^* :

$f_k = f(x_k),\: x_k = x^*+kh,\: k=-\frac{N-1}{2},\dots,\frac{N-1}{2}$ ,

where is some step.

Then we interpolate points $\{(x_k,f_k)\}$ by polynomial of (N-1)th degree: $P_{N-1}(x)=\sum_{j=0}^{N-1}{a_jx^j}$ . Its coefficients $\{a_j\}$ are found as a solution of system of linear equations:

$\left\{ P_{N-1}(x_k) = f_k\right\},\quad k=-\frac{N-1}{2},\dots,\frac{N-1}{2}$

By our assumption f’(x^*) can be approximated by the derivative of the constructed interpolating polynomial:

$f’(x^*)\approx P’_{N-1}(x^*)$

To illustrate the process let’s consider case when N=3 . Here interpolating parabola are defined by the system (we assume x^*=0 for simplicity):

$\left\{\begin{matrix} P_2(-h)&=&f_{-1}\\ P_2(0)&=&f_0\\ P_2(h)&=&f_1 \end{matrix}\right. \Rightarrow \left\{\begin{matrix} a_0-a_1h+a_2h^2&=&f_{-1}\\ a_0&=&f_0\\ a_0+a_1h+a_2h^2&=&f_1 \end{matrix}\right. $

After subtracting first equation from the last we get final result:

$f’(0)\approx P’_2(0)=a_1=\displaystyle\frac{f_{1}-f_{-1}}{2h}$

In the same way we can obtain expressions for any . Formulas for N = 3,5,7,9 are listed below:

$ \begin{tabular}{|r|c|} \hline $N$& $N$-point stencil Central Differences \\ \hline &\\ 3&$\displaystyle\frac{f_{1}-f_{-1}}{2h}$\\ &\\ 5&$\displaystyle\frac{f_{-2}-8f_{-1}+8f_1-f_2}{12h}$\\ &\\ 7&$\displaystyle\frac{-f_{-3}+9f_{-2}-45f_{-1}+45f_1-9f_2+f_3}{60h}$\\ &\\ 9&$\displaystyle\frac{3f_{-4}-32f_{-3}+168f_{-2}-672f_{-1}+672f_{1}-168f_{2}+32f_{3}-3f_{4}}{840h}$\\ &\\ \hline \end{tabular} $

These expressions are very widely used in numerical analysis and commonly refered as central(finite) differences. As it can be clearly seen they have simple anti-symmetric structure and in general difference of -th order can be written as:

${f’(x^*)\approx \frac{1}{h}\sum_{k=1}^{(N-1)/2}{c_k(f_k-f_{-k})}}$ ,

where c_k are coefficients derived by procedure described above.

It is easy to see that if f(x) is a polynomial of a degree N-1 , then central differences of order $\ge N$ give precise values for f(x) derivative at any point. This follows from the fact that central differences are result of approximating f(x) by polynomial. If f(x) is a polynomial itself then approximation is exact and differences give absolutely precise answer.

To differentiate a digital signal u_n we need to use h=1/SamplingRate and replace f_k by $u_{n+k}$ in the expressions above. In this case derivative of a signal u_n is found by

$d_n = \sum_{k=1}^{(N-1)/2}{c_k(u_{n+k}-u_{n-k})}$

Frequency response of central differences is:

$H(\omega)=2i\sum_{k=1}^{(N-1)/2}{c_k\sin(k\omega)}$

Magnitude responses for N = 3,5,7,9 are drawn below:

Red dashed line is the magnitude response of an ideal differentiator $H_d(\omega)=i\omega$ . As grows $H(\omega)$ moves closer to ideal differentiator response (by increasing tangency order with $H_d(\omega)$ at $\omega=0$ ).

In practice there is no need for ideal differentiators because usually signals contain noise at high frequencies which should be suppresed. This means that response of practical differentiator should be close to $H_d(\omega)$ on some passband interval $\omega\in[0,\omega_c]$ and should tend to zero in the upper part of Nyquist interval $\omega\in[\omega_c,\pi]$ .

From the plot we can see that central differences don’t resemble such behavior, all they care about is to get as closer as possible to the response of ideal differentiator, without supression of noisy high frequencies. As a consequence they perform well only on exact values, which contain no noise. Different technique is needed for robust derivative estimation of noisy signals. Such approaches are considered in the next sections:

(7 votes, average: 5.00)

Loading ...

Print This Page

17 Comments

Kem

Posted February 19, 2009 at 8:32 pm | Permalink

Hi Pavel, great tutorial!

How can this be done for the second derivatives f”(x) and also for the 2D case , especially how can we approximate d^2 f(x,y) / dxdy with different orders of approximation

Thank you!

Reply
Pavel Holoborodko

Posted February 20, 2009 at 5:42 pm | Permalink

Thanks, Kem.

1. $f^{\prime\prime}(x)$ is very simple to derive. We use the same interpolating polynomial and assume that $f^{\prime\prime}(x^*)\approx P^{\prime\prime}_m(x^*)$ . Second order central differences derived this way are:

$ \begin{tabular}{|r|c|} \hline $N$& Second order Central Differences \\ \hline &\\ 3&$\displaystyle\frac{f_{-1}-2f_{0}+f_1}{h^2}$\\ &\\ 5&$\displaystyle\frac{-f_{-2}+16f_{-1}-30f_0+16f_1-f_2}{12h^2}$\\ &\\ 7&$\displaystyle\frac{2f_{-3}-27f_{-2}+270f_{-1}-490f_0+270f_1-27f_2+2f_3}{180h^2}$\\ &\\ \hline \end{tabular} $

3. Third order central differences are:

$ \begin{tabular}{|r|c|c|} \hline $N$&Third order Central Differences&Error\\ \hline &&\\ 5&$\displaystyle\frac{-f_{-2}+2f_{-1}-2f_1+f_2}{2h^3}$&$O(h^2)$\\ &&\\ 7&$\displaystyle\frac{f_{-3}-8f_{-2}+13f_{-1}-13f_1+8f_2-f_3}{8h^3}$&$O(h^4)$\\ &&\\ \hline \end{tabular} $

2. Estimation of the mixed second order derivative $\textstyle{\frac{\partial^2{f}}{\partial{x}\partial{y}}}$ is a little more elaborate but still follows the same idea. In this case we should build interpolating polynomial over the 2D grid $\{f_{k,l}=f(x^*+kh,y^*+lh)\}$ and find its mixed derivative at the point assuming that it is approximately equal to $\textstyle{\frac{\partial^2{f}}{\partial{x}\partial{y}}}$ .

Here are reference formulas for mixed second order differences:

For the order of approximation:

$ \displaystyle{{\frac{\partial^2{f}}{\partial{x}\partial{y}}}\approx \frac{1}{4\,h^2}\left[f_{-1,-1}+f_{1,1}-f_{1,-1}-f_{-1,1}\right] }$

For the order of approximation:

$ \displaystyle{{\frac{\partial^2{f}}{\partial{x}\partial{y}}}\approx \frac{1}{600\,h^2} \left[\begin{matrix} -63(f_{1,-2}+f_{2,-1}+f_{-2,1}+f_{-1,2})+\\ 63(f_{-1,-2}+f_{-2,-1}+f_{1,2}+f_{2,1})+\\ 44(f_{2,-2}+f_{-2,2}-f_{-2,-2}-f_{2,2})+\\ 74(f_{-1,-1}+f_{1,1}-f_{1,-1}-f_{-1,1}) \end{matrix}\right] } $

Reply
shinji

Posted August 12, 2009 at 11:30 pm | Permalink

anyone know about the derivation of the third order central differences formula?

Reply
Chris Willy

Posted June 12, 2010 at 9:33 pm | Permalink

Hi Pavel –
Do you have a published reference for the O(h^4) approximation for the mixed 2nd partial derivatives? Thanks.

Reply
- Pavel Holoborodko
 
 Posted June 13, 2010 at 10:32 pm | Permalink
 
 Here is formula for the order of approximation:
 
 $ \displaystyle{{\frac{\partial^2{f}}{\partial{x}\partial{y}}}\approx \frac{1}{600\,h^2} \left[\begin{matrix} -63(f_{1,-2}+f_{2,-1}+f_{-2,1}+f_{-1,2})+\\ 63(f_{-1,-2}+f_{-2,-1}+f_{1,2}+f_{2,1})+\\ 44(f_{2,-2}+f_{-2,2}-f_{-2,-2}-f_{2,2})+\\ 74(f_{-1,-1}+f_{1,1}-f_{1,-1}-f_{-1,1}) \end{matrix}\right] } $
 
 Reply
 - Chris Willy
 
 Posted June 13, 2010 at 11:46 pm | Permalink
 
 Hi Pavel -
 Thanks for your quick reply. I read that you had posted that same formula earlier … in fact, it’s the earlier post that prompted me to write. I’m looking for a published reference (book or journal article) that contains the formula so that I can cite the reference in a paper I’m writing. Thanks.
 Best regards,
 Chris Willy
 
 Reply
 - Pavel Holoborodko
 
 Posted June 14, 2010 at 9:26 pm | Permalink
 
 Well, I’ve derived this formula by myself. I have no idea is it published somewhere or not.
 
 You can reference my site – it would be of huge support for me. Several published papers already include references to my site – so I think there is no problem.
 
 You can tell me more about your task, maybe I can derive more suitable filters for your conditions.
 
 Reply
Mercè Bruned

Posted July 1, 2010 at 7:15 pm | Permalink

Hello Pavel,
Congratulations for your good and well organised work. I’m not mathematic and I’m writing an algorithm to derivate a discrete function. My set of points are not equidistant, so they have x values that are no constant. I read your advanced work about derivatives for noisy functions and probably I will use it in a future.
For the moment, I will derivate with central differences method. I wonder if it’s easy to extend the method when points are not equidistant. Is it ok to divide by the x increment?
Thanks for your blog

Reply
- Pavel Holoborodko
 
 Posted July 2, 2010 at 3:41 pm | Permalink
 
 Thank you for your comment!
 
 Central difference extension for irregular spaced data is:
 
 $ \displaystyle {f’(x^*)\approx \sum_{k=1}^{(N-1)/2}{c_k\cdot\frac{f_{k}-f_{-k}}{ x_{k}-x_{-k}}\cdot 2k}} $
 
 So, yes, difference of function values should be divided by the increment with weight. Coefficients $\{c_k\}$ are from formula for uniform spacing. For instance, for $N=3,\,c_1=1/2$ and for $N=5,\,c_1=2/3,\,c_2=-1/12$ .
 
 Actually it is better to use $\{c_k\}$ from smooth noise-robust differentiators, since they are much more numerically stable even in the absence of noise. Plus they have similar approximation properties to central differences – so you won’t lose any precision and get more stability as a bonus.
 
 Coefficients are for $N=5,\,c_1=1/4,\,c_2=1/8$ and for $N=7,\,c_1=13/32,\,$ $c_2=1/8,\,c_3=-5/96$ . See referenced page for other $\{c_k\}$ .
 
 Let me know about the results please.
 
 Reply
Mercè Bruned

Posted July 2, 2010 at 8:40 pm | Permalink

Hello again Pavel,

Thank you very much for your fast and good answer.

Really, I have a noisy discrete function (points come from peak calculation of a laser stripe, acquired by a digital camera. I use the laser stripe to enhance profile of objects).

I need to smooth the function for display purpose and I need also to derivate the function in order to find inflexion points, local maxima and other features that I use to do measurements of the object that I’ve enhanced the profile.

At this moment I’m smoothing the function using Savitzky-Golay filter and I apply central differences to the smoothed function. I read your work and I realized that smooth noise-robust differentiators would be really useful to my work, but I don’t know how to use it to smooth the function for display purposes.

Reply
- Pavel Holoborodko
 
 Posted July 5, 2010 at 10:21 am | Permalink
 
 Thank you for the description of your project – it is really interesting.
 
 I think your application could benefit from using SNRD filters. I will help you to apply them for smoothing (for visual output) and for finding derivative.
 
 Just let me know what polynomial degree and number of points you are using now in Savitzky-Golay filter. Knowing that info I can build SNRD filters with better properties for your task.
 
 Reply
 - Mercè Bruned
 
 Posted July 12, 2010 at 7:45 pm | Permalink
 
 Hello Pavel, I appreciate your help very much.
 
 I’m smoothing my function with a third degree polynomial and I’m using 7 points. Once smoothed I’m applying central differences to calculate first and second derivative.
 
 I’ve tryed to calculate first derivative by using SNRD filters on the Savitzy-Golay smoothed functions, and I’m really loosing important function detail, I’m smoothing too much, I think this is not a good way to proceed.
 
 Thanks again for your interest and best regards from Spain.
 
 Reply
 - Pavel Holoborodko
 
 Posted July 13, 2010 at 12:30 pm | Permalink
 
 Congratulations on the Victory!!
 
 Thank you for the information about your test results.
 I have several comments though.
 
 I’ve checked properties of 3rd order 7 point SG and derive some special filters with better properties. I would appreciate if you would try them and give me feedback.
 
 1. For your application I recommend to use following filters for smoothing (for visualization):
 
 $ \displaystyle {\hat{f}(x_0)\approx a_0f_0+\sum_{k=1}^{(N-1)/2}{2a_k\cdot\frac{f_{k}(x_k-x_0)+f_{-k}(x_0-x_{-k})}{x_{k}-x_{-k}}}} $
 where coefficients $\{a_k\}$ are:
 $ \scriptsize{ \begin{tabular}{|l|l|} \hline &\\ N = 7&$a_0=1/2,\,a_1 = 9/32,\,a_2=0,\,a_3 = -1/32$\\ &\\ N = 9&$a_0=55/128,\,a_1= 9/32,\,a_2 = 3/64,\,a_3 = -1/32,\,a_4 = -3/256$\\ &\\ \hline \end{tabular} } $
 Both of them preserve important details better than SG and have soft and guaranteed noise suppression ( is stronger).
 
 2. For derivative estimation I recommend to use
 $ \displaystyle {f’(x^*)\approx \sum_{k=1}^{(N-1)/2}{c_k\cdot\frac{f_{k}-f_{-k}}{ x_{k}-x_{-k}}\cdot 2k}} $
 with
 $ \scriptsize{ \begin{tabular}{|l|l|} \hline &\\ N = 7&$c_1=13/32,\,c_2=1/8,\,c_3=-5/96$\\ &\\ N = 9&$c_1= 9/32,\,c_2 = 1/6,\,c_3 = -1/96,\,c_4 = -1/48$\\ &\\ \hline \end{tabular} } $
 
 Important: do not apply differentiation filters after smoothing filters. This could hurt approximation order (as you experienced with SNRD after SG smoothing).
 
 Apply derivative filter directly on noisy non-smoothed data.
 
 Reply
Jan Kadlec

Posted July 20, 2010 at 10:22 pm | Permalink

Awesome tutorial.

I’m writing a GPU raytracer and need to estimate the gradient (∂/∂x, ∂/∂y, ∂/∂z) of a distance function (for surface normals). I assume that the function is expensive to evaluate.

Now I’m using central differences – these should get me O(h²) error with seven function samples (the center is always evaluated). I was wondering if I could get the same error order with fewer samples – some kind of sample sharing between the three axes (I’ve tried the tetrahedral arrangement, but it’s even worse than forward differences because of ∂x∂y-like terms).

Reply
- Pavel Holoborodko
  
  Posted July 21, 2010 at 11:06 am | Permalink
  
  Thanks.
  
  Actually 7-points central difference has . Three point stencil is enough to achieve .
  
  To get the lowest error with the smallest number of points in 3D one should apply 1D stencil along the corresponding axis. Only points along the required axis play major role for approximation order, others can be used for other goals (like noise suppression).
  
  P.S.
  I like your site. Some of your projects made me nostalgic about old good DOS-hacking times. Trees were taller and grass was greener back then…….
  
  Reply
  - Jan Kadlec
    
    Posted July 21, 2010 at 7:02 pm | Permalink
    
    Thanks for the quick answer and for visiting my site .
    
    What I originally meant was using the three-point stencil along every axis (which makes 6 evaluations in total + one in the center for other unrelated purposes).
    
    I’ve finally settled for 4-point forward differences (1+9 evaluations in total with error O(h³)): $-2 f_{-1} -3 f_0 + 6 f_1 – f_2 \over 6h$ .
    
    Reply
    - Pavel Holoborodko
      
      Posted July 22, 2010 at 1:45 pm | Permalink
      
      Oh, sorry for misunderstanding. Glad you have found suitable solution.
      Maybe I should include one sided differences in the tutorial.
      
      Reply

Central Differences

17 Comments

Post a Comment

Search

Pages

Categories

Archives

Links

RSS Feeds

Meta