Finite Difference Methods for the Laplacian Equation

John S Butler john.s.butler@tudublin.ie Course Notes Github

Overview

This notebook will focus on numerically approximating a homogenous second order Poisson Equation which is the Laplacian Equation.

The Differential Equation

The general two dimensional Poisson Equation is of the form:

(648)$$\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial x^2}=f(x,y),(x,y)\in \mathrm{\Omega }=(0,1)\times (0,1),$$

with boundary conditions

(649)$$U(x,y)=g(x,y),(x,y)\in \delta \mathrm{\Omega }\text{ - boundary}.$$

Homogenous Poisson Equation

This notebook will implement a finite difference scheme to approximate the homogenous form of the Poisson Equation $f(x,y)=0$:

(650)$$\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial x^2}=0.$$

with the Boundary Conditions:

(651)$$u(x,0)=\sin (2\pi x),0\le x\le 1,\text{ lower},$$

(652)$$u(x,1)=\sin (2\pi x),0\le x\le 1,\text{ upper},$$

(653)$$u(0,y)=2\sin (2\pi y),0\le y\le 1,\text{ left},$$

(654)$$u(1,y)=2\sin (2\pi y),0\le y\le 1,\text{ right}.$$

# LIBRARY
# vector manipulation
import numpy as np
# math functions
import math 

# THIS IS FOR PLOTTING
%matplotlib inline
import matplotlib.pyplot as plt # side-stepping mpl backend
import warnings
warnings.filterwarnings("ignore")
from IPython.display import HTML
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt

Discete Grid

The region $\mathrm{\Omega }=(0,1)\times (0,1)$ is discretised into a uniform mesh ${\mathrm{\Omega }}_h$. In the $x$ and $y$ directions into $N$ steps giving a stepsize of

(655)$$h=\frac{1-0}{N},$$

resulting in

(656)$$x[i]=0+ih,i=0,1,...,N,\$$$

and

(657)$$x[j]=0+jh,j=0,1,...,N,$$

The Figure below shows the discrete grid points for $N=10$, the known boundary conditions (green), and the unknown values (red) of the Poisson Equation.

N=10
h=1/N
x=np.arange(0,1.0001,h)
y=np.arange(0,1.0001,h)
X, Y = np.meshgrid(x, y)
fig = plt.figure()
plt.plot(x[1],y[1],'ro',label='unknown');
plt.plot(X,Y,'ro');
plt.plot(np.ones(N+1),y,'go',label='Boundary Condition');
plt.plot(np.zeros(N+1),y,'go');
plt.plot(x,np.zeros(N+1),'go');
plt.plot(x, np.ones(N+1),'go');
plt.xlim((-0.1,1.1))
plt.ylim((-0.1,1.1))
plt.xlabel('x')
plt.ylabel('y')
plt.gca().set_aspect('equal', adjustable='box')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.title(r'Discrete Grid $\Omega_h,$ h= %s'%(h),fontsize=24,y=1.08)
plt.show();

Boundary Conditions

The discrete boundary conditions are

(658)$$w_{i0}=w[i,0]=\sin (2\pi x[i]),\text{ for }i=0,...,10,\text{ upper},$$

(659)$$w_{iN}=w[i,N]=\sin (2\pi x[i]),\text{ for }i=0,...,10,\text{ lower},$$

(660)$$w_{0j}=w[0,j]=2\sin (2\pi y[j]),\text{ for }j=0,...,10,\text{ left},$$

(661)$$w_{Nj}=w[N,j]=2\sin (2\pi y[j]),\text{ for }i=0,...,10,\text{ right}.$$

The Figure below plots the boundary values of $w[i,j]$.

w=np.zeros((N+1,N+1))

for i in range (0,N):
        w[i,0]=np.sin(2*np.pi*x[i]) #left Boundary
        w[i,N]=np.sin(2*np.pi*x[i]) #Right Boundary

for j in range (0,N):
        w[0,j]=2*np.sin(2*np.pi*y[j]) #Lower Boundary
        w[N,j]=2*np.sin(2*np.pi*y[j]) #Upper Boundary

        
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot a basic wireframe.
ax.plot_wireframe(X, Y, w,color='r', rstride=10, cstride=10)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('w')
plt.title(r'Boundary Values',fontsize=24,y=1.08)
plt.show()

Numerical Method

The Poisson Equation is discretised using ${\delta }_x^2$ is the central difference approximation of the second derivative in the $x$ direction

(662)$${\delta }_x^2=\frac{1}{h^2}(w_{i+1j}-2w_{ij}+w_{i-1j}),$$

