The most common way of computing numerical derivative of a function
at any point
is to approximate
by some polynomial
in the neighborhood of
. It is expected that if selected neighborhood of
is sufficiently small then
approximates
near
well and we can assume that
.
Let’s consider this approach in details (or go directly to the table of formulas).
At first, we sample
at the
(
is odd) equidistant points around
:
,
where
is some step.
Then we interpolate points
by polynomial of
degree:
. Its coefficients
are found as a solution of system of linear equations:

By our assumption
can be approximated by the derivative of the constructed interpolating polynomial:

To illustrate the process let’s consider case when
. Here interpolating parabola are defined by the system (we assume
for simplicity):
![Rendered by QuickLaTeX.com \[ \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. \]](http://www.holoborodko.com/pavel/wp-content/ql-cache/quicklatex.com-48a4d87b6a1200d709c56757da50e8d8_l3.png)
After subtracting first equation from the last we get final result:

In the same way we can obtain expressions for any
. Formulas for
are listed below:
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:
,
where
are coefficients derived by procedure described above.
It is easy to see that if
is a polynomial of a degree
, then central differences of order
give precise values for
derivative at any point. This follows from the fact that central differences are result of approximating
by polynomial. If
is a polynomial itself then approximation is exact and differences give absolutely precise answer.
To differentiate a digital signal
we need to use h=1/SamplingRate and replace
by
in the expressions above. In this case derivative of a signal
is found by
![Rendered by QuickLaTeX.com \[ d_n = \sum_{k=1}^{(N-1)/2}{c_k(u_{n+k}-u_{n-k})} \]](http://www.holoborodko.com/pavel/wp-content/ql-cache/quicklatex.com-88fd4f6671027cdff86f3c88f9c9478b_l3.png)
Frequency response of central differences is:

Magnitude responses for
are drawn below:
Red dashed line is the magnitude response of an ideal differentiator
. As
grows
moves closer to ideal differentiator response (by increasing tangency order with
at
).
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
on some passband interval
and should tend to zero in the upper part of Nyquist interval
.
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:



(13 votes, average: 4.69)
47 Comments
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!
Thanks, Kem.
1. Second order central difference is simple to derive. We use the same interpolating polynomial and assume that
. Final formulas are:
3. Third order central differences are:
2. Estimation of the mixed second order derivative
is a little more elaborate but still follows the same idea. In this case we should build interpolating polynomial over the 2D grid
and find its mixed derivative at the point
assuming that it is approximately equal to
.
Here are reference formulas for mixed second order differences:
For the
order of approximation:
For the
order of approximation:
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
anyone know about the derivation of the third order central differences formula?
Hi Pavel –
Do you have a published reference for the O(h^4) approximation for the mixed 2nd partial derivatives? Thanks.
Here is formula for the
order of approximation:
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
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.
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
Thank you for your comment!
Central difference extension for irregular spaced data is:
So, yes, difference of function values should be divided by the
increment with
weight. Coefficients
are from formula for uniform spacing. For instance, for
and for
.
Actually it is better to use
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
and for
. See referenced page for other
.
Let me know about the results please.
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.
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.
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.
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):
where coefficients
are:
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
with
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.
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.
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).
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…….
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³)):
.
Oh, sorry for misunderstanding. Glad you have found suitable solution.
Maybe I should include one sided differences in the tutorial.
What would be the two dimensional 7 and 9 point formulas for (d^2/dx^2 + d^2/dy^2))f(x,y)?
What do you need them for?
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
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)
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.
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.
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!
What is the 9 point formula for the second derivative in one dimension?
You can try 9 point stencils from here Noise Robust Second Derivative.
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
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
without solving linear system.
Generation of Finite Difference Formulas for arbitrary Spaced Grids. Bengt Fornberg
I think you might find it interesting.
Thanks for the reference!
Pavel,
thanks for the nice introduction!
For mixed second derivatives it may be more efficient to use an alternative formula:
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
:
Only 4 extra values are needed here.
Cheers!
Ed, thank you for the nice formulas!!
Just checked, first one has
, second
.
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.
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
I think results of 3 and 5-points methods won’t differ very much after smoothing you’re using.
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.
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
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?
great tutorial! It helps me with a better understanding of how an approx of derivative to an arbitrary accuracy is achieved.
Your tutorial is greadfull.I want to know how explicit central differance method is apply to the wave equation
please can you help with the fomular for the order of o(h)^3
thanks
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
HI
i want fifth order central difference (3 or 5 point).can anyone help me,my email:ghanbari.ashkan@yahoo.com
Hi. I am looking for a formulation of backward finite difference in 2d. Can anyone help?
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!
Yes, it can, check this comment for formula.
do you have some kind of history or concept for this central difference method?