# How to evaluate the points near/at the boundary when using Richardson extrapolation for improved accuracy of a derivative

If we want to improve the accuracy of our numerical estimation of a derivative, we can use Richardson extrapolation. The method is very beneficial when using a centered difference scheme and the stencil can be derived the following way. Let $$h$$ be the step size. Using Taylor's expansion, we can write:

$$f(x + h) = f(x) + f'(x) h+{f''(x) \over 2}h^2 + {f'''(x) \over 6} h^3 \tag 1$$ $$f(x - h) = f(x) - f'(x) h+{f''(x) \over 2}h^2 - {f'''(x) \over 6} h^3 \tag 2$$

$$f(x + 2h) = f(x) + f'(x) 2h+{f''(x) \over 2}4h^2 + {f'''(x) \over 6} 8h^3 \tag 3$$ $$f(x - 2h) = f(x) - f'(x) 2h+{f''(x) \over 2}4h^2 - {f'''(x) \over 6} 8h^3 \tag 4$$

Subtracting $$(2)$$ from $$(1)$$ and $$(4)$$ from $$(3)$$, and by regrouping some terms yields:

$$f'(x) = {f(x + h) - f(x - h) \over 2h} - {f'''(x) \over 6} h^2 \tag 5$$ $$f'(x) = {f(x + 2h) - f(x - 2h) \over 4h} - {f'''(x) \over 6} 4h^2 \tag 6$$

By reformulating $$(5)$$ as:

$${f'''(x) \over 6} h^2 = {f(x + h) - f(x - h) \over 2h} - f'(x) \tag 7$$

and by inserting this expression into equation $$(6)$$, we get the following estimation of the first derivative:

$$f'(x) = {-f(x+2h) + 8f(x+h) - 8f(x-h)+f(x-2h)\over 12h} \tag 8$$

Similarly, we can derive the estimation for the second derivative:

$$f''(x) = {-f(x+2h) + 16f(x+h)-30f(x) + 16f(x-h)-f(x-2h)\over 12h^2} \tag 9$$

These expressions are fourth-order accurate.

If we are trying to solve a second-order PDE with second-order spatial derivatives using the finite difference method, it seems to me that the boundary points cannot be calculated with fourth-order accuracy (I am hoping that I am wrong). Let's say that we have the boundary condition:

$$f'(x_{max}) = B \tag {10}$$

where $$B$$ is a known constant. From $$(5)$$ and $$(6)$$ we can write:

$$f(x_{max} + h) = f(x_{max} - h) + {f'''(x_{max}) \over 6} 2h^3 + 2hB \tag{11}$$ $$f(x_{max} + 2h) = f(x_{max} - 2h) + {f'''(x_{max}) \over 6} 16h^3 +4hB \tag{12}$$

Inserting these expressions into equation $$(9)$$ gives:

$$f''(x_{max})= { {f'''(x_{max}) \over 6} 16h^3 -30f(x_{max}) + 32f(x_{max} - h) - 2 f(x_{max} - 2h) \over 12h^2 } \tag {13}$$

which has the error term in it!

So it seems the error doesn't cancel out when the boundary condition is applied which leads me to believe I am doing something wrong. Also, evaluating $$f(x_{max} - h)$$ seems to also contain the error term if I just use the boundary condition $$(11)$$ to eliminate the ghost point. On top of that, I don't even understand how to evaluate the point $$f(x_{max} - h)$$ if I have a non-derivative boundary condition.

My question is, how am I supposed to evaluate the boundary points $$f(x_{max})$$ and $$f(x_{max} - h)$$ when I'm dealing with derivative boundary conditions and how am I supposed to evaluate the $$f(x_{max} - h)$$ point when I'm dealing with a non-derivative boundary condition? I am hoping that a way exists that preserves the accuracy of the Richardson extrapolation.

