Solving 2D Poisson equation with nonhomogeneous boundary conditions (Dirichlet) and a source

Question

I am a physicist who is fairly new to numerical analysis, currently, I am trying to simulate a non-linear paraxial equation, and part of my calculation involves solving a 2D Poisson equation with Dirichlet boundary conditions and a source function.

In order to solve the 2D Poisson equation, I tried using the iterative relaxation method, meaning I repeated the process:

$$\varphi_\mathrm{new}(x) = \frac{1}{4}(\varphi(x+h,y) + \varphi(x-h,y) + \varphi(x,y+h) + \varphi(x,y-h)-h^2\cdot I(x,y))$$

With a uniformly spaced grid and $N=2048$ points on each axis. $I$ is a source function given (in my case a Gaussian distribution), and there are non-homogeneous boundary conditions. However, this method converges very slowly at this resolution.

Are there any other recommended methods that will be applicable to my case? It will be extremely helpful if they are either relatively simple or are available in Python libraries.

Answer 1

Here is a basic implementation of solving the problem $-\Delta u = f$ using matrix-free methods, which are prefereble for large PDE discretizations. The efficiency could be improved in the Laplacian evaluation and a preconditioner would save lots of time, but this runs in $\approx$ 7 minutes on my machine.

The method used to solve the linear system is conjugate-gradient, which only requires a function that compute matrix-vector products $Av$ for symmetric-positive definite matrices $A$. That is why I had to flip the sign to solve $-\Delta u = f$ to ensure positivity. The nonhomogeneous boundary conditions then have to be moved into the RHS vector. This can be derived from the finite-difference equations for the nodes near the edges of the domain, but it is easier to implement by reshaping u to be a matrix then just applying the loads on the edges. Also note that the plot doesn't display the boundary conditions exactly since it is only plotting the function values at the interior nodes.

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import LinearOperator, cg


def lap(u):
    # apply the discrete negative Laplacian with homogeneous Dirichlet boundary conditions

    N_dof = len(u)
    N = int(np.sqrt(N_dof))
    h2 = (N+1)*(N+1)
    u_mat = np.reshape(u, (N, N))

    Lu = 4*u_mat

    Lu[:-1, :] = Lu[:-1, :] - u_mat[1:, :]
    Lu[1:, :] = Lu[1:, :] - u_mat[:-1, :]
    Lu[:, 1:] = Lu[:, 1:] - u_mat[:, :-1]
    Lu[:, :-1] = Lu[:, :-1] - u_mat[:, 1:]

    Lu *= h2
    Lu = np.reshape(Lu, (N_dof))
    return Lu


def load(F, f1, f2, f3, f4):
    # construct load vector for system -\Delta u = f
    N = len(f1)
    b = F
    h2 = (N+1)*(N+1)

    # y=0
    b[0, :] += f1*h2

    # x=0
    b[:, 0] += f2*h2

    # y=L
    b[-1, :] += f3*h2

    # x=L
    b[:, -1] += f4*h2

    b = np.reshape(b, (N ** 2))
    return b


N = 2 ** 11
L = 1
h = L/(N+1)
N_dof = N ** 2

x = np.linspace(h, L-h, N)
xv = np.reshape(x, (1, N))
y = np.linspace(h, L-h, N)
yv = np.reshape(y, (N, 1))

# source
sig = 0.1
F = 100*np.exp((-(yv-0.5) ** 2)/sig ** 2) * \
    np.exp(-((xv-0.5) ** 2)/sig ** 2) / (2*np.pi*sig)

# y=0
f1 = 1*np.sin(x*np.pi)

# x=0
f2 = 2*np.sin(y*np.pi)

# y=L
f3 = 3*np.sin(x*np.pi)

# x=L
f4 = 4*np.sin(y*np.pi)


A = LinearOperator((N_dof, N_dof), matvec=lap)
b = load(F, f1, f2, f3, f4)

u, exitcode = cg(A, b)

print(exitcode)

(xg, yg) = np.meshgrid(x, y)
u_mat = np.reshape(u, (N, N))

ax = plt.axes(projection="3d")
ax.plot_surface(xg, yg, u_mat)
plt.show()

Answer 2

Here is a more complex implementation of the explicit SOR iteration with the compact 8-point stencil for higher accuracy. It is also Multithreaded in Cython if you wish to benchmark your CPU and post runtime here :)

Using the mask method of numbers you can put any arbitrary shape for dirichlet boundary conditions. You can even make them in MsPaint and import it as an array from numpy to asarray and use it as the mask.

Automatic calculation of the ideal omega or alpha value remains hard for the explicit method with custom b.c

MP_Gauss_Seidel_SOR