# Matrix representation of a finite difference with Neumann boundary conditions

Given 1D data $$[c_1,c_2,c_3,\cdots,c_N]$$ I can represent the derivative operation as a matrix product. For example, using the central difference

$$\left.\frac{\partial c}{\partial x}\right|_k \approx \frac{c_{k+1}-c_{k-1}}{2h}$$

With this I can express the derivative operator as a matrix,

$$A = \frac{1}{2h} \left[ {\begin{array}{rrrrr} * & * & 0 & 0 & \cdots \\ -1 & 0 & 1 & 0 & \\ 0 & -1 & 0 & 1 & \\ & & & & \ddots \\ \end{array} } \right] \;\;\Rightarrow\;\; \frac{\partial}{\partial x} \vec{c} \;\approx\; A\,\vec{c}$$

where the * * are set by the boundary conditions. The boundary condition that I want is that $$\mathbf{c' = 0}$$ at the boundary (i.e. Neumann boundary condition). The problem I'm having is that there are two ways to interpret what I just wrote, that lead to different matrices. I need some help in understanding the implications of choosing one or the other.

OPTION 1: Let $$x_0$$ be the left boundary value (i.e. $$x_0 = x_1 - h$$). I can interpret the boundary condition to mean that the slope between $$x_0$$ and $$x_1$$ is zero. Therefore, $$c_0 = c_1$$, and we get:

$$\left.\frac{\partial c}{\partial x}\right|_1 \approx \frac{c_2-c_0}{2h} = \frac{c_2-c_1}{2h}$$

$$\Rightarrow\;\; A = \frac{1}{2h} \left[ {\begin{array}{rrrrr} -1 & 1 & 0 & 0 & \cdots \\ -1 & 0 & 1 & 0 & \\ 0 & -1 & 0 & 1 & \\ & & & & \ddots \\ \end{array} } \right]$$

OPTION 2: I can interpret the boundary condition to mean that the slope at $$x_1$$ is zero:

$$\left.\frac{\partial c}{\partial x}\right|_1 = 0 \;\;\;\Rightarrow\;\;\; A = \frac{1}{2h} \left[ {\begin{array}{rrrrr} 0 & 0 & 0 & 0 & \cdots \\ -1 & 0 & 1 & 0 & \\ 0 & -1 & 0 & 1 & \\ & & & & \ddots \\ \end{array} } \right]$$

As far as I can see, both options are entirely defensible. I'm pretty sure I've seen both online. A similar conundrum appears when you derive the matrix for the second order finite difference. Just to add some context, the problem I'm trying to solve is that I want to evolve the advection-diffusion equation numerically:

$$\frac{\partial c}{\partial t} = - v \frac{\partial c}{\partial x} + D \frac{\partial^2 c}{\partial x^2}$$

So I have a 1D tube with some material in it and I want to see how it evolves and propagates. Neumann boundary conditions are equivalent to a closed tube where material can't escape.

Thanks in advance for the help.

• Well... there haven't been any responses. So for posterity I'll just say that I ran a toy experiment for the second derivative case, and option (1) did a much better job at conserving the amount of "material" in the "tube". So I have an empirical reason to choose (1), but I still wish I understood the choice better. – DanielC Mar 12 at 17:45