fourth-order finite difference for $(a(x)u'(x))'$

Question

Previously I asked here about constructing a symmetric matrix for doing finite difference for $(a(x)u'(x))'$ where the (diffusion) coefficient $a(x)$ is spatially varying. The answer provided there works for getting a second order accurate method. What about getting a fourth-order accurate method? If I follow the same idea and apply the following fourth-order accurate formula for first derivative

$ u'_i = \dfrac{u_{i-1} - 8u_{i-1/2} + 8u_{i+1/2} - u_{i+1}}{6\Delta x}$ (1)

in succession, then I end up getting a formula for $(a u')_i'$ which involves $u_{i-2}, u_{i-3/2}, u_{i-1}, u_{i-1/2}, u_i, u_{i+1/2}, u_{i+1}, u_{i+3/2}, u_2$. It involves the mid point values $u_{i+n/2}$ and it depends on nine neighbouring values of $u_i$, which doesn't sound right for fourth order accurate scheme. For the special case of $a(x)=1$ it also doesn't reduce to the formula

$ u'' = \dfrac{-u_{i-2} + 16u_{i-1} - 30u_i + 16u_{i+1} - u_{i+2}}{12\Delta x^2}$. (2)

So what's wrong with applying (1) in succession? What's the correct approach to get a finite difference formula for $(au')'$ with fourth order accuracy?

Edit 1: After reading this post I managed to derive a symmetric four-order centered difference scheme for $[a(x)u'(x)]'$ by applying

$ u'(x) = \dfrac{u_{i-3/2} - 27u_{i-1/2} + 27u_{i+1/2} - u_{i+3/2}}{24\Delta x} + \mathcal{O}(\Delta x^4)$

twice in succession. The resulting formula for $[a(x)u'(x)]'$ is

$ (au')'_i = \dfrac{1}{576\Delta x^2} \left( a_{i-\frac{3}{2}}u_{i-3} - 27(a_{i-\frac{3}{2}} + a_{i-\frac{1}{2}})u_{i-2} + (27a_{i-\frac{3}{2}} + 729a_{i-\frac{1}{2}}+27a_{i+\frac{1}{2}})u_{i-1} - (a_{i-\frac{3}{2}} + 729a_{i-\frac{1}{2}} + 729a_{i+\frac{1}{2}} + a_{i+\frac{3}{2}})u_i + (27a_{i-\frac{1}{2}} + 729a_{i+\frac{1}{2}}+27a_{i+\frac{3}{2}})u_{i+1} - 27(a_{i+\frac{1}{2}} + a_{i+\frac{3}{2}})u_{i+2} + a_{i+\frac{3}{2}}u_{i+3} \right)$

which involves the mid point values of $a(x)$ and has seven stencils. Is there a formula that involves even less computations / stencils(e.g. five stencils)?

Answer 1

Let's use $$ f'(x) = \frac{f_{i-3/2} - 27 f_{i-1/2} + 27 f_{i+1/2} - f_{i+3/2}}{24 \Delta x} + O(\Delta x^4) $$ as an approximation for the outer derivative. To obtain a fourth-order approximation for the $(au_x)_x$ term we need to obtain $a u_x$ at each of the $x_{i-3/2}, x_{i-1/2}, x_{i+1/2}, x_{i+3/2}$ points with fourth-order accuracy.

Assuming that we know $a$ everywhere, the problem reduces to computing $u_x$ at those points. Taking our stencil to be $(x_{i-2}, x_{i-1}, x_{i}, x_{i+1}, x_{i+2})$ we can compute fourth degree interpolating polynomial $P(x)$ for $u(x)$: $$u(x) = P(x) + O(\Delta x^5).$$ Differentiating that polynomial gives fourth-order approximation for $u_x$: $$ u_x(x) = P'(x) + O(\Delta x^4). $$

Omitting straightforward computations (which I've done using Mathics) we get the following approximations for $u_x(x_{i-3/2}), \dots, u_x(x_{i+3/2})$ $$ \begin{pmatrix} u_x(x_{i-3/2})\\ u_x(x_{i-1/2})\\ u_x(x_{i+1/2})\\ u_x(x_{i+3/2}) \end{pmatrix} = \frac{1}{24 \Delta x} \begin{pmatrix} -22 & 17 & 9 & -5 & 1\\ 1 & -27 & 27 & -1 & 0 \\ 0 & 1 & -27 & 27 & -1\\ -1 & 5 & -9 & -17 & 22 \end{pmatrix} \begin{pmatrix} u_{i-2}\\ u_{i+1}\\ u_{i}\\ u_{i+1}\\ u_{i+2} \end{pmatrix} + O(\Delta x^4). $$

Combining everything to a single formula we get $$ (a u_x)_x = \frac{1}{24 \Delta x} \times \left[\\ a_{i-3/2} \frac{u_{i+2} - 5 u_{i+1} + 9 u_i + 17 u_{i-1} - 22 u_{i-2}}{24 \Delta x} -27a_{i-1/2} \frac{-u_{i+1} + 27 u_i - 27 u_{i-1} + u_{i-2}}{24 \Delta x} +\\ +27a_{i+1/2} \frac{-u_{i+2} + 27 u_{i+1} - 27 u_{i} + u_{i-1}}{24 \Delta x} -a_{i+3/2} \frac{22u_{i+2} - 17 u_{i+1} - 9 u_i + 5 u_{i-1} - u_{i-2}}{24 \Delta x}\\ \right] + O(\Delta x^4). $$

Answer 2

A bit of work with five-point Lagrange interpolants based on a uniform grid for $a$ and $u$ shows, with \begin{align} w_{-2} &:= \tfrac{1}{144}a_{-2} - \tfrac{1}{18} a_{-1} - \tfrac{1}{12}a_0 + \tfrac{1}{18} a_1 - \tfrac{1}{144} a_2, \\ w_{-1} &:= -\tfrac{1}{18}a_{-2} + \tfrac{4}{9} a_{-1} + \tfrac{4}{3} a_0 - \tfrac{4}{9}a_1 + \tfrac{1}{18}a_2 \\ w_0 &:= -\tfrac{5}{2}a_0, \\ w_{1} &:= \tfrac{1}{18}a_{-2} - \tfrac{4}{9}a_{-1} + \tfrac{4}{3}a_0 + \tfrac{4}{9}a_1 - \tfrac{1}{18}a_2, \\ w_{2} &:= -\tfrac{1}{144}a_{-2} + \tfrac{1}{18}a_{-1} - \tfrac{1}{12}a_0 - \tfrac{1}{18}a_1 + \tfrac{1}{144}a_2 \end{align} that $(au')'(0) = (w_{-2}u_{-2}+ w_{-1}u_{-1} + w_0u_0 + w_1 u_1 + w_2 u_2)/h^2 + O(h^4)$. It typically won't be a symmetric rule.

ADDED EDIT

Since I was thinking in terms of using grid values I started with the Lagrange interpolant for $a$ \begin{align*} A(x) &= \\ & \quad a_{-2}\, \frac{x - x_{-1}}{x_{-2} - x_{-1}} \frac{x - x_{0}}{x_{-2} - x_{0}} \frac{x - x_{1}}{x_{-2} - x_{1}} \frac{x - x_{2}}{x_{-2} - x_{2}} \\ &+ a_{-1}\, \frac{x - x_{-2}}{x_{-1} - x_{-2}} \frac{x - x_{0}}{x_{-1} - x_{0}} \frac{x - x_{1}}{x_{-1} - x_{1}} \frac{x - x_{2}}{x_{-1} - x_{2}} \\ &+ a_{0}\, \frac{x - x_{-2}}{x_{0} - x_{-2}} \frac{x - x_{-1}}{x_{0} - x_{-1}} \frac{x - x_{1}}{x_{0} - x_{1}} \frac{x - x_{2}}{x_{0} - x_{2}} \\ &+ a_{1}\, \frac{x - x_{-2}}{x_{1} - x_{-2}} \frac{x - x_{-1}}{x_{1} - x_{-1}} \frac{x - x_{0}}{x_{1} - x_{0}} \frac{x - x_{2}}{x_{1} - x_{2}} \\ &+ a_{2}\, \frac{x - x_{-2}}{x_{2} - x_{-2}} \frac{x - x_{-1}}{x_{2} - x_{-1}} \frac{x - x_{0}}{x_{2} - x_{0}} \frac{x - x_{1}}{x_{2} - x_{1}} \end{align*} and the corresponding $U$ for $u$. Setting $x_{k} = kh$ for a uniform mesh simplifies the expressions and makes the template symmetric about $x=0$. Computing $(AU')'$ and evaluating it at $x=0$ gives the result which can be verified to be $O(h^4)$. The symmetric template is key.