Central Differences

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 $N$ ( $N$ 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 $h$ 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 $N$ . 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 $N$ -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 $N$ 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:

(30 votes, average: 4.70)

80 Comments

Kem
Posted February 19, 2009 at 8:32 pm | #
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 | #
Thanks, Kem.
1. Second order central difference is simple to derive. We use the same interpolating polynomial and assume that $f^{\prime\prime}(x^*)\approx P^{\prime\prime}_m(x^*)$ . Final formulas 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 $(x^*,y^*)$ 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 $O(h^2)$ 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 $O(h^4)$ 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
- Thomas Cameron
  Posted April 12, 2011 at 10:27 am | #
  Pavel,
  I just wanted to say how much i enjoyed finding this resource as i am taking my first course in numerical differential equations. I am having some confusion based on the definitions for the central difference operator that i am given and the one you are using. Could you post some of the steps taken to arrive at the finite difference approximation of: d^2 f(x,y) / dxdy for O(h^2).
  Thank you,
  Thomas
  Reply
- Felix
  Posted October 29, 2014 at 1:59 am | #
  Hi,
  I wondered whether you followed just the same approach in deriving the $O(h^4)$ result for $\partial^2 f / \partial x \partial y$ . I once wrote a Mathematica script to compute central differences and get your results except for that case. I get:
  $\frac{\partial^2 f}{\partial x \partial y} \approx \frac{1}{144 h^2}\left[ \begin{array}{l} 8(f_{1,-2}+f_{2,-1}+f_{-2,1}+f_{-1,2})-\\ 8(f_{-1,-2}+f_{-2,-1}+f_{1,2}+f_{2,1})-\\ (f_{2,-2}+f_{-2,2}-f_{-2,-2}-f_{2,2})+\\ 64(f_{-1,-1}+f_{1,1}-f_{1,-1}-f_{-1,1}) \end{array}\right]$
  Both schemes are correct in $O(h^4)$ , you just have an additional $-\frac{289}{900}\frac{\partial^6 f}{\partial_x^3 \partial_y^3} h^4$ term. I’m just curious why the results are not the same (I can send you the Mathematica script if you like).
  Best,
  Felix
  Reply
  - Pavel Holoborodko
    Posted October 29, 2014 at 1:29 pm | #
    Felix,
    I have used polynomial of degree 5 – it is minimally required to get the $O(h^4)$ . Also I have used least-squares instead of interpolation.
    In your case – you have used polynomial of degree 6 (and not all of the combinations of $x^ny^m$ ). This doesn’t change the overall order of approximation, but error terms are (as they should) be different.
    Reply
    - Felix
      Posted October 29, 2014 at 7:12 pm | #
      Thanks! I actually used all combinations of $x^n y^m$ such that the polynomial is of degree 8. It’s probably overkill. Thanks for the advice.
      Best, Felix
      Reply
      - Pavel Holoborodko
        Posted October 29, 2014 at 7:55 pm | #
        Well, 5×5 grid has only 25 degrees of freedom. Full polynomial of 6th degree has 28 coefficients – which is already overkill:) Polynomial of 5th degree (21 coeffs) provides a good trade-off.
shinji
Posted August 12, 2009 at 11:30 pm | #
anyone know about the derivation of the third order central differences formula?
Reply
Chris Willy
Posted June 12, 2010 at 9:33 pm | #
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 | #
  Here is formula for the $O(h^4)$ 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 | #
    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 | #
      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 | #
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 | #
  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 $x$ increment with $2k$ 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 | #
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 | #
  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 | #
    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 | #
      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 ( $N=9$ 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
      - Susan
        Posted April 6, 2011 at 2:58 pm | #
        Thank you for your info on central differencing!! For a homework assignment I am given a table (similar to yours) which gives us 3 and 5 point central difference (along with 2,3, and 5 point forward and backward) formulas, but then we are asked to determine the 4 point central difference formula from the table. I’m confused because I cannot find information on this anywhere online, is it possible to do a 4 point central divided difference approximation for a first derivative?? Any help would be greatly appreciated.
Jan Kadlec
Posted July 20, 2010 at 10:22 pm | #
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 | #
  Thanks.
  Actually 7-points central difference has $O(h^6)$ . Three point stencil is enough to achieve $O(h^2)$ .
  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 | #
    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 | #
      Oh, sorry for misunderstanding. Glad you have found suitable solution.
      Maybe I should include one sided differences in the tutorial.
      Reply
Bob
Posted September 30, 2010 at 3:09 am | #
What would be the two dimensional 7 and 9 point formulas for (d^2/dx^2 + d^2/dy^2))f(x,y)?
Reply
- Pavel Holoborodko
  Posted September 30, 2010 at 10:42 am | #
  What do you need them for?
  Reply
