[Math] The rate of convergence for finite difference methods for Poisson’s equation with piecewise constant data

Question

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?

Answer 1

Ian made two excellent points in comments:

1. In the absence of nice smoothness properties in the equation, nice convergence properties of a numerical method can fail.

2. Your numerical method is actually an accurate solver for a Poisson problem where the right side is the indicator function of a small square whose side length has order $h$. I think the convergence rate of the solution to this problem to the solution of the original problem may already be only second order, in which case refining the method can't improve anything.

I will elaborate on the second point. The discrete values of $\rho$ you feed into any finite difference method do not represent pointwise values of $\rho$, but rather its averages on scale $h$. So, feeding in $\rho=2$ at $(0.5,0.5)$ (which I assume is one of the grid points) corresponds to prescribing the Laplacian of $2$ on the square with vertices $(0.5\pm h/2, 0.5\pm h/2)$. The analytic solution to the latter problem is of order $h^2$. Indeed, it is given by the integral of Green's function $G$: $$ u(x,y) = \int_{0.5-h/2}^{0.5+h/2}\int_{0.5-h/2}^{0.5+h/2} 2G(x,y;x',y')\,dx'\,dy' $$ which, at a point $(x,y)\ne (0.5,0.5)$, is asymptotic to $h^2 G(x,y,0.5,0.5)$.

At $(x,y)=(0.5,0.5)$ the analytic solution is actually slightly larger: $h^2\log(1/h)$, but this is probably not apparent in the numerical results.