• That is an interesting observation, but by the same argument with substitution of the derivative you can argue that the standard scheme is not even consistent at the boundaries, so I think there's something off about this. The way I see it with the standard scheme is that I have fit a quadratic polynomial passing through my 3 points and taken the second derivative of this. At the boundary I have just prescribed one of those 3 points as fixed (Dirichlet) or through the derivative to its neighbour (Neumann). The polynomial through those is still quadratic though. Commented Feb 11 at 0:34
• One could put the counterargument that the consistency order is based on the Taylor expansion and that I should check what happens with the 0-th, first, second, third, etc. derivatives of my interpolating polynomial and that's probably a valid point. I feel though that this is somehow ignoring the fact that the boundary condition already enters the game at this point. While we do treat this as an interior point to the domain I feel like it's something in between compared to the continuous setting. So I believe that requiring that this be the continuous second derivative is potentially wrong. Commented Feb 11 at 0:41

## 1 Answer

Your question was quite interesting - I haven't seen it addressed in any of the books on FDM that I have read (LeVeque, Morton, Lapidus, Thomas). I had never thought about this, so I tried to come up with some derivations myself and did a small literature review, both with positive outcomes.

The Problem

Let's consider the simplest possible case, of a regular 1D grid $$x_i = ih$$ for $$0\leq i\leq n+1$$ and $$h=\frac{1}{n+1}$$, where $$x_0=0$$ and $$x_{n+1}=1$$ are assumed to be at the boundary. Let us take the following standard discretisations: \begin{align} (d_{xx} f)(x_i) &= \frac{f(x_{i-1}) - 2f(x_i) + f(x_{i+1})}{h^2} + O(h^2) \\ (d_xf)(x_i) &= \frac{f(x_{i})-f(x_{i-1})}{h} + O(h). \end{align} Suppose you have the reflecting boundary condition $$(d_x f)(x_{n+1}) = 0$$, or upon using our discretisation $$f(x_{n+1}) = f(x_n) + O(h^2)$$, where the first derivative is zero because of the boundary condition. If we plug this directly into the discretisation for $$(d_{xx} f)(x_n)$$ we would get: \begin{align} (d_{xx} f)(x_n) &= \frac{f(x_{n-1}) - 2f(x_n) + f(x_{n+1})}{h^2} + O(h^2) \\ &= \frac{f(x_{n-1})-f(x_n) + O(h^2)}{h^2} + O(h^2) \\ &= \frac{f(x_{n-1})-f(x_n)}{h^2} +O(1), \end{align} which formally would be an inconsistent discretisation of $$d_{xx} f$$ at $$x_n$$ since $$\frac{O(h^2)}{h^2} = O(1)$$. I don't really think that is an issue because I would argue that it isn't exactly $$d_{xx}$$ that we are discretising here, since $$x_{n+1}$$ is a boundary point, but is included in the discretisation of an interior point $$x_n$$. Additionally, in virtually all FDM books that I know they just ignore the inconsistency at the boundary. You can find these "inconsistent" discretisations used in papers by reputed mathematicians such as Strang: Functions of Difference Matrices are Toeplitz plus Hankel. I feel like this is similar to how you can have e.g. non-conforming FEM discretisations, which work just fine in practice. The main difference between the formulations in the setting I derive below is that the new derivation enforces the boundary conditions in the continuous setting, formally yielding the desired consistency order, while the standard one enforces the boundary conditions with the discrete FDM operators. So while the standard one may be formally inconsistent, I believe it will still converge to the correct solution (i.e. it's not a practically inconsistent scheme), just more slowly near the boundaries.

My Derivation

