Numerical error of a spectral-domain Poisson solver

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?
• Your plots would be more useful if you included a color bar to show the size of the “error” – Brian Borchers Jun 4 '18 at 2:50

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.
• Thanks for the insight. Sorry for my blunt. To articulate, if I understand it correctly, by saying "solving the Poisson equation for two different RHS", you mean my $f$ is noisy so the solution must be different somehow. Yes, and I showed the artifact as some low frequency; and I suspect the ill-conditioning of DCT basis is the reason to blame. So it leads to my question: can I use other approaches to avoid such low-frequency artifacts? – WDC Jun 4 '18 at 22:12
• Yes, I say that if the $f$ is noisy that noise has to appear in the solution. To test your scheme you can do the following. First, create exact solution $u_{ext}$=whatever_you_want. Second, compute RHS $f=\nabla^2 u_{ext}$. Third, use your method to solve the Poisson problem $\nabla^2 u = f$ for the numerical solution $u$. Four, compute the difference $|u - u_{ext}|$, this difference is the error. Other discretizations are not going to be better, all discretizations are similar for low frequencies ($\sim k^2$). – floren Jun 5 '18 at 20:05
• Setting noise_level = 0 meets your suggestion here. And I got the difference as $10^{-13}$, the eps. I think my Poisson solver is theoretically fine. So the question is, when there is noise, how can I handle them in a more appropriate way. – WDC Jun 5 '18 at 22:20
• If you believe that your $f$ is contaminated with noise and would otherwise be smooth, then you can filter $f$ to smooth it before applying your solution algorithm. You'll have to decide how much smoothing would be appropriate. – Brian Borchers Aug 4 '18 at 1:30