dan
Posted October 28, 2010 at 4:10 pm | #
hey!
how can one prove that the third derivative of the central difference scheme of the function
f(x)=1/x with arguments a,b,c,d is given as
f”'(1/x)=(4ab(c-d)+4cd(a-b))/abcd
Reply
Catherine
Posted November 20, 2010 at 6:34 am | #
Thanks for putting up these solutions. Do you have any opinion on the usefulness of Dispersion-Relation-Preserving finite difference schemes* for time dependent problems? If so, do you know of any site where coefficients are listed for various N-point stencils?
* (ie C. K. W. Tam, and J. C. Webb, “Dispersion-Relation-Preserving Finite Difference Schemes for Computational Acoustics”, Journal of Computational Physics, Vol. 107, 1993, pp. 262-281)
Reply
- Pavel Holoborodko
  Posted November 20, 2010 at 11:14 am | #
  It seems that paper is not available online for free. Could you send me pdf file or scan if you have it?
  I found related papers of C.K.W. Tam:
  http://www.math.fsu.edu/~tam/Multi-size-mesh%20DRP.pdf
  http://www.math.fsu.edu/~tam/Tam_IJCFD_04.pdf
  They look interesting – I’ll read them first and then get back to your question.
  Thanks for the pointer.
  Reply
  - Pavel Holoborodko
    Posted November 29, 2010 at 2:57 pm | #
    Sorry for the late reply.
    After having read these papers I would say that DRP are almost the same to classical central differences.
    The criterion authors use to construct DRP (minimum deviation from frequency response of ideal differentiators) is actually (near-)equivalent to the exactness on polynomials condition (which is the basement of CD).
    I see no particular reason to choose DRP over CD. Actually I would suggest to use CD since they have rational coefficients, so rounding error is minimized.
    Reply
Jeramy Dickerson
Posted November 25, 2010 at 6:15 am | #
Hello,
I am trying to set up a FD equation of the form:
d/dx(a(x)*d/dx f(x))
I am looking for a five point stencil in particular.
This is basically the Poisson equation where the permitivity is not a constant.
Thanks for any help!
Reply
Daniel
Posted December 28, 2010 at 11:22 am | #
What is the 9 point formula for the second derivative in one dimension?
Reply
- Pavel Holoborodko
  Posted December 28, 2010 at 11:26 am | #
  You can try 9 point stencils from here Noise Robust Second Derivative.
  Reply
KS
Posted January 22, 2011 at 8:51 am | #
Hi Pavel,
Thanks for this post, it was very insightful. I have not seen an in depth discussion of this topic past the N = 3 case.
I am somewhat confused on the technique you used to solve for the coefficients {c_k}. You give a general expression for the central difference based derivative approximation in terms of {c_k} and say “where {c_k} are coefficients derived by the procedure described above.” But I do not see a general technique above for finding the coefficients. I can derive the coefficients for a specific h & N by setting up the linear system Ax = b and finding the 2nd row of Pinv(A). I was curious how you derived the general form in terms of the variable h.
Thanks
Reply
- Pavel Holoborodko
  Posted January 23, 2011 at 10:52 am | #
  Hi KS,
  Yes, “procedure described above” refers to the solution of Ax=b.
  Usually differentiators are considered to have anti-symmetric coefficients.
  That is what I meant under “general form” without mathematical rigor.
  There is a simple algorithm on how to derive $\{c_k\}$ without solving linear system.
  Generation of Finite Difference Formulas for arbitrary Spaced Grids. Bengt Fornberg
  I think you might find it interesting.
  Reply
  - KS
    Posted January 23, 2011 at 2:59 pm | #
    Thanks for the reference!
    Reply
