# Finite Difference Method Neumann Boundary Condition with Variable Coefficients

Disclaimer

In the process of typing up this question, I determine its solution. Since I went through the trouble of typing up the question in its entirety, I will post its answer as well. It may help out others who find themselves in the same predicament. Think of this as a sort of blog-post, if you will.

The Goal

Consider the mixed boundary value problem \begin{align} \frac{d}{dx}\left(k(x)\frac{du}{dx}\right)=f \text{ in }\Omega\\ u=P \text{ at } x=0\\ \frac{du}{dx}=T \text{ at } x=1 \end{align} where $P$ and $T$ are constants and $f$ is a source term.

I'm using finite differences and my goal is to impose the boundary condition at $x=1$ in such a way as to achieve second order accuracy. Assume that the grid has $N+1$ equispaced points (including the boundary points) given as $x_0,x_1,...,x_N$

My Approach

At the right boundary, use a 2nd order centered difference for the boundary condition $$\frac{u_{N+1}-u_{N-1}}{\Delta x}=T$$

and the 2nd derivative operator as: \begin{align} \frac{k_{N+\frac{1}{2}}\frac{u_{N+1}-u_{N}}{\Delta x} -k_{N-\frac{1}{2}}\frac{u_{N}-u_{N-1}}{\Delta x}}{\Delta x}=f_i\end{align}

We can solve for the ghost point $U_{N+1}$ in the 1st equation, substitute it into the 2nd equation and simplify.

The Problem

Using this discretization requires the evaluation of $k_{N+\frac{1}{2}}=k(x_N+\frac{\Delta x}{2})$, which is outside of the domain! In general, $k(x)$ is only defined within the domain and I can't/shouldn't use values outside of it. Thus, I don't think this is the right approach to achieve 2nd order accuracy. What else can I do in this case?

• The physically correct boundary conditions for this equation are $k \partial u/\partial x=T$ and this is how you should discretize it. – Wolfgang Bangerth Jun 17 '13 at 22:25
• @WolfgangBangerth: Very true. If the boundary condition were given as $k(x)\frac{\partial u}{\partial x}=T$, then would it be advisable to multiply $k(x_N)$ to the discretization below? – Paul Jun 17 '13 at 23:39
• You mean to the formula below "My approach"? If so, yes! – Wolfgang Bangerth Jun 18 '13 at 12:13
• @WolfgangBangerth: Actually, I was referring to the discretization technique in my answer, not the question :) – Paul Jun 18 '13 at 12:49
• Yes, there too. – Wolfgang Bangerth Jun 18 '13 at 19:23

Forget about discretizing the 2nd derivative operator since it requires evaluating at points outside the domain. Instead, just adjust the stencil for the boundary condition to use points $u_{N},u_{N-1},$ and $u_{N-2}$. By using standard method of undetermined coefficients, you can produce the coefficients necessary to achieve 2nd order accuracy. If i'm not mistaken, the discretized boundary condition in this case turn out to produce
$$\frac{du}{dx}\approx\frac{-\frac{1}{2}u_{N-2} +2u_{N-1} -\frac{3}{2}u_N}{\Delta x}=T$$