40 Poisson’s Equation in 2d

A possible topic for a presentation/poster/paper. Fill in the blanks….

Published

March 23, 2026

using LinearAlgebra, Plots, LaTeXStrings, SparseArrays

We consider Poisson’s equation with homogeneous Dirichlet boundary conditions on the unit square [0,1]^2:

\begin{align} -\Delta u(x,y) &= f(x,y) % &&\text{for} \quad (x,y) \in(0,1) \times (0,1) \nonumber\\ % u(x,y) &= 0 &&\text{for} \quad x \in \{0,1\} \text{ or } y \in\{0,1\} % \tag{Poisson} \end{align}

where \Delta u := \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}. We approximate the solution of this equation on the uniform mesh \{(x_{i},x_j)\} where x_{j} = \frac{j}n for 0 \leq j \leq n. We define the grid spacing h := \frac{1}n.

Using the Taylor expansion, we have

\begin{align} \Delta u(x,y) = \frac{u(x - h,y) + u(x + h,y) + u(x, y-h) + u(x,y+h) - 4u(x,y)}{h^2} + O(h^2).\nonumber \end{align}

Therefore, if u_{i,j} := u(x_{i},x_j), then we can approximate the Laplacian in 2d with the following 5-point finite difference:

\begin{align} \Delta u(x_{ij}) % &\approx \frac{1}{h^2} \big( u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} - 4 u_{i,j} \big) \end{align}

In order to approximate \Delta (in 2d) as a matrix, we consider the flattening of u_{i,j} to obtain a vector in \mathbb R^{(n-1)^2} (again, we only consider the internal mesh points). We use vec() which does this flattening:

u = [1 4 7; 2 5 8; 3 6 9]
display(u)
vec(u)

3×3 Matrix{Int64}:
 1  4  7
 2  5  8
 3  6  9

9-element Vector{Int64}:
 1
 2
 3
 4
 5
 6
 7
 8
 9

(vec() uses a column-major format).

We let L \in \mathbb R^{(n-1)^2 \times (n-1)^2} be the matrix that maps the vectorisation of u on to the vectorisation of the finite difference \text{(1)}. We can then approximate the solution of \text{(Poisson)} by solving a linear system of the form -L u = f (in the same way you did in Assignment 4). Here is an implementation for reasonably small n:

n = 1000
g(x) = -x[1]^3
x = [[i/n, j/n] for j∈1:n-1, i∈1:n-1]
L = spdiagm(-1 => ones(n-2), 0 => -2ones(n-1), 1 => ones(n-2))
Id = Matrix(I, n-1, n-1)
L = kron(Id, L) + kron(L, Id)
u = reshape(( -L \ g.(vec(x)) )/n^2, n-1, n-1)

surface((1:n-1)/n, (1:n-1)/n, u, 
    title="Solution to -Δu = g (direct)", 
    xlabel="x", ylabel="y")

direct method (with dense matrices) took 5.6133964 seconds

It turns out that if A\in\mathbb R^{m\times m} has band width w (that is, A_{ij} = 0 if |i-j|>w) then the LU decomposition involves O(m w^2) operations. Our matrix L \in \mathbb R^{(n-1)^2 \times (n-1)^2} has band width w = n (why?) and therefore solving the linear system involves O( mw^2 ) = O(n^4) operations.

Instead, for large system sizes one uses iterative methods (such as CG). Compare CG with the above code for larger n. When does it become necessary to use iterative methods (on your machine)?

Other Links