Ed
Posted February 3, 2011 at 11:06 pm | #
Pavel,
thanks for the nice introduction!
For mixed second derivatives it may be more efficient to use an alternative formula:
$\frac{\partial^2{f}}{\partial{x}\partial{y}}}\approx \frac{ - f_{+1,-1} - f_{-1,+1} + f_{1,0} + f_{-1,0} + f_{0,1} + f_{0,-1} - 2 f_{0,0}}{2\,h^2}$
Its advantage is that evaluating the mixed derivative along with the “diagonal” derivatives (d^2 f / dx^2 and d^2 f / dy^2) requires only 2 extra function values instead of 4.
Same idea can be used to derive formula accurate to $h^3$ :
$\frac{\partial^2{f}}{\partial{x}\partial{y}}}\approx \frac{1}{24\,h^2} \left[\begin{matrix} 16 (f_{0, -1} + f_{0,1} + f_{-1,0} + f_{1,0})\\ -(f_{2,0} + f_{-2,0} + f_{0,2} + f_{0,-2})\\ -16 (f_{-1, 1} + f_{1, -1})\\ + f_{-2, 2} + f_{2, -2}\\ -30 f_{0, 0} \end{matrix}\right]$
Only 4 extra values are needed here.
Cheers!
Reply
- Pavel Holoborodko
  Posted February 4, 2011 at 2:31 pm | #
  Ed, thank you for the nice formulas!!
  Just checked, first one has $O(h^2)$ , second $O(h^4)$ .
  Interesting how they behave in frequency domain – I’ll take a look on the transfer functions to check isotropy (important for 2D formulas).
  Maybe I should start collecting all the formulas in one big compendium.
  P.S.
  You can insert math formulas into comments by using LaTeX syntax:
  $ ... $ – for inline formulas
  \[ ... \] – for displayed.
  Reply
- Nana
  Posted July 7, 2014 at 2:33 pm | #
  Hi,
  I just had a comment.
  Isn’t the cross derivative Ed posted asymmetric (with respect to the diagonal terms)? It would seem to me that it might introduce some anisotropy into the calculation.
  cheers,
  Reply
  - Pavel Holoborodko
    Posted July 8, 2014 at 7:34 pm | #
    Mixed derivative should be symmetric along diagonal (depending on order of differentiation) – and Ed’s formulas are.
    If you change the order of differentiation (dxdy – > dydx), then you have to “rotate” the formula so that isotropicity (with respect to “x” and “y”) is preserved.
    Reply
Jono
Posted March 14, 2011 at 8:00 pm | #
Pavel
Thanks for this interesting intro. I am not mathematical in any way and am trying to argue the case over the advantages of a 5-point derivative method over say a 3-point method.
My application is in using the first and second derivative of angular displacement data to yield velocity and acceleration. I have smoothed the displacement data, collected at 100Hz, with a forward and reverse 5th order butterworth filter (1/50 ‘low’) prior to calculating velocity and acceleration. However I am confused as to the advantages and disadvantages of the 3 and 5-point methods.
Thanks
Reply
- Pavel Holoborodko
  Posted March 15, 2011 at 11:05 am | #
  I think results of 3 and 5-points methods won’t differ very much after smoothing you’re using.
  Reply
