<s>[INST] fourth-order finite difference for $(a(x)u'(x))'$ 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)?
 [/INST] 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.
</s>