The effect of the boundary condition on the convergence of the fdm scheme

I know the boundary condition is usually a tricky question. However, I am testing a finite-difference scheme for the equation of the form $$u_t=a(x)u_{xx}$$ that I know the analytical solution of. So if I specify the Dirichlet boundary condition, I can see the perfect convergence with a desired rate, but if I specify the value of the derivative I see the deterioration of accuracy of the numerical solution at the boundary(still have a good convergence in the middle of the grid). Both are exact from the analytical solution but looks like the derivative is not a good choice. What I think is that there should be only one solution that satisfies the initial condition and such a value of the "slope", or I am wrong? I might have a bug in my code but could not spot it so I wonder if imposing the partial derivative on the boundary instead of the value of the function could destroy uniqueness and this is why I am converging to a slightly different solution?

• What do you mean "both are exact from the analytical solution"? If both methods give correct results, what's your problem?
– Dan
Commented Jan 22, 2013 at 19:54
• It sounds like you might have a problem with the order of approximation of your neumann (derivative) boundary condition. If you're using a centered difference scheme for $u_{xx}$, then you're tempted to expect your method to be $O(\Delta x^2)$. But if your boundary condition is implemented with a forward or backward difference (i.e. $O(\Delta x)$ accuracy), then the whole scheme will only be $O(\Delta x)$ in space.
– Paul
Commented Jan 22, 2013 at 20:16
• Are you using neumann boundary conditions on all boundaries? If so, then you have an ill-posed problem with infinitely many solutions.
– Paul
Commented Jan 22, 2013 at 20:18
• I agree with @Paul that you might not discretize your Neumann boundary conditions correctly. You could use the ghost points outside the domain to obtain second-order convergence. Commented Jan 22, 2013 at 22:09
• @Dan: that means I can compute the value of $u_x$ and $u$ on the boundaries exactly and set it to those values, but I do find numerically the value inside of the domain. What I see that if I plot the error it has significantly larger value close to the boundaries if I set the derivative instead of the function. Commented Jan 23, 2013 at 1:59

As Hui pointed out above, to apply Neumann boundary conditions correctly you should utilize ghost points and extend your $n^{th}$ order discretisation stencil to your domain boundaries. Utilising forward/backward difference or extrapolation at the boundary will degrade your solution.

Assuming a $2^{nd}$ order central difference scheme the Neuman boundary condition, $u_x(t,x=0) = p_0(t)$, becomes (omitting $(t)$'s for convenience):

$$u_x\rvert_0 = p_0 \approx \frac{u_1 - u_{-1}}{2h}$$

$$\Rightarrow u_{-1} = u_1 - 2hp_0$$

Now, substitute for $u_{-1}$ into your discretised equation at the boundary, to obtain a Dirichelt BC $u_0$:

$$u_0(t) = a(0)\frac{u_1 - 2u_0 + u_{-1}}{h^2}$$

$$\Rightarrow u_0 = a(0)\frac{2u_1 - 2hp_0}{h^2 + 2a(0)}$$

Applying the Neumann boundary condition in this way is consistent and is strictly equivalent to applying the analytic Dirichelt bc in the limit $h\rightarrow 0$.

As for which is more accurate, for finite $h$, if you apply the Dirchlet condition at the boundary you impose exact agreement between the analytic and numerical solution at the boundary, but your estimation of the solutions derivative at this point will contain a $2^{nd}$ order error. Applying the Neumann boundary condition simply does the reverse, i.e. yields exact agreement between the analytic and numerical values of the derivative of your solution (estimated using your finite difference stencil) at the boundary, while estimations of $u$ at this point will contain a $2^{nd}$ order error. In the limit $h \rightarrow 0$ both methods should converge to the same unique solution.

If you are not utilising ghost points then you are not applying the Neumann conditions in a consistent way and your result will indeed simply be less accurate then when you solve using the Dirichelt boundary conditions.

• I agree with you that I must have a bug in the code as both approximations should be consistent. However, I see no problem approximating the derivative from one side by setting $u_x(0)=-3/(2h)u_0+2/hu_1-1/(2h)u_2=const$ and then getting $u_0$ from there. This approximation is of second order and no need for the ghost point! Commented Jan 24, 2013 at 1:05
• I know this approach well, having used it myself many times. This is indeed 2nd order but personally I have found it frequently detrimental when used in stiff time evolving systems and I've never seen it recommended in a text, though I do not have a rigorous argument as to why. Perhaps a good follow up question!
– dmon
Commented Jan 24, 2013 at 9:33
• Actually,if you approximate the derivative with a second order accuracy, then $u_x(0)=O(h^2)$ the value itself would be approximates as $u(0)=$ Commented Jan 25, 2013 at 1:20