Lee-Ping
Posted April 28, 2011 at 1:08 am | #
Hi there,
Your website was very helpful! Thank you. My scientific application is a least-squares minimization in a large parameter space where some of the parameteric derivatives are known analytically but others are unknown. I used finite difference derivatives to estimate the gradient and diagonal elements of the Hessian, and I fill in the rest of the Hessian elements using BFGS.
For my application, I checked the three-point difference result against the seven-point difference result and got agreement to within five decimal places. This means that the three-point difference is usable in the actual application. 🙂
Reply
Tal
Posted May 13, 2011 at 8:03 am | #
Hello all, i’ve wondered if anyone can help me with the following:
I have a poisson equation: (d^2/dx^2+d^2/dy^2)=-1
Let’s say that around point U(i,j) i have 8 points as follows:
U(i-1,j-1), U(i+1,j+1), U(i,j+1), U(i,j-1), U(i-1,j), U(i+1,j), U(i-1,j+1), U(i+1,j-1)
How would the solution for U(i,j) be if DeltaX and DeltaY (finite differences) are constant but different (DeltaX not equal DeltaY)?
Thanks for the help
Reply
razye
Posted May 26, 2011 at 1:36 pm | #
hello,what is the discretization of Hessian term ([(∂^2 f)/〖∂x〗^2 ][(∂^2 f)/〖∂y〗^2 ]-〖((∂^2 f)/∂x∂y)〗^2)with 9-point stencil?
Reply
BiLinear Approx
Posted June 17, 2011 at 12:30 am | #
great tutorial! It helps me with a better understanding of how an approx of derivative to an arbitrary accuracy is achieved.
Reply
Tharanga
Posted August 16, 2011 at 3:29 pm | #
Your tutorial is greadfull.I want to know how explicit central differance method is apply to the wave equation
Reply
ochi
Posted August 20, 2011 at 8:05 pm | #
please can you help with the fomular for the order of o(h)^3
thanks
Reply
Lutfi Mulyadi
Posted September 19, 2011 at 11:01 am | #
Dear Mr Pavel
Let me introduce myself. my name Lutfi mulyadi I come from Indonesia. I was very interested when I saw your blog. now I’m working on my BSc thesis in the field of nuclear reactor physics in Bandung Institut Of Technology, Indonesia. My thesis title is Fast Nuclear Reactor Burnup Analysis Using Fourth Order Runge Kutta Methods. in this thesis I solve the boundary value problem with the central finite difference and then the initial value problem with 4 Order Runge Kutta Methods. I have managed to create a program for the Runge Kutta method but I have not succeeded for the central difference to the neutron diffusion equation . if you can analyze the mistakes that I do? I really hope you can help me. a have sent to you my computer code in C programming language for solving neutron diffusion equations with central difference and I attach also my thesis to your email that include the numerical methods that I use in chapter 3 Simulation Methods (Numerical Methods) .
Thank you so much for your help
Reply
ashkan
Posted October 16, 2011 at 11:50 pm | #
HI
i want fifth order central difference (3 or 5 point).can anyone help me,my email:ghanbari.ashkan@yahoo.com
Reply
Mike
Posted November 18, 2011 at 1:21 pm | #
Hi. I am looking for a formulation of backward finite difference in 2d. Can anyone help?
Reply
Mitzey
Posted January 10, 2012 at 2:05 pm | #
Hi ! can the central difference formula be used even when the interval between x values is not constant? if not, what formula can be applied? thanks!
Reply
- Pavel Holoborodko
  Posted January 10, 2012 at 2:54 pm | #
  Yes, it can, check this comment for formula.
  Reply