and ${\delta }_y^2$ is the central difference approximation of the second derivative in the $y$ direction

(663)$${\delta }_y^2=\frac{1}{h^2}(w_{ij+1}-2w_{ij}+w_{ij-1}).$$

The gives the Poisson Difference Equation,

(664)$$-({\delta }_x^2{w}_{ij}+{\delta }_y^2{w}_{ij})=f_{ij}(x_i,y_j)\in {\mathrm{\Omega }}_h,$$

(665)$$w_{ij}=g_{ij}(x_i,y_j)\in \partial {\mathrm{\Omega }}_h,$$

where $w_{i}j$ is the numerical approximation of $U$ at $x_i$ and $y_j$. Expanding the the Poisson Difference Equation gives the five point method,

(666)$$-(w_{i-1j}+w_{ij-1}-4w_{ij}+w_{ij+1}+w_{i+1j})=h^2{f}_{ij}$$

for $i=1,...,N-1$ and $j=1,...,N-1$ which can be written

(667)$${\mathrm{\nabla }}_h^2{w}_{ij}=f_{ij}.$$

Matrix form

This can be written as a system of $(N-1)\times (N-1)$ equations can be arranged in matrix form

(668)$$A\mathbf{w}=\mathbf{r},$$

where $A$ is an $(N-1{)}^2\times (N-1{)}^2$ matrix made up of the following block tridiagonal structure

(669)$$\left(\begin{array}{ccccccc}T& I& 0& 0& .& .& .\\ I& T& I& 0& 0& .& .\\ .& .& .& .& .& .& .\\ .& .& .& 0& I& T& I\\ .& .& .& .& 0& I& T\end{array}\right),$$

where $I$ denotes an $N-1\times N-1$ identity matrix and $T$ is the tridiagonal matrix of the form:

(670)$$T=\left(\begin{array}{ccccccc}-4& 1& 0& 0& .& .& .\\ 1& -4& 1& 0& 0& .& .\\ .& .& .& .& .& .& .\\ .& .& .& 0& 1& -4& 1\\ .& .& .& .& 0& 1& -4\end{array}\right).$$

The plot below shows the matrix $A$ and its inverse $A^{-1}$ as a colourplot.

N2=(N-1)*(N-1)
A=np.zeros((N2,N2))
## Diagonal            
for i in range (0,N-1):
    for j in range (0,N-1):           
        A[i+(N-1)*j,i+(N-1)*j]=-4

# LOWER DIAGONAL        
for i in range (1,N-1):
    for j in range (0,N-1):           
        A[i+(N-1)*j,i+(N-1)*j-1]=1   
# UPPPER DIAGONAL        
for i in range (0,N-2):
    for j in range (0,N-1):           
        A[i+(N-1)*j,i+(N-1)*j+1]=1   

# LOWER IDENTITY MATRIX
for i in range (0,N-1):
    for j in range (1,N-1):           
        A[i+(N-1)*j,i+(N-1)*(j-1)]=1        
        
        
# UPPER IDENTITY MATRIX
for i in range (0,N-1):
    for j in range (0,N-2):           
        A[i+(N-1)*j,i+(N-1)*(j+1)]=1
Ainv=np.linalg.inv(A)   
fig = plt.figure(figsize=(12,4));
plt.subplot(121)
plt.imshow(A,interpolation='none');
clb=plt.colorbar();
clb.set_label('Matrix elements values');
plt.title('Matrix A ',fontsize=24)
plt.subplot(122)
plt.imshow(Ainv,interpolation='none');
clb=plt.colorbar();
clb.set_label('Matrix elements values');
plt.title(r'Matrix $A^{-1}$ ',fontsize=24)

fig.tight_layout()
plt.show();

The vector $\mathbf{w}$ is of length $(N-1)\times (N-1)$ which made up of $N-1$ subvectors ${\mathbf{w}}_j$ of length $N-1$ of the form

(671)$${\mathbf{w}}_j=\left(\begin{array}{c}{w}_{1j}\\ w_{2j}\\ .\\ .\\ w_{N-2j}\\ w_{N-1j}\end{array}\right).$$

The vector $\mathbf{r}$ is of length $(N-1)\times (N-1)$ which made up of $N-1$ subvectors of the form ${\mathbf{r}}_j=-h^2{\mathbf{f}}_j-{\mathbf{bx}}_j-{\mathbf{by}}_j$, where ${\mathbf{bx}}_j$ is the vector of left and right boundary conditions

