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.

  • $\begingroup$ Convert your problem to linear algebra, that's the most straightforward way to do it. $\endgroup$ Mar 12, 2023 at 20:16
  • $\begingroup$ I don't think there are any packages that will do everything for you, but this type of problem is very well-studied and there are a lot of resources out there. Given the size of your system, a matrix-free implementation of the conjugate-gradient method would probably work best, scipy.sparse.linalg.cg in Python. To do this, you need to write a function that applies the 2D Laplacian via the stencil above, then the RHS will be a vector containing your sources and boundary conditions. $\endgroup$
    – whpowell96
    Mar 13, 2023 at 4:01
  • $\begingroup$ More concretely, what @MaximUmansky is saying is that you can write the finite difference approximation of the Poisson equation in the form $A\varphi = h^2 I$ where $\varphi$ and $I$ are vectors. This is a standard system of linear equations for which there are many well-known algorithms (such as Conjugate Gradients). $\endgroup$ Mar 13, 2023 at 17:31
  • $\begingroup$ @WolfgangBangerth Thank you! Is it guaranteed that the solution will be a good approximation for the continuous solution? (since my lattice is discretized of course). More precisely - what can I say about this with regard to the accuracy the solution? $\endgroup$
    – Omer Paz
    Mar 13, 2023 at 19:27
  • $\begingroup$ @OmerPaz The scheme you outline is just one particular way of solving the linear system. Any other scheme is going to result in the same solution, with the same error guarantees, just possibly faster than the simple Jacobi iteration you are using. $\endgroup$ Mar 13, 2023 at 20:29

2 Answers 2


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)


(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)

enter image description here


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



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.