franzes
Posted February 16, 2012 at 12:24 am | #
do you have some kind of history or concept for this central difference method?
Reply
Amin
Posted March 30, 2012 at 10:56 pm | #
Very nice work Pavel!
I was just wondering whether you have formulated second order accurate central difference schemes for mixed third derivatives such as $\frac{\partial f^3}{\partial x^2 \partial y}$ and $\frac{\partial f^3}{\partial x \partial y^2}$ and possibly forth derivative $\frac{\partial f^4}{\partial x^2 \partial y^2}$
I would greatly appreciate your help
Reply
- Amin
  Posted March 31, 2012 at 1:07 am | #
  By applying the second order central difference schemes in turn, e.g. $\frac{\partial}{\partial x}(\frac{\partial ^2f}{\partial y^2})$ , I calculated the following formulae:
  $\frac{\partial^3f}{\partial x \partial y^2}=\frac{1}{2\Delta x (\Delta y)^2}(f_{i+1,j+1} - f_{i-1,j+1}-2f_{i+1,j} +2f_{i-1,j}+f_{i+1,j-1} - f_{i-1,j-1})$
  $\frac{\partial^3f}{\partial x^2 \partial y}=\frac{1}{2(\Delta x)^2 \Delta y}(f_{i+1,j+1} + f_{i-1,j+1}-2f_{i,j+1} +2f_{i,j-1}-f_{i+1,j-1} - f_{i-1,j-1})$
  $\frac{\partial^4f}{\partial x^2 \partial y^2}=$
  $\frac{1}{(\Delta x)^2 (\Delta y)^2}(f_{i+1,j+1} -2f_{i,j+1}+f_{i-1,j+1} -2f_{i+1,j}+4f_{i,j}-2f_{i-1,j} +f_{i+1,j-1}-2f_{i,j-1}+ f_{i-1,j-1})$
  where i and j refer to the nodes in the x and y directions respectively and $\Delta x$ and $\Delta y$ are the grid spacings .
  Are these correct? Are these second order accurate?
  Reply
  - Pavel Holoborodko
    Posted March 31, 2012 at 1:57 pm | #
    Hi Amin, at first glance your formulas seems to be all right, however I haven’t check their asymptotics yet. Where do you use such high order derivatives?
    Reply
    - Amin
      Posted April 2, 2012 at 1:08 am | #
      Thanks for the quick reply.
      I’m looking at implementing an ADI method to solve the unsteady convection-diffusion equation in 2D and came across a paper on the topic by Karaa and Zhang*. After expanding the formulae used I came across the following difference operators $\delta _x ^2\delta _y ^2$ , $\delta _x ^2\delta _y$ and $\delta _x\delta _y ^2$ where $\delta _x ^2$ and $\delta _y ^2$ are described as
      $\delta _x ^2 f=\frac{f_{i+1,j}-2f_{i,j}+f_{i-1,j}}{\Delta x^2}$
      $\delta _y ^2 f=\frac{f_{i,j+1}-2f_{i,j}+f_{i,j-1}}{\Delta y^2}$
      I assumed these terms represent central difference approximations to $\frac{\partial^4}{\partial x^2 \partial y^2}$ , $\frac{\partial^3}{\partial x^2 \partial y}$ and $\frac{\partial^3}{\partial x \partial y^2}$ respectively, hence I became interested in these high order derivatives.
      *High order ADI method for solving unsteady convection-diffusion problems, Journal of Computational Physics (2004)
      Reply
  - Pete Rodriguez
    Posted January 6, 2014 at 1:43 pm | #
    This is a great blog everyone. An exceptional reference book for finite difference formulas in two dimensions can be found in “modern methods of engineering computation” by Robert L, Ketter and Sgerwood P. Prawel, Jr.. This is a 1969 book but it is a jewel. I have derived the 4 partial difference “molecule” for (∂^4 w)/(∂x^3 ∂y) but I am getting some numerical overflow problems in my computer model. Does anyone have this derivation? I used the Taylor series for i+1 and i-1 and took the partial derivative wrt x and y to obtain my answer. I think it is ok but I can’t find anywhere that shows that derivation.
    Thanks, Pete
    Reply
    - Pavel Holoborodko
      Posted January 6, 2014 at 1:55 pm | #
      Pete, thank you for your feedback!
      Common reason for overflow (in differences) is cancellation error. You might try to increase h – sampling step.
      If this won’t help – post the coefficients, it will be easy to check asymptotic properties.
      Reply
  - Almir Tricic
    Posted September 26, 2020 at 5:28 am | #
    Hey!
    Can you inform me of the process you used to derive this form of the 3rd der? ie the stencils i guess. And is this 2nd order accurate?
    Thank you
    Reply
