# The rate of convergence for finite difference methods for Poisson's equation with piecewise constant data

I am solving the following PDE;

$$\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \rho,$$

where $\rho(0.5,0.5) = 2$ (zero elsewhere), $0\leq x,y\leq1$ and the boundaries will satisfy the condition $u = 0$.

I am solving by finite differences using both a 5 point second order stencil and a 9 point fourth order stencil. It is my understanding that the solution to this PDE is the zero function and this is what I am using to calculate rates of convergence. I have used various error measures including the average difference, max difference and the value at the point $(0.5,0.5)$.

For each error measure I am getting convergence rates of approximately order 2 for both stencils, is there a reason why the fourth order stencil would not converge with order 4?

This is a repeat of a question I asked on MSE - https://math.stackexchange.com/questions/1073945/the-rate-of-convergence-for-finite-difference-methods-for-poissons-equation-wit/1074255#1074255

Although the answers there have helped a lot I would appreciate anyone who could elaborate a little further?

• The answer you got on MSE looks good to me, why post a duplicate question? That said, have you absolutely made sure your implementation is correct? When you want to test a method's implementation, you should start by picking some known functions $u(x,y)$, finding the corresponding b.c. and r.h.s. $\rho$, and comparing with numerical approximations. – Kirill Dec 19 '14 at 18:38
• No I am not completely certain...mainly because it does not converge how I expected! I am trying to test it on a known solution but I am having issues deciding how to deal with boundary points since the nine point stencil extends out two nodes...Duplicate question is because it was suggested to do so... – Josh Greenhalgh Dec 19 '14 at 18:38
• See here about cross-posting; it's a slight breach of etiquette. For boundary points, if you know an exact solution, set the boundary points and the points just beyond the boundary to their exact known values; this gives you the linear equations you need to solve. An alternative is to use a smaller stencil at points next to the boundary, so that no extra ghost points are needed. – Kirill Dec 19 '14 at 18:40

Let's examine the one-dimensional three-point stencil case in detail, because I think it's important to be clear just how this behaviour arises, and what it means to set a point to a certain value in a finite-difference grid when the underlying function is discontinuous. The equation will be $$u''(x) = \rho(x).$$ Instead of using the interval $[-1,1]$ with discontinuity at $\frac12$, I will use the interval $[-\frac12,\frac12]$ with the discontinuity placed at $0$. The grid size will be $h$, and I will only have to consider the interval $[-\frac h2,\frac h2]$ around the grid point $0$.

First, in a finite-difference approximation, we approximate $$u''(0) = \frac{\hat u(-h) - 2\hat u(0) + \hat u(h)}{h^2}$$ and solve $$\frac{\hat u(-h) - 2\hat u(0) + \hat u(h)}{h^2} = \hat \rho(0).$$ (Here the variables with hats are the numerical approximations on the grid to the variables without hats.)

But this is a very bad specification of the problem, because your function $\rho(0) = 2$ at $x=0$, and $\rho(x) = 0$ everywhere else, is discontinuous. In particular, if we shift the grid point $0$ to either side by some tiny amount $\epsilon$, the numerical solution changes entirely and becomes exactly zero. This means that this is a terribly misspecified problem.

We can make sense of it by converting it to an equivalent finite-volume formulation, where it will make much more sense.

In a finite-volume method, we solve $$\int_{-h/2}^{h/2} u''(x)\,dx = \int_{-h/2}^{h/2}\rho(x)\,dx,$$ by picking $u(x)$ to be a suitable approximation to the unknown function. Let's pick, on the interval $[-\frac h2,\frac h2]$ the approximation $$u(x) = \hat u(-h) \phi(x+h) + \hat u(0) \phi(x) + \hat u(h) \phi(x-h),$$ where $\phi(x)$ is the basis function $$\phi(x) = \max\left(1-\frac{|x|}{h}, 0\right),$$ (it's a piecewise linear function that goes from $0$ at $-h$ to $1$ at $0$ to $0$ at $+h$, thus interpolating between grid points.). The approximation $u(x)$ is a weighted sum of three basis functions that look like this:

We can then compute $$\phi''(x) = \frac1h\delta(h-|x|) - \frac2h \delta(x),$$ so that the approximation to the integral is $$\int_{-h/2}^{h/2} u''(x) = \frac{\hat u(-h) -2\hat u(0)+\hat u(h)}{h},$$ and the finite-volume approximation to our equation becomes $$\frac{\hat u(-h) -2\hat u(0) + \hat u(h)}{h} = \int_{-h/2}^{h/2} \rho(x)\,dx = h \hat \rho(0).$$

