# Trouble implementing Neumann boundary conditions because the ghost points cannot be eliminated

Neumann boundary conditions are implemented by introducing ghost points outside the domain and then using the boundary conditions to eliminate the ghost points. For example, see this question.

I have a set of three equations in variables $(u,v,w)$ which are one dimensional functions of $x$, and I cannot see how to eliminate the ghost points for the equation for the $w$ variable.

Background

The problem is something like this,

$\frac{\partial u}{\partial t} = a_u\frac{\partial^2 u}{\partial x^2} + b_u\frac{\partial u}{\partial x} + f_u(x,u,v,w) \\ \frac{\partial v}{\partial t} = a_v\frac{\partial^2 v}{\partial x^2} + b_v\frac{\partial v}{\partial x} + f_v(x,u,v,w) \\ \frac{\partial w}{\partial t} = a_u\frac{\partial u}{\partial x} +a_v\frac{\partial v}{\partial x} + f_w(x,u,v,w) \\$

The set of equations describe a two species advection-diffusion problem where the third equation couples to the other two. I am having problems applying Neumann boundary conditions to the third equation. Notice that the third equation does not contain a differential of $w$, but it does contain differentials of both $u$ and $v$.

In discretizated form (only showing the $w$ equation, where $\Delta t$ is the time step and $s_u$ and $s_v$ are the coefficients of the equations after discretization),

$w_j^{n+1} - s_u\beta(u_{j+1}^{n+1} - u_{j-1}^{n+1}) - s_v\beta(v_{j+1}^{n+1} - v_{j-1}^{n+1}) = s_u(1-\beta)(u_{j+1}^{n} - u_{j-1}^{n}) + s_v(1-\beta)(v_{j+1}^{n} - v_{j-1}^{n}) + w_j^n + \Delta tf_{j}^n$

The problem

The above is a Crank-Nicolson discretization, but for the boundary points let's use an implicit scheme (setting $\beta=1$),

$w_1^{n+1} - s_u(u_{2}^{n+1} - u_{0}^{n+1}) - s_v(v_{2}^{n+1} - v_{0}^{n+1}) = w_1^n + \Delta tf_{1}^n$

The ghost points, $u_0^{n+1}$ and $v_0^{n+1}$ can be eliminated, by using the boundary conditions for the $u$ and $v$ variables. However, applying Neumann boundary conditions to the $w$ variable gives,

$\frac{\partial w}{\partial x} = \sigma |_{j=1} \\ \frac{1}{2 \Delta x} (w_{2}^{n+1} - w_{0}^{n+1} ) = \sigma$

The problem is that the boundary equation cannot be substituted into the equation for $j=1$ because it does not contain at term for $w_{0}^{n+1}$.

Question

Do you have any suggestions on how I can apply Neumann boundary conditions to this equation?

As there are no spatial derivatives of $w$ in the PDE system, asserting Neumann boundary conditions on $w$ overdetermines the problem. Provided you have an appropriate initial condition for $w$, and thus for $w^0_j$, the discretized evolution equation for $w^{n+1}_j$ (including for $w^{n+1}_1$) you derived above fully determines $w$ for all time.
To put the above in slightly different terms, consider that, in your discretization, no boundary condition on $w$ is needed to evolve the solution forward in time. Given an initial condition $w^0_j$ defined on the boundary and interior, the evolution equations for $u$, $v$, and $w$ can all be solved without reference to values of $w$ across the boundary (i.e., ghost cells). The key is that, because there are no spatial derivatives of $w$, the ghost term $w^n_0$ does not appear in the discretized evolution equations for $u$, $v$, or $w$.
• Maybe this is one of the places where physics meets maths... the boundary condition on $w$ are dependent, meaning that $dw/dx$ is completely determined by the values of $u$ and $v$ are the boundaries. I think this gets round the system being over-constrained because I am not enforcing this boundary to take a specific value. I am thinking of applying second order accurate one-sided differences instead of second-order two-sided difference for the boundary conditions. This doesn't introduce ghost points. Sounds reasonable? May 13, 2013 at 4:58
• @boyfarrell You cannot in general add Neumann boundary conditions (e.g., constraints on $dw/dx$ at the boundary) to a fully determined system without overdetermining that system. Consider the Poisson equation $-u_{xx} = f$ on $[0,1]$ with $u(0) = u(1) = 0$. This is a fully determined system. I cannot also stipulate that $u_x(0) = 1$ without accepting the possibility that no solutions exist. This is the same situation as your problem except no boundary condition is needed to begin with on $w$.