Arman
Posted December 23, 2012 at 5:11 am | #
hi,i want to solve this differential equation with finite differences and cetered on bu i do not have any idea how to pleas help me
U”(x)=U^2(x) U(0)=0 and U(1)=1 (0,1)
u”=(u(i+1)-2*u(i)+u(i-1))/(h^2)
u’=(u(i+1)-u(i-1)/(2*h))
but i dont know how to solve and where i should start from please help me
Reply
George
Posted May 12, 2013 at 10:10 pm | #
Hi! I just wish to note that for high order precision formula derivation you can use some polynomial expansions, which requires only integrability of the function.
e.g. : $f=\lim\sum\limits_{i=0}^n{|P_n>}$ , where, for instance, ${P_n}$ is normalized Legendre’s polynomials, scaled over interval in question. If f is in $L^2_{(a,b),\mu}$ it should be noted that $f'=\lim\sum\limits_{i=0}^n{|P_n'>}$ . The inner products may be obtained from Gauss Legendre quadratures.
Reply
- Pavel Holoborodko
  Posted December 4, 2013 at 4:49 pm | #
  Hi George!
  Thank you for your comment.
  There are many possible types of discretization to compute the derivatives.
  We consider uniform sampling above.
  However the best one are non-uniform (e.g. derived from Chebyshev or Legendre polynomials).
  I love the Trefethen’s book on the subject: Spectral Methods in MATLAB.
  Another useful read: Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations
  Reply
Simon
Posted December 1, 2013 at 6:29 pm | #
I have a smooth function of n variables. At a certain point I have the value and all the first order order derivatives (df/dxi) and all the non-mixed second order derivatives (d^2f/dxi^2) for i = 1,2,…n. My function is expensive to calculate. Are there any methods to estimate which mixed second order derivatives are likely to be the most material so I can focus on them? thank you!
Reply
Pete Rodriguez
Posted January 7, 2014 at 2:06 pm | #
Pavel: Thank you so much for your reply. Here are the coefficients I came up with. They seem logical to me but I am not sure how to verify.
for (∂^4 w)/(∂x ∂y^3) I get
[-w(j-2,k+1)+2w(j-1,k+1)-2w(j+1.k+1)+w(j+2,k+1)+w(j-2,k-1)-2w(j-1,k-1)+2w(j+1,k-1)-w(j+2,k-1)]/4∆x∆y^3
where the j are along the x axis and the k are along the y axis. the same process was repeated for (∂^4 w)/(∂x^3 ∂y) and the form was the same with obvious diffeences in the location of the grid points. Is this correct?
The reason I am using this higher order difference is for the numerical solution of a composite anisotropic plate where the bending and twisting coupling terms of the stiffness matrix are not equal to zero. Those stiffness coefficients are mutliplied by these two derivatives and go into the final set of equations [A]{w}={q}. Where the [A] matrix has these coefficients embedded in the solution. I have solved this same problem with the same number of grid points but assuming isotropic properties and I get within 1% of the analytical solution. But, of course, these two derivatives never come into play since for isotropic materials there is no bending -twisting coupling and those terms are zero.
Thanks so much.
Pete
Reply
- Pavel Holoborodko
  Posted January 7, 2014 at 2:52 pm | #
  Hi Pete,
  This formula looks a little strange to me. It is rectangular along the x (I rather expect it to be rectangular along y).
  So I have derived my own formula for the case. Here it is:
  (1) $\begin{multline*} \frac{\partial^4f}{\partial x \partial y^3} = \frac{1}{20\,{h}^{4}}\,({2\,f_{{-2,-2}}-4\,f_{{-2,-1}}+4\,f_{{-2,1}}-2\,f_{{-2,2} }+f_{{-1,-2}}-2\,f_{{-1,-1}}+2\,f_{{-1,1}} \\ -f_{{-1,2}}-f_{{1,-2}}+2\,f_ {{1,-1}}-2\,f_{{1,1}}+f_{{1,2}}-2\,f_{{2,-2}}+4\,f_{{2,-1}}-4\,f_{{2,1 }}+2\,f_{{2,2}}}) \end{multline*}$
  It assumes h=∆x=∆y. Also I have omitted (j,k) to save the space.
  This formula is of $O(h^2)$ and looks more logical. Could you try it for your problem and let me know how it works?
  Reply