Nevertheless, formally the consistency order is $$0$$, so let's try to improve upon that by directly using the fact that $$(d_x f)(x_{n+1}) = 0$$ and consider only the two points $$x_{n-1}$$ and $$x_n$$. Our goal would be to construct a linear combination $$\alpha_{-1}f(x_{n-1}) + \alpha_0 f(x_n)$$ such that it equals $$(d_{xx} f)(x_n) + O(h^p)$$ where $$p$$ must be greater than zero for this to be considered consistent. To do so let us do the Taylor expansion around $$x_n$$: \begin{align} f(x_n) &= f(x_n) \\ f(x_{n-1}) &= f(x_n) - (d_xf)(x_n)h + (d^2_xf)(x_n) \frac{h^2}{2} + O(h^3). \end{align} To use $$(d_x f)(x_{n+1}) = 0$$ directly we will need to expand $$(d_x f)(x_n)$$ around $$x_{n+1}$$: $$(d_x f)(x_n) = (d_x f)(x_{n+1}) - (d^2_x f)(x_{n+1}) h + O(h^2) = - (d^2_x f)(x_{n+1}) h + O(h^2),$$ where I have dropped $$(d_x f)(x_{n+1})$$ since we know $$(d_x f)(x_{n+1}) = 0$$ from the boundary condition. Now I need to reexpress the remaining terms around $$x_n$$: $$(d^2_x f)(x_{n+1}) = (d^2_x f)(x_n) + O(h).$$ Ultimately I have: $$(d_x f)(x_n) = - (d^2_x f)(x_{n}) h + O(h^2).$$ Plugging this into the expansion for $$f(x_{n-1})$$ I get: $$f(x_{n-1}) = f(x_n) + (d^2_xf)(x_n) \left(\frac{h^2}{2}+h^2\right) + O(h^3).$$ I can now set up a system for $$\alpha_{-1}f(x_{n-1}) + \alpha_0 f(x_n) = (d^2_x f)(x_n) + O(h^p)$$: $$\begin{bmatrix} 1 & 1 \\ \frac{3h^2}{2} & 0 \end{bmatrix} \begin{bmatrix} \alpha_{-1} \\ \alpha_0 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \implies \alpha_{-1} = \frac{2}{3h^2}, \,\, \alpha_0 = -\frac{2}{3h^2}, \,\, O(\alpha_{-1}h^3) = O(h).$$ Then the consistent discretisation is: $$(d_{xx} f)(x_n) = 2\frac{f(x_{n-1}) - f(x_n)}{3h^2} + O(h).$$ Typically for the $$l$$-th derivative with a non-symmetric stencil containing $$m$$ terms one can get only a consistency order of $$O(h^{m-l})$$ (for symmetric stencils that can be one order higher). However, since we have one boundary condition, it takes care of one equation allowing us to gain an additional order of consistency.

We can also do a sanity check by using the fact that we are essentially computing the derivative of a polynomial interpolating the points divided by $$l$$ factorial (for the $$l$$-th derivative). A convenient basis is the Newton basis, then we can write the interpolating polynomial as: $$p(x) = f[x_{n-1}] + f[x_{n-1},x_n] (x-x_{n-1}) + f[x_{n-1},x_n,x_{n+1}](x-x_{n-1})(x-x_n)$$. The second derivative is $$(d^2_x p)(x) = f[x_{n-1},x_n,x_{n+1}]$$ and we have $$(d^2_x f)(x_n) = 2(d^2_x p)(x_n) + O(h^2)$$. We can also compute the first derivative at $$x_{n+1}$$: $$(d_x p)(x) = f[x_{n-1}, x_n] + f[x_{n-1}, x_n, x_{n+1}]((x-x_n) + (x-x_{n-1})),$$ and require that $$(d_{x} f)(x_{n+1}) = (d_x p)(x_{n+1}) + O(h) = 0$$, which yields: $$f[x_{n-1},x_n,x_{n+1}] = -\frac{f[x_{n-1}, x_n]}{3h} \implies (d^2_xf)(x_n) = 2\frac{f(x_{n-1})-f(x_n)}{3h^2} + O(h).$$

