Poisson Solver Example

num::pde::poisson2d solves the Dirichlet problem

\[ -\Delta u = f,\qquad u|_{\partial\Omega}=0 \]

on an \(N\times N\) interior grid by DST-I diagonalization.

Solve on a Square Grid

constexpr int N = 63;
constexpr double h = 1.0 / (N + 1);
 
num::Matrix f(N, N, 0.0);
for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
        double x = (i + 1) * h;
        double y = (j + 1) * h;
        f(i, j) = 2.0 * std::numbers::pi * std::numbers::pi
                * std::sin(std::numbers::pi * x)
                * std::sin(std::numbers::pi * y);
    }
}
 
num::Matrix u = num::pde::poisson2d(f, N);

The exact solution in this example is

\[ u(x,y)=\sin(\pi x)\sin(\pi y). \]

Finite-Difference Reference

num::Matrix u_fd = num::pde::poisson2d_fd(f, N);

Use poisson2d_fd as a direct finite-difference reference for small grids.

Size Constraint

The DST implementation uses FFTs of length \(2(N+1)\). The current built-in radix-2 path requires \(N+1\) to be a power of two, for example \(N=7,15,31,63,127\).