Pete Rodriguez
Posted January 11, 2014 at 10:40 pm | #
Thank you Pavel. I will try your solution and see what happens. It seems to be more than O(h^2) though. Which may be exactly what I need. Also my apologies when I submitted my derivation previously. What I submitted was (∂^4 w)/(∂x^3 ∂y) which is, of course, rectangular along the x. I tried my derivation a couple of different ways and I keep coming up with the same answer so you are probably correct in that I need to change my grid size. The simplest approach I tried was simply to take the derivative of (∂^3 w)/(∂x^3 ) wrt y or ∂/∂y ((∂^3 w)/(∂x^3 )). I know that my derivation for (∂^3 w)/(∂x^3 ) is correct so simply taking the derivative of it wrt y should give me the right solution and it seems to work. I hadn’t thought about encountering truncation errors so soon though. Thanks so much for your insight.
Best Regards, Pete
Reply
Vladimir
Posted January 12, 2014 at 11:12 pm | #
For second order vibrational perturbationtheory (PT2) as discribed in http://arxiv.org/pdf/1306.4212.pdf are required the calculation of third and forth order numerical-derivatives. Namely dxdxdy, dxdydz, dxdxdydy.
Can anyone propose a formulas for calculation this derivatives.
Reply
Pete Rodriguez
Posted January 19, 2014 at 10:35 pm | #
Pavel:
I just wanted to thaank you for this forum and for your responses to so many other researchers. Tthere are many questions that just can’t be answered in books or papers. Many times difficulties with the mechancs of the detail derivations can put bring your work to a halt. It is always of benefit to be able to talk to colleagues with like interests. I used your derivation and I was still getting the overflow problem. I was at a complete loss and then I decided to normalize my equation solver (Gaus-Jordan Reduction) with repespect to my largest stiffness coefficient. The problem works beautifully now. Once again, thank you professor Holoborodko!
Pete
Reply
- Pavel Holoborodko
  Posted January 23, 2014 at 3:11 pm | #
  Dear Pete,
  Thank you very much for your supportive words and feedback!
  And I am glad you managed to overcome numerical difficulties in your solver.
  I don’t know, but maybe you can use automatic differentiation? It is usual choice when stencils lead to overflow.
  (btw I have no professor degree ;)).
  Reply
sami ullah khan
Posted April 2, 2014 at 8:57 pm | #
Hi….
i want to derive formula of 5th order derivative for finite difference scheme..there are some issues …help me in finding this..any body knowing this please share with me…thanks for supporting me..
Best Regards,
Sami
Reply
Amr
Posted November 4, 2014 at 11:09 pm | #
Dear Pavel Holoborodko
Can you find a representations for higher derivatives for the function f
f_xxxx, f_yyyy, f_xxxy, f_xyyy and f_xxyy
Reply
- Pavel Holoborodko
  Posted November 5, 2014 at 10:00 am | #
  Dear Amr,
  Many people asked the same question in the comments. I think it is time to do this and collect everything in one reference.
  I plan to derive & publish the formulas over next few days.
  Reply