(672)$${\mathbf{bx}}_j=\left(\begin{array}{c}{w}_{0j}\\ 0\\ .\\ .\\ 0\\ w_{Nj}\end{array}\right),$$

for $j=1,..,N-1$, where ${\mathbf{by}}_j$ is the vector of the lower boundary condition for $j=1$,

(673)$${\mathbf{by}}_1=\left(\begin{array}{c}{w}_{10}\\ w_{20}\\ .\\ .\\ w_{N-20}\\ w_{N-10}\end{array}\right),$$

upper boundary condition for $j=N-1$

(674)$${\mathbf{by}}_{N-1}=\left(\begin{array}{c}{w}_{1N}\\ w_{2N}\\ .\\ .\\ w_{N-2N}\\ w_{N-1N}\end{array}\right),$$

for $j=2,...,N-2$

(675)$${\mathbf{by}}_j=0,$$

and

(676)$${\mathbf{f}}_j=-\left(\begin{array}{c}0\\ 0\\ .\\ .\\ 0\\ 0\end{array}\right)$$

for $j=1,...,N-1$.

r=np.zeros(N2)

# vector r      
for i in range (0,N-1):
    for j in range (0,N-1):           
        r[i+(N-1)*j]=-h*h*0      
# Boundary        
b_bottom_top=np.zeros(N2)
for i in range (0,N-1):
    b_bottom_top[i]=np.sin(2*np.pi*x[i+1]) #Bottom Boundary
    b_bottom_top[i+(N-1)*(N-2)]=np.sin(2*np.pi*x[i+1])# Top Boundary
      
b_left_right=np.zeros(N2)
for j in range (0,N-1):
    b_left_right[(N-1)*j]=2*np.sin(2*np.pi*y[j+1]) # Left Boundary
    b_left_right[N-2+(N-1)*j]=2*np.sin(2*np.pi*y[j+1])# Right Boundary
    
b=b_left_right+b_bottom_top

Results

To solve the system for $\mathbf{w}$ invert the matrix $A$

(677)$$A\mathbf{w}=\mathbf{r},$$

such that

(678)$$\mathbf{w}=A^{-1}\mathbf{r}.$$

Lastly, as $\mathbf{w}$ is in vector it has to be reshaped into grid form to plot.

The figure below shows the numerical approximation of the homogeneous Equation.

C=np.dot(Ainv,r-b)
w[1:N,1:N]=C.reshape((N-1,N-1))

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d');
# Plot a basic wireframe.
ax.plot_wireframe(X, Y, w,color='r');
ax.set_xlabel('x');
ax.set_ylabel('y');
ax.set_zlabel('w');
plt.title(r'Numerical Approximation of the Poisson Equation',fontsize=24,y=1.08);
plt.show();

Consistency and Convergence

We now ask how well the grid function determined by the five point scheme approximates the exact solution of the Poisson problem.

Consistency

Consitency (Definition)

Let

(679)$${\mathrm{\nabla }}_h^2(\phi )=-({\phi }_{i-1j}+{\phi }_{ij-1}-4{\phi }_{ij}+{\phi }_{ij+1}+{\phi }_{i+1j})$$

denote the finite difference approximation associated with the grid ${\mathrm{\Omega }}_h$ having the mesh size $h$, to a partial differential operator

(680)$${\mathrm{\nabla }}^2(\phi )=\frac{{\partial }^2\phi }{\partial x^2}+\frac{{\partial }^2\phi }{\partial y^2}$$

defined on a simply connected, open set $\mathrm{\Omega }\subset R^2$. For a given function $\phi \in C^{\mathrm{\infty }}(\mathrm{\Omega })$, the truncation error of ${\mathrm{\nabla }}_h^2$ is

(681)$${\tau }_h(\mathbf{x})=({\mathrm{\nabla }}^2-{\mathrm{\nabla }}_h^2)\phi (\mathbf{x})$$

The approximation ${\mathrm{\nabla }}_h^2$ is consistent with ${\mathrm{\nabla }}^2$ if

(682)$$\underset{h\to 0}{lim}{\tau }_h(\mathbf{x})=0,$$

for all $\mathbf{x}\in D$ and all $\phi \in C^{\mathrm{\infty }}(\mathrm{\Omega })$. The approximation is consistent to order $p$ if ${\tau }_h(\mathbf{x})=O(h^p)$.

In other words a method is consistent with the differential equation it is approximating.

Proof of Consistency

The five-point difference analog ${\mathrm{\nabla }}_h^2$ is consistent to order 2 with ${\mathrm{\nabla }}^2$.