There are two important things here. First, this is equivalent to the finite-difference formulation in that we end up solving the same equations. Second, the discontinuity in $\rho$ is given a very precise meaning: when we use the value $2$ for $\hat \rho(0)$, we are saying that this is the average value of $\rho$ on the interval $[-h/2,h/2]$: $$2 = \hat\rho(0) = \frac1h \int_{-h/2}^{h/2} \rho(x)\,dx.$$ This interpretation is not available in the finite-difference formulation. It is also not so sensitive to the location of the discontinuity: if the discontinuity were at some small distance $\epsilon$ away from $0$, the average value would be almost the same, but the value at $0$ might be completely different.

But if we say that $2$ is the average value of $\rho$ near the grid point, we can then go back and compute the exact solution of the equation with the right-hand side given by $$\tilde \rho(x) = 2[-h/2 < x < h/2].$$ (We pick a function of our own choice that gives the right average.) In this case, the exact solution will be $$\tilde u(x) = 2 \int_{-h/2}^{h/2} G(x; u)\,du \approx 2h G(x; 0),$$ in terms of the Green function for the Poisson equation. In the two-dimensional case it will be $$\approx 2h^2 G(x, y; 0, 0),$$ as in the other (correct) answer on MSE.

Finally, the outcome of all this is that when you say that you compare your numerical approximate solution with the exact solution $u(x)=0$, this is wrong. The exact solution should not be zero, it should be $$\approx 2h^2 G(x, y; 0,0).$$ It therefore should make perfect sense that the fourth-order solution does not converge to zero with order $4$: it should converge to the correct solution, which is not zero but has magnitude of order $O(h^2)$.

If you do want to get zero as the numerical solution and compare with the mathematical solution, you should set $\hat\rho(0)$ to be the average value of $\rho$, which is $0$, not $2$.

The fact that the exact solution depends on the chosen grid size indicates that this is not a good way to check whether you implemented the method correctly. A very straightforward technique known as the method of manufactured solutions is better for this.

• I really appreciate the work you have put into this answer, it has made things much clearer for me. The method I am using does converge as expected (order 4) for a manufactured solution with no discontinuities - so I am happy that my implementation is correct. In terms of this particular source term - the reason it seems to be converging to the zero function is because of the scaling of the greens function by $h^2$ which in the limit is zero? I still have a bit of work to do to fully understand the convergence I am seeing but I am much closer to being able to do this - thanks again!! – Josh Greenhalgh Dec 19 '14 at 21:13
• @JoshGreenhalgh Upvote answers that you find helpful. Yes, if the exact solution is $O(h^2)$, it will converge to zero as $h\to0$, making it look like you see your numerical method converges to the "right" solution at order 2. – Kirill Dec 19 '14 at 21:17
• I am not yet worthy to upvote...but it is done in spirit! – Josh Greenhalgh Dec 19 '14 at 21:42
• @Kirill Beautiful ! +1 – Sektor Dec 20 '14 at 10:52

Think about what "second order" and "fourth order" approximations actually mean. A first order method can exactly reproduce a first order (linear) function, a second order can reproduce a quadratic, etc.

So say if the exact solution to a given problem is linear - if so, a 100th order method won't give you any benefit over a 1st order method, because both can exactly represent the solution. I suspect this might be a similar case with your zero function solution. Try something with a sinusoidal forcing term or boundary condition.

• I don't see how this answers the question, which is about a very specific right-hand side; the answer shouldn't be "try a different r.h.s.". Can you explain the observed results directly? – Kirill Dec 19 '14 at 18:44
• The mathexchange answer is effectively the same thing, just with some rigor and without my hand-waving. He showed the actual answer to the problem is itself $O(h^2)$. If the analytical answer to a problem is $O(h^2)$ then an $O(h^4)$ approximation will never be any better than an $O(h^2)$ approximation. The answer really is "try a different RHS" or "try different BCs", because the example given won't show a difference between 2nd and higher order schemes. Truncation errors will be identically 0. – Aurelius Dec 19 '14 at 18:55
• The question is "is there a reason why the fourth order stencil would not converge with order 4?", which I interpret as asking about this precise r.h.s. I agree that if the question were "how should I test these methods", the answer would be "try a better r.h.s.". – Kirill Dec 19 '14 at 19:00
• "is there a reason why the fourth order stencil would not converge with order 4?" - yes, and the reason is that the actual solution is order 2, so higher convergence is not possible. One can write out a mathematical proof, but it should be intuitively obvious. If you're trying to approximate a function which is a straight line, a tenth order approximation cannot be more accurate than a first order one, since the first order would be exact. – Aurelius Dec 19 '14 at 19:11
• The fact that if the solution is order 2 a method of order 4 wound not show any improvement is clear to me...however why the 'solution' is order 2 is not! But thanks anyway...my method does converge with order 4 for a problem involving a trigonometric source term just not this one – Josh Greenhalgh Dec 19 '14 at 19:29