Amr
Posted November 5, 2014 at 7:03 pm | #
Thank you Pavel Holoborodko
I try to get an expansions for the derivatives $f_{xx} , f_{yy}_.... f_{xxxx}$ and $f_{yyyy}$ by using Taylor series for a function in two variables $f(x,y)$ .
I solve a matrix equation $A_{8*14} X_{14*1}=B_{8*1}$ . But the problem is $A$ is a singular matrix. I don’t know what I do.
Reply
Mike S
Posted December 14, 2014 at 10:34 am | #
Hi Pavel,
I’ve been working with the Transonic Small Disturbance Equation using a non-uniform grid, and am having issues deriving a second order central difference equation to use. Any guidance on where to go would be fantastic.
I understand that a first order is
$\dfrac{u_i_+_1-u_i_-_1}{\Delta x_i - \Delta x_i_-_1}$
so would that make a second order simply
$\dfrac{u_i_+_1-2u_i+u_i_-_1}{(\Delta x_i - \Delta x_i_-_1)^2}$
Thank you very much!
Reply
L
Posted August 23, 2015 at 5:17 pm | #
Thank you so much for all formulas!
How can I find the formulae with lowest accuracy for derivative of a function of three variables, d^3f(x,y,z)/dxdydz
Thank you in advance.
Reply
Jorge
Posted May 22, 2016 at 4:29 pm | #
Hello Pavel,
Currently I am working on the development of schemes to solve parabolic equations. I have one scheme where I define my unknown at the nodes of my grid and its gradient at integration points (using Gauss-Legendre quadrature) of a given stencil, e.g. for a 2 node stencil I compute the gradient using the values at the nodes, which is in fact the central difference for one integration point (since it resides at the middle). The problem is that in the case where I have 3 (equidistant) nodes, I must use 3 integration points, this rises me some troubles since the integration points are not equidistant therefore I can not use exactly the Central Differences to compute the gradient at those points. I need some kind of “skew-Central Differences” I am trying to come out with a way of using a convex combination at the nodes to overcome this issue but I still have not figure it out how. While reading your post it is evident your expertise in Discrete computations and besides I have enjoyed going trough all the comments and replies, it is inspiring finding this sites like yours. Anyways before diverging more, I would like to know if you have encounter an alike problem before or if you could direct me towards an idea, site or book to find a way over my problem. By the way in my case N is the number of nodes in the stencil.
Thanks Pavel
Jorge
Reply
Michael Schwegel
Posted June 20, 2017 at 9:32 pm | #
Hello Pavel!
Thank you for your nice explanations.
I wonder how you solve for the frequency response of the finite difference scheme:
$H(j\omega) := \sum_{n \in N} h[n] e^{-j\omega n}$
writing out the expression for the derivatives in $h[n]$ i get a double sum, one over n and an inner one over k from the polynomial. I see you use the relation
$2i \sin(\omega k) = e^{j\omega k} - e^{-j\omega k}$
but I can’t figure out how you deal with the terms $u[n+k]$ and how the k index ends up in the exponent. Please explain to me 1 or 2 lines how you end up at that result!
Best regards
Mike
Reply
- Pavel Holoborodko
  Posted June 21, 2017 at 12:04 pm | #
  To derive the frequency response, just insert $e^{i\omega\,t}$ into your formula and see what multiplier it will have after simplifications.
  This page has example for smoothing filter: http://www.holoborodko.com/pavel/2011/01/16/richard-hamming-digital-filters/
  Reply
Alejandro
Posted August 3, 2018 at 7:19 am | #
Pavel, Thank you so much for sharing this information about Finite Differences. I have one question for you regarding high order finite differences:
Based on the following linkhttps://my.siam.org/books/ot98/sample/OT98Chapter1.pdf I would like to obtain a fourth order approximation of the second order derivative using a forward and then a backward first order finite difference scheme.
If I define:
$U_0^*= (-3/2U_0+2U_1-1/2U_2)/h$
$U_{-1}^*= (-3/2U_{-1}+2U_0-1/2U_1)/h$
$U_{-2}^*= (-3/2U_{-2}+2U_{-1}-1/2U_0)/h$
Then I can calculate the second derivative applying a first order backward derivative:
$d^{2}u/dx^2= (3/2U_0^*-2U_{-1}^*+1/2U_{-2}^*)/h$
Using the previous results this leads to:
$d^{2}u/dx^2= (3/2(-3/2U_0+2U_1-1/2U_2)-2(-3/2U_{-1}+2U_0-1/2U_1)+1/2(-3/2U_{-2}+2U_{-1}-1/2U_0))/h^2$
$d^{2}u/dx^2= (-3/4U_{-2}+4U_{-1}-26/4U_0+4U_{1}-3/4U_{2})/h^2$
I don’t understand why the results is not the correct approximation:
$d^{2}u/dx^2= (-1/12U_{-2}+4/3U_{-1}-5/2U_0+4/3U_{1}-1/12U_{2})/h^2$
Maybe you can help me,
Thank you
Reply