Proof

Pick $\phi \in C^{\mathrm{\infty }}(D)$, and let $(x,y)\in \mathrm{\Omega }$ be a point such that $(x\pm h,y),(x,y\pm h)\in \mathrm{\Omega }\bigcup \partial \mathrm{\Omega }$. By the Taylor Theorem

$$\begin{array}{rcl}\phi (x\pm h,y)& =& \phi (x,y)\pm h\frac{\partial \phi }{\partial x}(x,y)+\frac{h^2}{2!}\frac{{\partial }^2\phi }{\partial x^2}(x,y)\pm \frac{h^3}{3!}\frac{{\partial }^3\phi }{\partial x^3}(x,y)+\frac{h^4}{4!}\frac{{\partial }^4\phi }{\partial x^4}({\zeta }^{\pm },y)\end{array}$$

where ${\zeta }^{\pm }\in (x-h,x+h)$. Adding this pair of equation together and rearranging , we get

(683)$$\frac{1}{h^2}[\phi (x+h,y)-2\phi (x,y)+\phi (x-h,y)]-\frac{{\partial }^2\phi }{\partial x^2}(x,y)=\frac{h^2}{4!}[\frac{{\partial }^4\phi }{\partial x^4}({\zeta }^{+},y)+\frac{{\partial }^4\phi }{\partial x^4}({\zeta }^{-},y)]$$

By the intermediate value theorem

(684)$$[\frac{{\partial }^4\phi }{\partial x^4}({\zeta }^{+},y)+\frac{{\partial }^4\phi }{\partial x^4}({\zeta }^{-},y)]=2\frac{{\partial }^4\phi }{\partial x^4}(\zeta ,y),$$

for some $\zeta \in (x-h,x+h)$. Therefore,

(685)$${\delta }_x^2(x,y)=\frac{{\partial }^2\phi }{\partial x^2}(x,y)+\frac{h^2}{2!}\frac{{\partial }^4\phi }{\partial x^4}(\zeta ,y)$$

Similar reasoning shows that

(686)$${\delta }_y^2(x,y)=\frac{{\partial }^2\phi }{\partial y^2}(x,y)+\frac{h^2}{2!}\frac{{\partial }^4\phi }{\partial y^4}(x,\eta )$$

for some $\eta \in (y-h,y+h)$. We conclude that ${\tau }_h(x,y)=(\mathrm{\nabla }-{\mathrm{\nabla }}_h)\phi (x,y)=O(h^2).$

Convergence

Definition

Let ${\mathrm{\nabla }}_h^{2}w({\mathbf{x}}_j)=f({\mathbf{x}}_j)$ be a finite difference approximation, defined on a grid mesh size $h$, to a PDE ${\mathrm{\nabla }}^{2}U(\mathbf{x})=f(\mathbf{x})$ on a simply connected set $D\subset R^n$. Assume that $w(x,y)=U(x,y)$ at all points $(x,y)$ on the boundary $\partial \mathrm{\Omega }$. The finite difference scheme converges (or is convergent) if

(687)$$\underset{j}{max}|U({\mathbf{x}}_j)-w({\mathbf{x}}_j)|\to 0\text{ as }h\to 0.$$

Theorem (DISCRETE MAXIMUM PRINCIPLE).

If ${\mathrm{\nabla }}_h^2{V}_{ij}\ge 0$ for all points $(x_i,y_j)\in {\mathrm{\Omega }}_h$, then

(688)$$\underset{(x_i,y_j)\in {\mathrm{\Omega }}_h}{max}{V}_{ij}\le \underset{(x_i,y_j)\in \partial {\mathrm{\Omega }}_h}{max}{V}_{ij},$$

If ${\mathrm{\nabla }}_h^2{V}_{ij}\le 0$ for all points $(x_i,y_j)\in {\mathrm{\Omega }}_h$, then

(689)$$\underset{(x_i,y_j)\in {\mathrm{\Omega }}_h}{min}{V}_{ij}\ge \underset{(x_i,y_j)\in \partial {\mathrm{\Omega }}_h}{min}{V}_{ij}.$$

Propositions

1. The zero grid function for which $U_{ij}=0$ for all $(x_i,y_j)\in {\mathrm{\Omega }}_h\bigcup \partial {\mathrm{\Omega }}_h$ is the only solution to the finite difference problem

(690)$${\mathrm{\nabla }}_h^2{U}_{ij}=0\text{ for }(x_i,y_j)\in {\mathrm{\Omega }}_h,$$

