Fast computation of 2D-convolution using "kernel" in fourier space

Question

I want to numerically evaluate the following integral:

If the function f is given in discrete form on a grid, I want to be able to very efficiently compute the integral. The idea is to understand it as a multiplication in the Fourier space and use FFT:

if xi and eta are the Fourier variables. I tried to implement this in Matlab, but dont get exact results (probably because I dont quite understand the FFT).

I have the following minimal working example in 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 is the result of the integral, f is the input function and G is the fraction in the Fourier-space; the singularity at xi = 0 and eta = 0 of G is dealt with by taking the average value of the function in the surroundings of the point. I compare the final result to an analytical solution, which in this case exists. But it doesn't match perfectly (and doesn't get better, when I increase N? Why is that? And I also don't understand why the chess-matrix is necessary (I copied it from an example I found)?

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)

this is f(x,y) with a=.5

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

And f1 is the integrant, the actual function to integrate 2D

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

Setting x=.5;y=.5; as parameters.

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

resulting in

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) 

function integral2 needs some fields changed from default to return a result

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

Z =
   NaN

too strict

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

Now it works

I am submitting the above lines as answer while one of my machines is busy calculating the following

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)

The narrow gap is again MATLAB getting results but error out of specs.

Comment:

Double loop on integral2 is not time efficient for N>10 , at least for the machine I used to answer this question.

There are other ways to get a result for large N but since

understand that would be another question, on how to optimize an already working solution to the above question.