MATLAB: 2D Integration with infinite limits (Fourier Transform)

2d integralanalytical fourier transforminfinite limitsquadgk

Hi All,
I have to take a Fourier transform of a rather complicated 2D function, however my impression is that creating large arrays and doing it with fft would not be a good idea (I need to get as close to infinity with both both kx and ky –see below– as possible).
I have used Mike Hosea's mydblquad function but it is throwing a lot of tolerance warnings and taking forever to compute. Here's what I'm doing right now:
%%Constants (not all unity, but easier for demonstration purposes)
Ly = 1;
Wx = 1;
x = 1;
y = 1;
%%Variables which I eventually want to loop over (nested for-loop)
c1 = 1;
c2 = 1;
%%Function:
L = @(kx,ky) sqrt(kx.^2+ky.^2);
u = @(kx,ky) (sqrt(L(kx,ky).^2+1i*c1)./c2);
RL = @(kx,ky) (L(kx,ky)- u(kx,ky))./(L(kx,ky)+u(kx,ky));
fun = @(kx,ky) (RL(kx,ky)./L(kx,ky)).*sinc(kx.*Wx).*sinc(ky.*Ly).*exp(1i*(x.*kx+y.*ky));
Result = mydblquad(fun,-inf,inf,-inf,inf);
As you can see, minus some factors of 2*pi floating around, this just the Fourier transform of
(RL(kx,ky)./L(kx,ky)).*sinc(kx.*Wx).*sinc(ky.*Ly)
Does anyone have suggestions on how to either speed this up, get rid of the tolerance warnings or compute it in a completely different way?
Thanks in advance!
-Bradford

Best Answer

You're using L as both a function and as a variable there. Must have happened when you were trying to simplify the presentation. Anyway, I just changed sinc(ky.*L) to sinc(ky.*W) in the defintion of fun and got the example to run.
You have a singularity at the origin. Assuming it is integrable, you should at the very least split this up into 4 integrals over the respective quadrants (singularities need to be on the boundary for QUADGK). You can also experiment by splitting it up into more pieces. The finite (proper) sub-integrals might be possible for QUAD2D to do, and it's usually faster than iterated QUADGK. Again, any singularities should be positioned on the boundaries of sub-integrations.
Also, you can set 'AbsTol' and 'RelTol' for QUADGK. If you can live with only a few digits of accuracy, increasing 'AbsTol' and 'RelTol' should make things faster in general.
Related Question