(691)$$U_{ij}=0\text{ for }(x_i,y_j)\in \partial {\mathrm{\Omega }}_h.$$

2. For prescribed grid functions $f_{ij}$ and $g_{ij}$, there exists a unique solution to the problem

(692)$${\mathrm{\nabla }}_h^2{U}_{ij}=f_{ij}\text{ for }(x_i,y_j)\in {\mathrm{\Omega }}_h,$$

(693)$$U_{ij}=g_{ij}\text{ for }(x_i,y_j)\in \partial {\mathrm{\Omega }}_h.$$

Definition

For any grid function $V:{\mathrm{\Omega }}_h\bigcup \partial {\mathrm{\Omega }}_h\to R$,

(694)$$||V|{|}_{\mathrm{\Omega }}=\underset{(x_i,y_j)\in {\mathrm{\Omega }}_h}{max}|V_{ij}|,$$

(695)$$||V|{|}_{\partial \mathrm{\Omega }}=\underset{(x_i,y_j)\in \partial {\mathrm{\Omega }}_h}{max}|V_{ij}|.$$

Lemma

If the grid function $V:{\mathrm{\Omega }}_h\bigcup \partial {\mathrm{\Omega }}_h\to R$ satisfies the boundary condition $V_{ij}=0$ for $(x_i,y_j)\in \partial {\mathrm{\Omega }}_h$, then

(696)$$||V_{|}{|}_{\mathrm{\Omega }}\le \frac{1}{8}||{\mathrm{\nabla }}_h^{2}V|{|}_{\mathrm{\Omega }}.$$

Given these Lemmas and Propositions, we can now prove that the solution to the five point scheme ${\mathrm{\nabla }}_h^2$ is convergent to the exact solution of the Poisson Equation ${\mathrm{\nabla }}^2$.

Convergence Theorem

Let $U$ be a solution to the Poisson equation and let $w$ be the grid function that satisfies the discrete analog

(697)$$-{\mathrm{\nabla }}_h^2{w}_{ij}=f_{ij}\text{ for }(x_i,y_j)\in {\mathrm{\Omega }}_h,$$

(698)$$w_{ij}=g_{ij}\text{ for }(x_i,y_j)\in \partial {\mathrm{\Omega }}_h.$$

Then there exists a positive constant $K$ such that

(699)$$||U-w|{|}_{\mathrm{\Omega }}\le KM{h}^2,$$

where

(700)$$M=\{{\left|\left|\frac{{\partial }^{4}U}{\partial x^4}\right|\right|}_{\mathrm{\infty }},{\left|\left|\frac{{\partial }^{4}U}{\partial y^4}\right|\right|}_{\mathrm{\infty }}\}$$

Proof

The statement of the theorem assumes that $U\in C^4(\overline{\mathrm{\Omega }})$. This assumption holds if $f$ and $g$ are smooth enough. \begin{proof} Following from the proof of the Proposition we have

(701)$$({\mathrm{\nabla }}_h^2-{\mathrm{\nabla }}^2)U_{ij}=\frac{h^2}{12}[\frac{{\partial }^{4}U}{\partial x^4}({\zeta }_i,y_j)+\frac{{\partial }^{4}U}{\partial y^4}(x_i,{\eta }_j)],$$

for some $\zeta \in (x_{i-1},x_{i+1})$ and ${\eta }_j\in (y_{j-1},y_{j+1})$. Therefore,

(702)$$-{\mathrm{\nabla }}_h^2{U}_{ij}=f_{ij}-\frac{h^2}{12}[\frac{{\partial }^{4}U}{\partial x^4}({\zeta }_i,y_j)+\frac{{\partial }^{4}U}{\partial y^4}(x_i,{\eta }_j)].$$

If we subtract from this the identity equation $-{\mathrm{\nabla }}_h^2{w}_{ij}=f_{ij}$ and note that $U-w$ vanishes on $\partial {\mathrm{\Omega }}_h$, we find that

(703)$${\mathrm{\nabla }}_h^2(U_{ij}-w_{ij})=\frac{h^2}{12}[\frac{{\partial }^{4}U}{\partial x^4}({\zeta }_i,y_j)+\frac{{\partial }^{4}U}{\partial y^4}(x_i,{\eta }_j)].$$

It follows that

(704)$$||U-w|{|}_{\mathrm{\Omega }}\le \frac{1}{8}||{\mathrm{\nabla }}_h^2(U-w)|{|}_{\mathrm{\Omega }}\le KM{h}^2.$$