# 9-point stencil finite difference Laplacian with variable diffusion coefficients

So I'm trying to implement a 9-point stencil discretization to the 2D difussion equation. The stencil is here.

However, most of the literature deals with a Laplacian that has a constant diffusion coefficient. i.e.

$$f_t=d\Delta f(x,y).$$

However the problem I'm dealing with has a variable diffusion coefficient, i.e.

$$f_t=\nabla\cdot (d(x,y)\nabla f(x,y)).$$

How would you implement that in a 9-point stencil? I've seen in the literature the 5-point stencil version for variable coefficient but not the 9-point one.

• What about using the product rule, $$f_t = \nabla\cdot (d(x,y) \nabla f(x,y)) = d(x,y) \Delta f(x,y) + \nabla f(x,y) \cdot \nabla d(x,y) \enspace ,$$ and then applying your finite differences? – nicoguaro Jul 31 '15 at 21:51
• @nicoguaro The usual 9-point stencil is an approximation to $\nabla^2+\frac{h^2}{12}\nabla^4$. For the Poisson equation $\nabla^2u=f$ this leads to a fourth-order-accurate scheme after using $\nabla^4u=\nabla^2f$ and so modifying $f\mapsto f+\frac{h^2}{12}\nabla^2 f$ (Section 3.4 of LeVeque's FD Methods). If you know the appropriate modification for the homogeneous heat equation here, that would answer the question. I don't see how to do it in const-coef case because $d\nabla^4u = \nabla^2u_t$ is not known. – Kirill Jul 31 '15 at 22:15
• For $u_t=\nabla^2u$ it leads to $(1+\frac{h^2}{12}L_5)u_t = d L_9 u$. If discretized with Crank-Nicolson-type method, it looks stable in the von Neumann sense. – Kirill Jul 31 '15 at 22:45

Interesting question. I too would like to know if this is mentioned anywhere explicitly. My guess is maybe this is a little too hard and not too useful to be used/taught. Perhaps someone can give a proper reference.

My calculations (see code below) suggest, if I am not wrong, that for the equation $$u_t = \mathcal{L}u, \qquad \mathcal{L}u = \nabla\cdot(D\nabla u),$$ there is such a nine-point stencil $L_9$. I am going to assume the central point of the stencil is $(0,0)$ and examine that point.

It's a quick and dirty calculation, which is why it (unusually) uses derivatives of $D$ in the stencil coefficients. Here are its coefficients (sorry about the formatting; also the matrices were rotated $90^\circ$ clockwise by accident (sorry!), but they are a pain to retype): \begin{aligned} & \frac{D}{6h^2} \begin{pmatrix} 1 & 4 & 1 \\ 4 & -20 & 4 \\ 1 & 4 & 1 \end{pmatrix} \\&+ \frac{1}{12h} \begin{pmatrix} -D_y-D_x & -4D_x & D_y-D_x \\ -4D_y & 0 & 4D_y \\ -D_y + D_x & 4D_x & D_y+D_x \end{pmatrix} \\&+ \frac{1}{24D} \begin{pmatrix} -D_x D_y & -4 D_x^2-2 D_y^2 & D_x D_y \\ -2 D_x^2-4 D_y^2 & 12 D_x^2+12 D_y^2 & -2 D_x^2-4 D_y^2 \\ D_x D_y & -4 D_x^2-2 D_y^2 & -D_x D_y \end{pmatrix} \\&+ \frac{1}{12} \begin{pmatrix} D_{{xy}} & 3 D_{{xx}}+D_{{yy}} & -D_{{xy}} \\ D_{{xx}}+3 D_{{yy}} & -8 D_{{xx}}-8 D_{{yy}} & D_{{xx}}+3 D_{{yy}} \\ -D_{{xy}} & 3 D_{{xx}}+D_{{yy}} & D_{{xy}} \end{pmatrix} \\&+ \frac{h}{24D} \begin{pmatrix} 0 & D_x D_{{xx}}+D_{{xy}} D_y & 0 \\ D_x D_{{xy}}+D_y D_{{yy}} & 0 & -D_x D_{{xy}}-D_y D_{{yy}} \\ 0 & -D_x D_{{xx}}-D_{{xy}} D_y & 0 \end{pmatrix} \\&+ \frac{h}{24} \begin{pmatrix} 0 & -D_{xyy} - D_{xxx} & 0 \\ -D_{yyy} - D_{xxy} & 0 & D_{yyy} + D_{xxy} \\ 0 & D_{xyy} + D_{xxx} & 0 \end{pmatrix} \end{aligned} It's much like the original stencil, but with some corrections. Unfortunately, they look a bit awkward to me, so I can't say much about them beyond saying that they satisfy the equations I need them to satisfy.

I calculated the coefficient by brute force: by writing down the linear equations satisfied by the stencil coefficients, such that the approximation at the point $(0,0)$ would satisfy $$L_9u = \mathcal{L}u + \beta h^2\nabla^2\mathcal{L}u + \gamma h^2 \mathcal{L}\mathcal{L}u + O(h^4).$$ The point is that this kind of method of deferred corrections only works if in the equation $\mathcal{L}u = \cdots$ the first error term of order $h^2$ looks like $h^2\mathcal{M}\mathcal{L}u$ for some suitable operator $\mathcal{M} = \beta \nabla^2 + \gamma \mathcal{L}$.

I then calculated the error term coefficients $$\beta=\frac{1}{6}, \qquad \gamma = \frac{-1}{12D(0,0)}.$$ The corresponding modified equation after space discretization is $$\left(1 + \frac{h^2}6\nabla^2 - \frac{h^2}{12D(0,0)}\mathcal{L} \right) u_t = L_9 u.$$

Here is the code that I used:

(* VarCoefNinePointStencil.nb *)

Clear[ncdiff, ncdiff2]
ncdiff[d_, u_] :=
D[d[x, y] D[u[x, y], x], x] + D[d[x, y] D[u[x, y], y], y]
ncdiff2[d_, u_] := D[d[x, y] D[u, x], x] + D[d[x, y] D[u, y], y]

Module[{L9, Lapprox, trunc, conds, sol, target, errTerm1, errTerm2,
uvars, remainder, remCond, remSol},
L9 = Table[Subscript[\[Alpha], i, j], {i, -1, 1}, {j, -1, 1}];
Lapprox = Normal@Series[
Sum[L9[[i + 2, j + 2]] u[i h, j h], {i, -1, 1}, {j, -1, 1}]
, {h, 0, 4}];
uvars = {Derivative[a_, b_][u][0, 0] :> Subscript[u, a, b],
u[0, 0] -> Subscript[u, 0, 0]};
Lapprox = Lapprox /. uvars;
errTerm1 = ncdiff[d, u];
errTerm2 = ncdiff[d, u];
errTerm1 = D[errTerm1, x, x] + D[errTerm1, y, y];
errTerm2 = ncdiff2[d, errTerm2];
target =
ncdiff[d, u] + h^2 \[Beta] errTerm1 +
h^2 \[Gamma] errTerm2 /. {x -> 0, y -> 0};
target = target /. uvars;
trunc = Lapprox - target /. {x -> 0, y -> 0};
trunc = trunc /. uvars;
conds = Flatten@
Table[Coefficient[trunc, Subscript[u, a, b]], {a, 0, 2}, {b, 0, 2}];
Print[conds // Column];
sol = Solve[Thread[conds == 0], Flatten@L9];
If[Head[sol] === Solve || sol === {}, Print["No solution."]; Abort[]];
sol = sol[[1]];
Print[sol];
Print[Lapprox - target /. sol // Expand // Collect[#, h] &];
remainder = Coefficient[Lapprox - target /. sol, h^2];
remCond =
Coefficient[remainder, #] & /@ {Subscript[u, 4, 0], Subscript[u, 3,
0]};
remSol = Solve[Thread[0 == remCond], {\[Beta], \[Gamma]}];
Print[remSol];
Print[(1/(24 h^2)) MatrixForm[
24 h^2 L9 //. Join[sol, remSol[[1]]] // Expand]];
Print[Lapprox - target //. Join[sol, remSol[[1]]] //
Collect[#, h, Factor] &];
Collect[remainder,
Union@Cases[remainder, Subscript[u, _, _], \[Infinity]], Factor]
]


I derived my variable 9-point stencil but I'm not sure if its right.

Hundsurfer in his book has a variable coefficient case for the 5-point stencil:

$$f_t = \nabla(\cdot(d(x,y)\nabla) f) = \frac{d_{i-1/2,j}(f_{i-1,j}-f_{i,j})-d_{i-1/2,j}(f_{i,j}-f_{i+1,j})}{(\Delta h^2)}+ \frac{d_{i,j-1/2}(f_{i,j-1}-f_{i,j})-d_{i,j+1/2}(f_{i,j}-f_{i,j+1})}{(\Delta h^2)}$$

It seems he just assumes four "diffusion" fluxes coming from the four edges of the cell, where the diffusion coefficient is calculated at the edge of the cell we're calculating.

At the limit of $d_{i,j} \rightarrow d$, i.e. where d is constant, we get the standard 5-point stencil that is in books, the internet etc:

$$\frac{d}{h^2} \begin{pmatrix} & 1 & \\ 1 & -4 & 1 \\ & 1 & \end{pmatrix}$$

So I naively tried a similar interpolation for the 9-point stencil, where diffusion fluxes are entering the cell from 8 directions i.e. two horizontal, two vertical, and four diagonal. I used the same weights as the constant coefficient 9-point stencil:

$$f_t = \nabla(\cdot(d(x,y)\nabla) f) = 4\times \frac{d_{i-1/2,j}(f_{i-1,j}-f_{i,j})-d_{i+1/2,j}(f_{i,j}-f_{i+1,j})}{(6\Delta h^2)}+ 4\times\frac{d_{i,j-1/2}(f_{i,j-1}-f_{i,j})-d_{i,j+1/2}(f_{i,j}-f_{i,j+1})}{(6\Delta h^2)}+\frac{d_{i-1/2,j-1/2}(f_{i-1,j-1}-f_{i,j})}{(6\Delta h^2)}+\frac{d_{i+1/2,j-1/2}(f_{i+1,j-1}-f_{i,j})}{(6\Delta h^2)}+\frac{d_{i-1/2,j+1/2}(f_{i-1,j+1}-f_{i,j})}{(6\Delta h^2)}+\frac{d_{i+1/2,j+1/2}(f_{i+1,j+1}-f_{i,j})}{(6\Delta h^2)}$$

where I use a midpoint rule to calculate the diffusion at edges e.g. $d_{i+1/2,j+1/2}=0.5(d_{i+1,j+1}+d_{i,j})$.

If you take the limit where $d(x,y) \rightarrow d$ you recover your classical 9-point stencil for constant coefficient:

$$\frac{D}{6h^2} \begin{pmatrix} 1 & 4 & 1 \\ 4 & -20 & 4 \\ 1 & 4 & 1 \end{pmatrix}$$

I'm not sure if my naive reasoning is right, but it gave me sensible simulation results, e.g. a circular profile "diffuses" symmetrically as a sphere without deformations. Can anyone comfirm if my reasoning is right or wrong? I'm still a noob.

• I tried your FD formula in Mathematica, and its truncation error is $O(h^2)$. So by itself it is not better than the 5-point stencil. Further, in the case of the original 9-point stencil the truncation error was $O(h^2)$, which required an adjustment to be made to the equation to get the desired $O(h^4)$ error, but the adjustment is missing in this answer. So I'm not sure this is right. – Kirill Aug 1 '15 at 1:54
• Also, the function $d$ is usually assumed to be known, so I don't think you need to evaluate $d_{1/2,1/2}$ using a midpoint rule, you can just evaluate it directly (less confusing, more accurate). – Kirill Aug 1 '15 at 1:56
• I think the original 9-point stencil has a truncation of $O(h^2)$ which is the same to the original 5-point stencil. But the truncation is "invariant" to rotations so it shows less of the underlying grid .e.g a circular density profile diffusing in a cartesian grid will look more square for the 5-point stencil than 9-point stencil. – mathdummy Aug 1 '15 at 1:57
• Not saying it's wrong, though: I haven't checked. Only that it is not clear from what you wrote to establish that this is correct. "Sensible simulation results" is a poor way to test convergence order of a method. It's not identical to mine: for one thing all the $D_x D_y/D$ terms are missing. – Kirill Aug 1 '15 at 1:59
• Yes. I know. You might be correct. It was a naive derivation. Your post was very helpful. – mathdummy Aug 1 '15 at 2:00