Finite differences second derivative as successive application of the first derivative

Question

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!

Answer 1

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}$.

Answer 2

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

Answer 3

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.