The result is the same. If instead of reflecting boundary conditions you set $$(d_x f)(x_{n+1}) = B$$, then this changes only slightly and you get: $$(d^2_xf)(x_n) = \frac{2B}{3h} + 2\frac{f(x_{n-1})-f(x_n)}{3h^2} + O(h),$$ and you can move the term with $$B$$ to the right-hand side in an ODE/PDE.

If you wanted to have Dirichlet boundary conditions $$f(x_{n+1}) = B$$ then it is even simpler as there is no truncation error whatsoever. You just replace $$f(x_{n+1})$$ with $$B$$ in the stencil and you do not introduce any error: $$(d^2_x f)(x_n) = \frac{f(x_{n-1})-2f(x_n) + f(x_{n+1})}{h^2} + O(h^2) = \frac{f(x_{n-1})-2f(x_n)}{h^2} + \frac{B}{h^2} + O(h^2).$$

Neumann Boundary with $$O(h^2)$$ for $$d^2_x f$$

If you wanted to recover $$O(h^2)$$ at the Neumann boundary you would have to provide three interior points in your stencil. That is, consider $$x_{n-2}, x_{n-1}, x_n$$. For brevity will use $$g_i$$ to refer to $$g(x_i)$$. Let us try to approximate the second derivative as $$\alpha_{-2}f_{n-2} + \alpha_{-1} f_{n-1} + \alpha_0 f_n$$. To do so you can do the Taylor expansion again (now until $$O(h^4)$$): \begin{align} f_n &= f_n \\ f_{n-1} &= f_n - (d_x f)_n h + (d^2_x f)_n \frac{h^2}{2} - (d^3_x f)_n \frac{h^3}{6} + O(h^4) \\ f_{n-2} &= f_n - (d_x f)_n 2h + (d^2_x f)_n \frac{4h^2}{2} - (d^3_x f)_n \frac{8h^3}{6} + O(h^4). \end{align} We know $$(d_x f)_{n+1} = B$$, so we can expand $$(d_x f)_n$$ around $$x_{n+1}$$: \begin{align} (d_x f)_{n} &= (d_x f)_{n+1} - (d^2_x f)_{n+1} h + (d^3_x f)_{n+1} \frac{h^2}{2} + O(h^3)\\ &= B - (d^2_x f)_{n+1} h + (d^3_x f)_{n+1} \frac{h^2}{2} + O(h^3). \end{align} Now we need to expand the other terms around $$x_n$$: \begin{align} (d^2_x f)_{n+1} &= (d^2_x f)_n + (d^3_x f)_n h + O(h^2) \\ (d^3_x f)_{n+1} &= (d^3_x f)_n + O(h). \end{align} Plugging those in $$(d_x f)_{n}$$ you get: \begin{align} (d_x f)_{n} &= B - (d^2_x f)_{n} h + (d^3_x f)_{n} \left(\frac{h^2}{2}-h^2\right) + O(h^3). \end{align} Plugging this back into $$f_{n-1}$$ and $$f_{n-2}$$ yields: \begin{align} f_n &= f_n \\ f_{n-1} &= f_n - hB + (d^2_x f)_n \left(\frac{h^2}{2}+h^2\right) + (d^3_x f)_n \left(\frac{h^3}{2}-\frac{h^3}{6}\right) + O(h^4) \\ f_{n-2} &= f_n - 2hB + (d^2_x f)_n \left(\frac{4h^2}{2}+h^2\right) + (d^3_x f)_n \left(\frac{h^3}{2}-\frac{8h^3}{6}\right) + O(h^4). \end{align} The term $$\alpha_{-2}2hB + \alpha_{-1}hB$$ would go to the right-hand side so we can ignore it for now. By requiring that $$\alpha_{-2} f_{n-2} + \alpha_{-1} f_{n-1} + \alpha_0 f_n = (d^2_x f)_n + O(h^p)$$ you get: $$\begin{bmatrix} 1 & 1 & 1 \\ \frac{6h^2}{2} & \frac{3h^2}{2} & 0 \\ -\frac{5h^3}{6} & \frac{2h^3}{6} & 0\end{bmatrix}\begin{bmatrix} \alpha_{-2} \\ \alpha_{-1} \\ \alpha_0 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} \implies \alpha_{-2} = \frac{4}{27h^2}, \alpha_{-1} = \frac{10}{27h^2}, \alpha_{0} = -\frac{14}{27h^2}.$$ We get a consistency order of two as desired $$O((\alpha_{-2}+\alpha_{-1})h^4) = O(h^2)$$ for the stencil: $$(d^2_x f)_n = \frac{4f_{n-2} + 10f_{n-1} - 14 f_n}{27h^2} - \frac{18}{27 h}B.$$ The last term ought to be moved to the right-hand side when you're solving your ODE or PDE.

