The finite difference expressions for the first, second and higher derivatives in the first, second or higher order of accuracy can be easily derived from Taylor's expansions. But, numerically, the successive application of the first derivative, in general, is not same as application of the second derivative.

First, a case where it works. Let's say that we want to compute second derivative of function $f$ given on 3-points stencil $(i-1, i, i+1)$. The finite difference formula is: $$\left(\frac{\partial^2 f}{\partial x^2}\right)_i = \frac{1}{h^2}(f_{i-1} - 2f_i + f_{i+1})$$

This result is derived from Taylor's expansions, but it can also be interpreted in the following way. The first derivatives of the first order accuracy at the intervals $(i-1, i)$ and $(i, i+1)$ are: $$\left(\frac{\partial f}{\partial x}\right)_{i-1/2} = \frac{1}{h}(f_i - f_{i-1})$$ and $$\left(\frac{\partial f}{\partial x}\right)_{i+1/2} = \frac{1}{h}(f_{i+1} - f_{i})$$ where I use $i-1/2$ and $i+1/2$ because these derivatives are representative for the cell faces (In the first order I have actually approximated my function as piece-wise linear between the grid points $x_i$. Therefore, in every grid point the slope on the left and the right hand side of it is not the same.) The second derivative in point $i$ is now: $$\left(\frac{\partial^2 f}{\partial x^2}\right)_{i} = \frac{1}{h}(f'_{i+1/2} - f'_{i-1/2}) = \frac{1}{h^2}(f_{i+1} - f_{i} - (f_i - f_{i-1})) $$ And this is identical to the finite difference expression for the second derivative in the second order of accuracy.

I wonder if there is a similar procedure to represent the second derivative in the 4th order accuracy (on 5-points stencil) as successive application of two first order derivative of the lower accuracy (on shorter stencils)?

A naive approach would be to apply first derivatives of the second order accuracy to the stencils $(i-2, i-i, i)$ and $(i, i+1, i+2)$: $$\left(\frac{\partial u}{\partial x}\right)_{i-1} = \frac{1}{2h}(u_i - u_{i-2})$$ and $$\left(\frac{\partial u}{\partial x}\right)_{i+1} = \frac{1}{2h}(u_{i+2} - u_{i})$$ and then to find the second derivative as the first derivative of the previous two: $$\left(\frac{\partial^2 u}{\partial x^2}\right)_{i} = \frac{1}{4h^2}(u_{i+2} - 2u_{i} - u_{i-2})$$ This is obviously not correct or, at least, not the same as application of the second derivative of the 4th order straight away: $$\left(\frac{\partial^2 u}{\partial x^2}\right)_{i} = \frac{1}{12h^2}(-u_{i-2} + 16u_{i-1} + 30 u_i + 16 u_{i+1} - u_{i+2})$$

So, is there a way to reproduce the last equation as a successive combination of first derivatives of the lower accuracy order? If not, why not?

Many thanks for help! This is driving me crazy!

  • 73
  • 1
  • 4
  • Note that your finite differences for $i+1/2$ and $i-1/2$ are second-order approximations. Not at $x_i, x_{i+1}$ but at $x_{i-1/2}$ and $x_{i+1/2}$ respectively. Thus your question is not quite correct, you're combining two second-order first-derivative forumlas to obtain a second-order second-derivative formula – uranix Jul 15 '20 at 18:33

3 Answers3


Let's use the unknown coefficient approach to your problem. Assume that $$ (\Delta f)(x) = \frac{-f(x-2h)+16f(x-h)-30f(x)+16f(x+h)-f(x+2h)}{12h^2} $$ is a composition of two first order finite difference derivative formulas $$ \Delta f = \Delta_2(\Delta_1 f) $$ Each of the formulas has the form $$ \Delta_1 f = \frac{a_{-1} f(x-h) + a_0 f(x) + a_1 f(x+h)}{h}\\ \Delta_2 f = \frac{b_{-1} f(x-h) + b_0 f(x) + b_1 f(x+h)}{h}\\ $$ These formulas need to approximate first derivatives, so the following order conditions should hold: $$ a_{-1} + a_0 + a_1 = b_{-1} + b_0 + b_1 = 0\\ a_1 - a_{-1} = b_1 - b_{-1} = 1\\ $$ Composing these two formulas gives $$ (\Delta_2(\Delta_1 f))(x) = \frac{ (a_{-1} b_{-1}) f(x-2h) + (a_{-1} b_0 + a_0 b_{-1}) f(x-h) + (a_{-1} b_1 + a_0 b_0 + a_1 b_{-1}) f(x) + \dots }{h^2}\\ \frac{\dots + (a_0 b_1 + a_1 b_0) f(x+h) + (a_1 b_1) f(x+2h) }{h^2} $$ So we've arrived to a system of quadratic equations for $a_k, b_k$.

