Solving Poisson's Equation via Finite Difference¶
Constructing the Finite Difference Matrix¶
Let $\Omega = [0,1]\times[0,1]$ be the unit square domain, with boundary $\partial\Omega$. We will be using a Finite Difference method to obtain a numerial solution to Poisson's equation on $\Omega$, subject to homogeneous Dirichlet boundary conditions:
$$ \begin{align} \nabla^{2}u(\mathbf{x}) &= f(\mathbf{x}) &\quad \mathbf{x}\in\Omega \\ u\vert_{\partial\Omega} &= 0. \end{align} $$The function $u$ is the solution that we seek and $f$ is a known source term.
To begin the Finite Difference implementation, we first need to choose a mesh for $\Omega$, on which we will approximate the solution $u$.
Taking $N\in\mathbb{N}$ we uniformly discretise in the $x$- and $y$- directions with $N$ interior points, so we obtain $N+2$ points in each direction:
$$ \begin{align} x_{i} = \frac{i}{N+1}, &\quad y_{j} = \frac{j}{N+1} &\quad i,j\in\{0,...,N+1\} \end{align} $$and the set of all pairings $(x_{i}, y_{j})$ details the points on the mesh.
It's also useful to introduce the mesh diameter or parameter, $h = \frac{1}{N+1}$; as well as some convenient notation:
$$ \begin{align} u_{i,j} := u\left(x_{i},y_{j}\right) \text{ and } f_{i,j} := f\left(x_{i},y_{j}\right) \end{align} $$Next we write down the difference equations to approximate the Laplacian in each direction:
$$ \begin{align} \frac{\partial^{2} u_{i,j}}{\partial x^{2}} &\approx \frac{u_{i-1,j} -2u_{i,j} + u_{i-1,j} }{h^{2}} \\ \frac{\partial^{2} u_{i,j}}{\partial y^{2}} &\approx \frac{u_{i,j-1} -2u_{i,j} + u_{i,j-1} }{h^{2}} \end{align} $$and so our approximate system of equations, which approximate Poisson's equation on our mesh, is
$$ \begin{align}\label{eq:fd} \frac{1}{h^{2}}\left(-u_{i-1,j} -u_{i,j-1} +4u_{i,j} -u_{i+1,j+1} -u_{i,j+1} \right) &= f\left(x_{i}, y_{j}\right), &\quad \forall i,j\in\{1,...,N\}. \end{align} $$We only need to determine $u_{i,j}$ at interior mesh points; because of the homogeneous Dirichlet conditions we know that $u_{0,j} = u_{N+1,j} = u_{i,0} = u_{i,N+1} = 0$. Therefore the previous system of equations is a system of $N^{2}$ equations in $N^{2}$ unknowns, and we can represent the system in matrix form as
$$ \begin{align} M\mathbf{U} = -\mathbf{F}. \end{align} $$Here $\mathbf{U}$ is a vector of the $u_{i,j}$, $\mathbf{F}$ is a vector of the $f_{i,j}$, and $M$ is the finite difference matrix
$$ M = \frac{1}{h^{2}}\left[\left(\mathrm{tridiag}_{N}(-1,2,-1) \otimes I_{N}\right) + \left(I_{N} \otimes\mathrm{tridiag}_{N}(-1,2,-1)\right)\right] $$where $\otimes$ denotes the Kronecker product and $\mathrm{tridiag}_{N}(-1,2,-1)$ denotes the $N\times N$ matrix with $2$ on the leading diagonal and $-1$ on the first super- and sub- diagonals.
The elements of $\mathbf{U}$ and $\mathbf{F}$ are ordered in $x$ then $y$ - this ordering is automatically compatable with MATLAB's reshape command). That is, $\mathbf{U} = \left(u_{1,1}, u_{2,1}, ..., u_{N,1}, u_{1,2}, u_{2,2}, ... , u_{N-1,N}, u_{N,N} \right)^{\top}$, and the elements of $\mathbf{F}$ are similar.
Finding an approximation to the true solution $u$ then amounts to solving the linear system $M\mathbf{U} = -\mathbf{F}$ for $\mathbf{U}$, which can be done with MATLAB's \ operator.
Solving the Poisson Equation¶
We are interested in the Poisson equation with source term $$ \begin{align} f(x,y) &= -4\pi\sin(2\pi x)\left[\pi(1+(2y-1)^2)\cos(2\pi(y-0.5)^2)+\sin(2\pi(y-0.5)^2)\right] \end{align} $$ which results in the analytic solution $$ \begin{align} u(x,y) &= \sin(2\pi x)\cos(2\pi(y-0.5)^2). \end{align} $$
Examining the Error¶
As $h$ decreases, our approximate solution uApprox converges to the analytic solution $u$.
A common point of interest is how fast
this convergence is; one way we can examine this is by looking at how the error in our solution decays with the mesh size.
We hypothesize that $ || u - $ uApprox $ ||_2 = C h^{\alpha}$ for some constants $C$ and $\alpha$, where $\alpha$ is the rate that we want to estimate.
If we denote $E = || u - $ uApprox $ ||_2 $, then taking logarithms yields $ \log E = \alpha \log h + \log C $, and so if we can obtain a set of values for $E$ and the corresponding mesh gize $h$, we can predict $\alpha$ via a least-squares (regression) fit.