I haven't done a sanity check since this is rather cumbersome to do with Newton, so there may be typos.

Your Scheme

Your scheme considered 5 points. There are two cases: the first is if the stencil includes one ghost cell, the second is if the stencil includes two ghost cells. Let's consider the second case. We have the points $$x_{n-2}, x_{n-1}, x_{n}$$ and we want to compute $$(d_{xx} f)_n$$. But then this is exactly the extended case I wrote before and the scheme I derived was of consistency $$O(h^2)$$.

Generally, if you have $$m$$ non-symmetrically distributed points around the node you are trying to approximate the $$l$$-th derivative at, you will get a consistency of $$O(h^{m-l})$$. However, if you are near the boundary and you have $$q = \lfloor \frac{l}{2}\rfloor$$ boundary conditions you can get $$O(h^{n-l+q})$$ as illustrated here. You are of course free to make the stencil larger, e.g. with $$5$$ points as in the interior (although here they will go to the left: $$x_{n-4},\ldots,x_n$$), then in your case you will get $$p=5-2+1 = 4$$. So you can regain the order of consistency by allowing as many points as for your discretisation in the interior. It's a bit trickier for the stencil that we get for the node at $$x_{n-1}$$ with points $$x_{n-3}, x_{n-2}, x_{n-1}, x_n$$. I have left it out since it doesn't really provide any new insights - the same result that I stated above holds, it's just that the computation of the stencil is more cumbersome (it results in a $$4\times 4$$ system).

Some Literature

I found the following thesis: Analysis of Boundary and Interface Closures for Finite Difference Methods for the Wave Equation. On page 14 they reference this effect, namely: $$(D_1 v)_j = \begin{cases} (d_x u)(x_j)+O(h^{2p}) &\text{if } j = n_p,\ldots, n-n_p,\\ (d_xu)(x_j) + O(h^p) &\text{else}.\end{cases},$$ where $$n_p$$ are the number of stencils that include ghost nodes. They also cite the following reference for some such stencils: Summation by Parts for Finite Difference Approximations for d/dx. They are considering specific finite differences with additional constraints though, so I am not sure how relevant this is beyond the derivation I provided. I am just linking it since it confirms the idea that the consistency near the boundary decreases.

• What a fantastic answer. Thank you very much for the explanation and the references. It never occurred to me to calculate the boundary nodes without changing the number of nodes for calculation, but shifting the stencil to the right (or left). That's a very nice insight. I think even the WENO method uses this approach to avoid calculations at discontinuity. It's good to know that the order of accuracy can be preserved in the end, even though some mathematics must be done before the implementation of the numerical solver. Commented Feb 15 at 14:51
• @NikolaRistic It was interesting enough for me too - I never thought about this, although judging by the fact that I haven't seen it used much elsewhere it's probably not really problematic. Using non-symmetric stencils is pretty standard e.g. for transport equations - there if you take the central difference you would get a singular problem. Also as you noted, you can use non-symmetric stencils to account for discontinuities. E.g. forward, backward, and central differences converge to the same thing for a differentiable function, but differ if the function is non-differentiable at the point. Commented Feb 15 at 16:34