The problem is exactly the same as factorizing $$ p(x) = \frac{-x^4 + 16 x^3 - 30 x^2 + 16x - 1}{12} $$ into a product of $$ q(x) = a_{-1} x^2 + a_0 x + a_1\\ r(x) = b_{-1} x^2 + b_0 x + b_1 $$

Factoring the polynomial $p(x) = q(x) r(x)$ means that the roots of $p(x)$ are the union of the roots of $q(x)$ and the roots of $r(x)$ (including the multiplicity).

It is easy to see that $p(x)$ has root $x = 1$ with multiplicity 2 (this is a direct consequence of $\Delta$ being a second-order derivative approximation) and $q(x)$ and $r(x)$ also have the root $x = 1$ due to order conditions. $$ \frac{p(x)}{(x-1)^2} = \frac{x^2 - 14x + 1}{12}. $$ The polynomial in the right-hand side does not have real roots. This means that there is no factorization into $q(x) r(x)$ product with coefficients $a_k, b_k$ being real. There's no representation of the formula as a composition of two first-order three-points formulas.

Let's try some other forms of $\Delta_1$ and $\Delta_2$. $$ \Delta_1 f = \frac{a_{-1} f(x-h) + a_0 f(x)}{h}\\ \Delta_2 f = \frac{b_{-1} f(x-h) + b_0 f(x) + b_1 f(x+h) + b_2 f(x+2)}{h}\\ $$ Now $$ q(x) = a_{-1} x + a_0\\ r(x) = b_{-1} x^3 + b_0 x^2 + b_1 x + b_2 $$ The order conditions immediately give the solution for $q(x)$: $a_0 = 1, a_{-1} = -1$. Thus $\Delta_1$ is simply the left divided difference approximation. $$ (\Delta_1 f)(x) = \frac{f(x) - f(x-h)}{h}. $$ Finding $\Delta_2$ is straightforward: $$ r(x) = \frac{p(x)}{1 - x} = \frac{x^3 - 15 x^2 + 15 x - 1}{12}\\ (\Delta_2 f)(x) = \frac{f(x-h) - 15 f(x) + 15 f(x+h) - f(x+2h)}{12h}. $$ Verifying that $(\Delta_2 f)(x)$ actually approximates $f'(x)$ is left as an exercise.

Another solution may be obtained by taking $\Delta_1$ as right divided difference. It is pretty much the same solution with opposite signs and reflected nodes.

Another exercise: show that any finite difference formula of order $p$ can be represented as composition of $p-1$ order finite difference with $\frac{f(x) - f(x-h)}{h}$.

  • 6,759
  • 1
  • 16
  • 45
  • That's excellent! I still have to work out the details of your solution. Many thanks for the effort! The conclusion is actually very important and I am a bit surprised that I haven't find discussion of this problem in the standard textbook on numerical maths. – CFDIAC Jul 16 '20 at 19:39
  • In the context of diffusive fluxes in the heat equation or convection-diffusion equation I would interpret the $(f(x-h) - 15 f(x) + 15 f(x+h) - f(x+2h))/(12h)$ as a numerical flux at $x = x_i + h/2$ and the $(f(x) - f(x-h))/h$ as the divergence if the numerical fluxes for the $i$th cell. Everything is centered and symmetric. – Vladimir F Героям слава May 15 '22 at 19:21
  • Or let's just stay at the derivatives and say that the former approximates the first derivative at $x = x_i + h/2$ up to $-h^2 f''' / 24$. – Vladimir F Героям слава May 15 '22 at 19:32

Using these four relations:

$$f(x+\delta) = f(x) +\delta f'(x)+ \delta^2 \frac{1}{2!} f^{''}(x) + \delta^3 \frac{1}{3!} f^{(3)}(x) + \delta^4 \frac{1}{4!} f^{(4)}(x) + \delta^5 \frac{1}{5!} f^{(5)}(a) $$

$$f(x-\delta) = f(x) -\delta f'(x)+ \delta^2 \frac{1}{2!} f^{''}(x) - \delta^3 \frac{1}{3!} f^{(3)}(x) + \delta^4 \frac{1}{4!} f^{(4)}(x) - \delta^5 \frac{1}{5!} f^{(5)}(b) $$

