在傅立叶空间中使用“核”快速计算 2D 卷积

Question

我想对以下积分进行数值计算：

如果函数 f 在网格上以离散形式给出，我希望能够非常有效地计算积分。这个想法是将其理解为傅立叶空间中的乘法并使用 FFT：

如果 xi 和 eta 是傅里叶变量。我尝试在 Matlab 中实现这一点，但没有得到确切的结果（可能是因为我不太了解 FFT）。

我在 Matlab 中有以下最小工作示例：

N = 101;
x = linspace(-1, 1, N); hx = x(2)-x(1);
[X, Y] = meshgrid(x, x);

a = 1/2;
f = 2/pi*a*real(sqrt(1-(X.^2+Y.^2)/a^2));
f_pad = zeros(2*N, 2*N); f_pad(1:N, 1:N) = f;

k = linspace(0, pi/hx, N+1); k = [k k(end-1:-1:2)]; hk = k(3)-k(2);
[KX, KY] = meshgrid(k, k);
K = sqrt(KX.^2 + KY.^2);

chess = mod((1:2*N)+(1:2*N)',2); chess = -chess; chess(chess==0) = 1;

G = 1./K;
G(1,1) = 4*integral2(@(x,y) 1./(sqrt(x.^2+y.^2)),0,0.5*hk,0,0.5*hk)/hk^2;

R = ifft2(G.*2.*chess.*fft2(f_pad), 'symmetric'); 
R = R(N+1:end, N+1:end);

Rana = (a^2-x.^2/2) .* (abs(x)<=a) + ...
       real(a^2/pi*((2-x.^2/a^2).*asin(a./abs(x))+sqrt(x.^2-a^2)/a)) .* (abs(x)>a);

plot(Rana);
hold on;
plot(R((N-1)/2+1,:));

R 是积分结果，f 是输入函数，G 是傅立叶空间中的分数； G 的 xi = 0 和 eta = 0 处的奇点通过取该点周围函数的平均值来处理。 我将最终结果与解析解进行比较，在本例中解析解是存在的。但它并不完美匹配（并且当我增加 N 时并没有变得更好？为什么会这样？而且我也不明白为什么国际象棋矩阵是必要的（我从我找到的示例中复制了它）？

Answer 1

N = 101;
x = linspace(-1, 1, N); hx = x(2)-x(1);
y = linspace(-1, 1, N); hy = y(2)-y(1);
[X, Y] = meshgrid(x, x);

a = 1/2;
F1 = 2/pi*a*real(sqrt(1-(X.^2+Y.^2)/a^2));   % f(x,y)

这是

f(x,y)

和

a=.5

figure(1)
hs=surf(X,Y,F1)
hs.EdgeColor='none';

并且

f1

是积分，实际积分 2D 的函数

f1 = @(x_,y_) 2/pi*a*real(sqrt(1-(x_.^2+y_.^2)/a^2))

将

x=.5;y=.5;

设置为参数。

x0=.5;y0=.5;
F2=@(x1,y1) f1(x1,y1)./((x0-x1).^2+(y0-y1).^2).^.5;  

Z=integral2(F2,x(1),x(end),y(1),y(end));

导致

Warning: Reached the maximum number of function evaluations
(10000). The result fails the global error test. 
> In integral2Calc>integral2t (line 129)
In integral2Calc (line 9)
In integral2 (line 105)
In double_integral_example (line 23) 
Warning: Reached the limit on the maximum number of intervals
in use. Approximate bound on error is   1.6e-15. The integral
may not exist, or it may be difficult to approximate
numerically to the requested accuracy. 
> In integralCalc/iterateScalarValued (line 372)
In integralCalc/vadapt (line 132)
In integralCalc (line 75)
In integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)
In integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) (line 17)
In integralCalc/iterateScalarValued (line 323)
In integralCalc/vadapt (line 132)
In integralCalc (line 75)
In integral2Calc>integral2i (line 20)
In integral2Calc (line 7)
In integral2 (line 105)
In double_integral_example (line 26) 
Warning: The integration was unsuccessful. 
> In integral2 (line 108)
In double_integral_example (line 26) 

函数

integral2

需要更改一些默认字段才能返回结果

Z=integral2(F2,x(1),x(end),y(1),y(end),'Method','iterated','AbsTol',0,'RelTol',1e-10)

Z =
   NaN

太严格了

Z=integral2(F2,x(1),x(end),y(1),y(end),'Method','iterated','AbsTol',1e-5,'RelTol',1e-8)
    
    Z =
       0.250008992645364

现在可以了

我正在提交以上几行作为答案，而我的一台机器正忙于计算以下内容

F3=zeros(N,N);

for k1=1:1:numel(x)
    for k2=1:1:numel(y)
        
        x0=x(k1);y0=y(k2);
        F2=@(x1,y1) f1(x1,y1)./((x0-x1).^2+(y0-y1).^2).^.5; 
        
        F3(k1,k2)=integral2(F2,x(1),x(end),y(1),y(end),'Method','iterated','AbsTol',1e-5,'RelTol',1e-8);
        
    end
end

figure(2)
hs3=surf(X,Y,F3)

我计算了 2.7 小时，完成了约 40 分钟。一旦

figure(2)

图表可用，我就会将其包含在此处。

评论：

i

ntegral2

上的双循环对于

N>10

来说时间效率不高，至少对于我用来回答这个问题的机器来说是这样。

还有其他方法可以在不浪费 2.7 小时的机器的情况下获得所有分数，但我知道这将是另一个问题，即如何优化上述问题的已经有效的解决方案。