Numerical error of a spectral-domain Poisson solver

Question

In $\mathbb{R}^n$, I would like to solve a Poisson equation (given $f$, solve for $u$):

$$\nabla^2 u = f$$

assuming Neumann boundary condition (i.e. $\partial u = 0$ at boundaries).

I solved it in spectral domain. Using the 3x3 Laplacian stencil (so Neumann boundary = symmetric boundary), DCT & IDCT will be good. So the direct solver is

$$ u_{\text{ours}} = \mathcal{IDCT}\{ \frac{\mathcal{DCT}\{ f \}}{\mathcal{DCT\{\nabla^2\} }} \}$$

Here is a test MATLAB code.

clc;close all;
% Solve for u from f: nabla^2 u = f 
% (under Neumann boundary condition)

% for reproduction
rng(0);

% generate u
Udim = [256 256];
coeff = 0.2;
[~, X] = meshgrid(coeff*(1:Udim(2)),coeff*(1:Udim(1)));
u = sin(X);

% generate f
lambda = 2e1;
f = lambda*imfilter(u, fspecial('laplacian',0), 'symmetric');

% add Gaussian noise to f
noise_level = 1e-2;
f = f + noise_level*randn(size(f));


%%% Now solve for u from noisy f
% generate the DCT of f
numerator = dct2(f);

% generate the DCT denominator
[H,W] = size(f);
[x,y] = meshgrid(0:W-1,0:H-1);
denominator = 2*cos(pi*x/W) + 2*cos(pi*y/H) - 4;
denominator(1) = 1;     % set DC to be 1; does not matter

% inversion in DCT domain
numerator = idct2( numerator./( lambda*denominator) );

% zero-mean normalization
u = u - min(u(:));
u_ours = numerator - min(numerator(:));

% show results
figure;
subplot(131);imshow(u, [],'i','f');colormap(jet);title('u (real)');
subplot(132);imshow(u_ours, [],'i','f');colormap(jet);title('u (ours)');
subplot(133);imshow(u - u_ours, [],'i','f');
colormap(jet);   title('error');

Issues:

• The denominator $\mathrm{DCT}\{\nabla^2\}$ is ill-condition. Below figure shows that near the origin, the values are close to zero.

• Increasing noise level (the noise_level variable) leads to reconstruction artifacts. You may notice the artifacts are of low frequency: i.e. the consequence of almost-zero-division.

Figure: Left: noise_level = 1e-2; Right: noise_level = 1e-1

Figure: Error plots

Questions:

• Why this happens? Is it because of the ill conditioning of the spectral method?
• If so, how to avoid it? What are other alternative methods for this issue?
• Any possible regularization suggestions?

Answer 1

The Laplacian is $\sim k^2$ for low wavenumbers; it is a general property and it does not depend of your discretization. You are solving the Poisson equation for two different RHS so the solutions are going to be different. You cannot call that difference "error", the solutions have to be different.