2D Integration with infinite limits (Fourier Transform)

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

Answer 1

[Mike Hosea]

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.

[Bradford]

You're right, putting the singularities at the edge of integration and splitting up the origin into quadrants using QUAD2D did a good job of dealing with the singularity. Still having some warnings and larger errors with QUADGK as the limits of integration increase, but overall much better.

Thanks!