$$f(x+2\delta) = f(x) +2\delta f'(x)+ 4\delta^2 \frac{1}{2!} f^{''}(x) + 8\delta^3 \frac{1}{3!} f^{(3)}(x) + 16 \delta^4 \frac{1}{4!} f^{(4)}(x) + 32 \delta^5 \frac{1}{5!} f^{(5)}(c) $$

$$f(x-2\delta) = f(x) -2\delta f'(x)+ 4\delta^2 \frac{1}{2!} f^{''}(x) - 8\delta^3 \frac{1}{3!} f^{(3)}(x) + 16 \delta^4 \frac{1}{4!} f^{(4)}(x) - 32 \delta^5 \frac{1}{5!} f^{(5)}(d) $$

we can show that:

$$\frac{- f(x- 2\delta) + 16f(x-\delta)-30 f(x) + 16 f(x+\delta) -f(x+2\delta) }{12\delta^2} = f^{''}(x) +O(\delta^4) $$

We note that:

$$\frac{- f(x- 2\delta) + 16f(x-\delta)-30 f(x) + 16 f(x+\delta) -f(x+2\delta) }{12\delta^2} = $$

$$ = \frac{- \frac{f(x+2\delta)-f(x+\delta)}{\delta} + \frac{ f(x-\delta) -f(x-2\delta)}{\delta} -15\frac{f(x)-f(x-\delta)}{\delta} +15\frac{f(x+\delta)-f(x)}{\delta} }{12\delta}$$

  • 6,121
  • 1
  • 14
  • 18
  • Thanks! That's correct as the way to derive the equation, but not the answer to my question. I am looking for an alternative derivation that would reproduce the same formula by taking first derivative of a first derivative. Imagine that, for whatever reason, you have available code only for computing first derivatives (in any order of accuracy) and that you have to compute the second derivative only using some combination of them. – CFDIAC Jul 14 '20 at 16:53
  • I see. You get the central second derivative (second-order approximation) from forward/backward first derivatives (first order approximation). – ir7 Jul 14 '20 at 17:02
  • @CFDIAC I added a note that might help. – ir7 Jul 14 '20 at 17:47
  • Thanks! That's smart and may be the way to go. I still don't see how to combine elementary derivatives to get factor 15 from your expression, but it looks like a promising first step. I'll try to work it out. – CFDIAC Jul 14 '20 at 17:57
  • This might help two (finite difference relations): https://en.wikipedia.org/wiki/Finite_difference – ir7 Jul 14 '20 at 18:00
  • You can also get that formula from the Richardson extrapolation of the lower order formula. – Lutz Lehmann Jul 15 '20 at 06:57

A fourth order second derivative difference scheme can be derived using successive application of a first derivative difference scheme. However, it uses a 7-point stencil, not the 5-point one.

$$ \left( \frac{\partial f}{\partial x}\right)_{x_{i+1/2}} = \frac{9}{8}\left(\frac{f(x+h)-f(x)}{h} \right) - \frac{1}{8}\left(\frac{f(x+2h)-f(x-h)}{3h} \right)\\ \left( \frac{\partial^2 f}{\partial x^2}\right)_{x_{i}} = \frac{9}{8}\left(\frac{\left( \frac{\partial f}{\partial x}\right)(x+h/2)-\left( \frac{\partial f}{\partial x}\right)(x-h/2)}{h} \right) - \frac{1}{8}\left(\frac{\left( \frac{\partial f}{\partial x}\right)(x+\frac{3}{2}h)-\left( \frac{\partial f}{\partial x}\right)(x-\frac{3}{2}h)}{3h} \right)\\ = \frac{f(x-3h) - 54 f(x-2h) + 783 f(x-h) - 1460 f(x) + 783 f(x+h) - 54 f(x+2h) + f(x+3h)}{576 h^2} $$

with the leading discretization error term according to Mathematica is

$$ -\frac{3}{320}h^4f^{(6)} $$

The five point scheme can be derived as the second derivative of the Lagrange interpolating polynomial of the five point values and the leading error term according to Mathematica is $$ -\frac{1}{90}h^4f^{(6)} $$

This difference arises, for example, in the numerical discretization of incompressible flow, where the correct Laplacian to be used for the pressure-Poisson equation is the one consistent with the discrete gradient operator used for the pressure gradient and the discrete divergence used for the continuity equation.