Elliptic PDE: Proving that a 2nd order accurate discretization of the 2nd derivative of the unknown is 2nd order accurate for the unknown itself

Question

For an ODE:

$\frac{dy}{dt}=f(y(t),t)$

The Euler Explicit scheme reads:

$y_{n+1}=y_{n}+\Delta tf_n$

and it can be easily shown with a Taylor expension that:

$y_{n+1}=y_{n}+\Delta t \frac{dy}{dt}|_n+\Delta t^2 \frac{d^2y}{dt^2}|_n$

and therefore by comparing the scheme with the Taylor expansion we have:

$LTE= O(\Delta t^2)$

i.e. the local truncation error LTE for Euler explicit scales with $\Delta t^2$. Since we repeat this for $n$ times, with $n$ the number of grid points, the global truncation error is $n\Delta t^2$. Since $n=L/\Delta t$, with L the domain length, then the global truncation error scales with Δt, so Euler explicit is first order accurate. For an elliptic PDE:

$\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}=0$

if I use the central difference scheme:

$\frac{\partial^2\phi}{\partial x^2}= \frac{\phi_{i-1,j}-2\phi_{i,j}+\phi_{i+1,j}}{dx^2}+O(\Delta x^2)$

the error is second order as it can be shown by Taylor expanding $\phi_{i+1}$ and $\phi_{i-1}$ in the neighborhood of $\phi_{i}$. Similar reasoning in $y$ where:

$\frac{\partial^2\phi}{\partial y^2}= \frac{\phi_{i,j-1}-2\phi_{i,j}+\phi_{i,j+1}}{dy^2}+O(\Delta y^2)$

And here is were I am lost. The error in the derivative is second order, but what about the error in the function $\phi$ I am looking for? If I replace the discretization into the equation (let's assume $\Delta y = \Delta x$ for now) we have:

$\frac{\phi_{i,j-1}-2\phi_{i,j}+\phi_{i,j+1}}{dx^2}+ \frac{\phi_{i-$1,j}-2\phi_{i,j}+\phi_{i+1,j}}{dx^2}+2O(\Delta x^2)=0$

which, multiplying by $\Delta x^2$ reads:

$\phi_{i,j-1}-4\phi_{i,j}+\phi_{i,j+1}+ \phi_{i,i-1,j}+\phi_{i+1,j}+2O(\Delta x^4)=0$

So $O(\Delta x^4)$ shows up. From here, how can I prove that $\phi$ is obtained with second order accuracy? I guess the Local error for $\phi_{j}$ is 4th order (if all the 4 neighbours were known exactly) , but the global error is second error, I just do not know how to prove it. Any suggestion?

Answer 1

Consider as a simple case the 1D Poisson problem (it should be simple to extend to the 2D case you describe): $$ \tag{$*$}\label{eq:poisson} \frac{d^2 \phi}{dx^2}(x) = f(x), \quad 0<x<1, \qquad \phi(0) = \phi(1) = 0, $$ where $f$ is given. We approximate $$ \frac{d^2 \phi}{dx^2}(x_i) = \frac{\phi_{i-1} - 2\phi_i + \phi_{i+1}}{h^2} + \mathcal{O}(h^2). $$ Let $\phi(x)$ denote the exact solution to the problem $(*)$. Suppose $\tilde{\phi}$ is the numerical solution obtained by solving the system of equations $$ \frac{\tilde\phi_{i-1} - 2\tilde\phi_i + \tilde\phi_{i+1}}{h^2} = f_i. $$ We now write everything in terms of linear algebra. Let $\boldsymbol\phi$ denote the vector of values $\phi(x_i)$, let $\tilde{\boldsymbol\phi}$ denote the vector $\tilde\phi_i$, and let $\boldsymbol f$ denote the vector $f(x_i)$. Finally, let $A$ denote the matrix corresponding to the finite-difference Laplacian. For the discrete approximate solution, we have $$ A \tilde{\boldsymbol\phi} = \boldsymbol f. $$ Using the truncation error for the finite-difference approximation, we have for the exact solution $$ A \boldsymbol\phi = \boldsymbol f + \boldsymbol \tau, \qquad \tau_i = \mathcal{O}(h^2). $$ We want to measure the error in our approximate solution, which is given by $\| \tilde{\boldsymbol\phi} - \boldsymbol\phi \|$ in some norm. We see that the difference $\boldsymbol\varepsilon = \tilde{\boldsymbol\phi} - \boldsymbol\phi$ satisfies $$ A \boldsymbol\varepsilon = A (\tilde{\boldsymbol\phi} -\boldsymbol\phi)= -\boldsymbol \tau $$ and so $$ \boldsymbol\varepsilon = -A^{-1}\boldsymbol \tau. $$ In particular $$ \| \boldsymbol\varepsilon \| \leq \|A^{-1}\| \|\boldsymbol\tau\| $$ where $\|A^{-1}\|$ is the corresponding matrix norm. If we can bound the norm of $A^{-1}$ by a constant, then we can conclude that $\|\boldsymbol\varepsilon\| = \mathcal{O}(h^2)$. If we use the 2-norm and using that $A$ is symmetric, we see that $$ \| A^{-1} \|_2 = \frac{1}{\min |\lambda|} $$ where $\lambda$ is an eigenvalue of $A$. Using an explicit expression for the eigenvalues of $A$, the smallest in magnitude is given by $\lambda = -\pi^2 + \mathcal{O}(h^2)$ and the result follows.