计算科学 - 谱域泊松求解器的数值误差 - 吾爱随笔录

在 $\mathbb{R}^n$ ，我想解一个泊松方程（给定 $f$ , 求解 $u$ ):

\nabla^{2} u = f

$\nabla^2 u = f$

假设 Neumann 边界条件（即 $\partial u = 0$ 在边界处）。

我在光谱域中解决了它。使用 3x3 拉普拉斯模板（因此 Neumann 边界 = 对称边界），DCT 和 IDCT 会很好。所以直接求解器是

u_{ours} = I D C T {\frac{D C T {f}}{D C T {\nabla^{2}}}}

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

这是一个测试 MATLAB 代码。

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');

问题：

分母 $\mathrm{DCT}\{\nabla^2\}$ 病态。下图显示在原点附近，值接近于零。

增加噪声水平（noise_level变量）会导致重建伪影。您可能会注意到伪影的频率较低：即几乎零除的结果。

图：左：noise_level = 1e-2；右：noise_level = 1e-1

图：误差图

问题：

为什么会发生这种情况？是因为光谱法的病态吗？
如果有，如何避免？这个问题还有哪些其他替代方法？
任何可能的正则化建议？