# Computing second derivatives with Neumann boundary condition

I am implementing a finite difference method for a PDE with a Neumann boundary condition. I will simplify my question to a single dimension.

Suppose I have a PDE

$$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \text{ (other terms)}$$

on the interval $$[a, b]$$. I impose on the PDE the Neumann boundary conditions

$$\frac{\partial u}{\partial x} = 0$$

at $$x=a$$ and $$b$$.

In my finite difference scheme, I need to approximate $$\frac{\partial^2 u}{\partial x^2}$$ at the points $$x = a$$ and $$x = b$$. How should I do this?

One solution I found online is to introduce the virtual points $$a' < a$$ and $$b'>b$$ to use the central schemes at $$a$$ and $$b$$. The resulting approximations tend to vary strongly depending on how these virtual points are calculated. Are there better ways to construct these virtual points?

Let's do it for the left boundary point $$x=a$$, and for simplicity assume $$a=0$$. To implement $$u_x$$=0 at $$x$$=0, assume the function $$u(x)$$ is even in the vicinity of $$x$$=0, and let the grid points be $$x_0=0$$, $$x_1=h$$, etc. Then, using the symmetry of the function, $$u(-x_1)=u(x_1)$$, the second derivative evaluated with the standard second-order accurate central difference at $$x=0$$ becomes
$$u''(0) = \frac{1}{h^2} (u(x_1) + u(-x_1) - 2 u(x_0)) = \frac{2}{h^2} (u_1 - u_0)$$
It is easy to see that this expression is just the Taylor expansion at $$x$$=0, with $$u'(0)$$=0. Another way to look at it, we are creating a virtual (ghost) grid point $$x_{-1} = - x_1$$. For a higher order scheme, with this approach one would need to introduce two or more ghost points. But there is no arbitrariness here in the choice of the ghost points (they are just mirror images of the actual grid points), and as long as the solution is resolved on the grid the result should converge to the exact answer according to the formal accuracy order of the scheme.
• To estimate $u_{xx}$ at $x=0$ you can use the non-centered rule $u_{xx}(0) = (2u(0) - 5u(h) + 4 u(2h) - u(3h))/h^2 + O(h^2)$. Sep 12 at 16:44
• That's right - but that's not enough since we also need to enforce $u'(0)=0$. Sep 12 at 17:26
• Boundary conditions are taken care of in the formulation of the problem. The four-point rule is an a posteriori estimator for $u_{xx}(x=0)$. I might be answering the wrong problem. Sep 12 at 18:20
• How would you enforce $u^{'}(0)=h$ using your method? Sep 13 at 4:17
• @Nachiket If you mean $u'(0)=C$, where C is a given constant, then I would use a ghost grid point $x_{-1}$ and set $u_{-1} = u_1 - 2 h C$, where $h$ is the grid spacing $h=x_1-x_0$. Sep 13 at 12:12