Numerical computation of double integral

[hunt_mat]

Hi,

I wrote a piece of MATLAB code to compute a double integral of the form:
<br /> \int_{a}^{b}\int_{c}^{d}f(x,y)dxdy<br />
I went about it using the trapezium rule, so what I did was apply the rule to the x variable first to obtain:
<br /> \int_{a}^{b}\int_{c}^{d}f(x,y)dxdy\approx\frac{\delta x}{2}\int_{c}^{d}2\sum_{i=1}^{N}f(x_{i},y)-f(a,y)-f(b,y)dy<br />
I then went on to apply the trapezium rule in the y variable to obtain:
<br /> \begin{array}{rcl}<br /> \int_{a}^{b}\int_{c}^{d}f(x,y)dxdy &amp; = &amp; \left(\sum_{i=1}^{N}\sum_{j=1}^{M}f(x_{i},y_{j})-\frac{1}{2}\sum_{i=1}^{N}f(x_{i},c)-\frac{1}{2}\sum_{i=1}^{N}f(x_{i},d)\right)\delta x\delta y \\<br /> &amp; - &amp; \frac{1}{2}\left(\sum_{j=1}^{M}f(a,y_{j})-\sum_{j=1}^{M}f(b,y_{j})\right)\delta x\delta y \\<br /> &amp; - &amp; \frac{1}{2}(f(a,c)+f(b,c)+f(a,d)+f(b,d))\delta x\delta y<br /> \end{array}<br />
I tried using this to compute the inverse Fourier transform of the 2d Gaussian and I got some rather bad results.
I include the results. Not what I expected.

 

[hunt_mat]

The picture of the results, On one exis it seems to work whereas on the other it doesn't.

Thoughts?

 

[I like Serena]

There seems to be missing a factor (1/2) before sum f(xi,d).
And the minus sign before f(b,yj) should probably be a plus sign.
Should I assume those are just typos?

Beyond that, your formula is symmetrical in x and y.
So if it yields an asymmetric result, you probably made a typing mistake.
Perhaps typing an i somewhere where there should be a j, or something like that.
Or perhaps one of the typos I already mentioned.

 

[hunt_mat]

There is the code I used:
function y=trap_2d(A,dx,dy)
N=length(A(:,1));
M=length(A(1,;

a=A(1,1)+A(1,M)+A(N,1)+A(N,M);
b=sum(A(1,+sum(A(N,;
c=sum(A(:,1))+sum(A(:,M));
u=zeros(1,N);
for i=1:N
u(i)=sum(A(i,;
end

d=sum(u);

y=(d-0.5*c-0.5*b-0.25*a)*dx*dy;

I first input the integrand as a matrix and then just do the sums in the equation as I said. Another thing that I am worried about is that as I am taking the4 inverse Fourier transform of the Gaussian, it shouldn't be coming out to that high a number.

 

[hunt_mat]

It turned out that it wasn't the algorithm for the 2d trapezium rule that was the problem, but the definition of the function itself. I am still off with the scaling though.

 

[I like Serena]

Sorry, but I've got no clue how to give a meaningful response to your last 2 replies.

 

[hunt_mat]

The integral algorithm wasn't the problem as it turns out, but another problem with the definition, I defined the 2D Gauss function as
<br /> \pi e^{(k^{2}+l^{2})/4}<br />
The algorithm I used to compute the inverse Fourier transform is:
const=1/(4*pi^2);
%Now compute the inverse tranform.

for i=1:n
for j=1:m
A=gauss_2d(K,L).*exp(sqrt(-1)*(K*x(i)+L*y(j)));
eta(i,j)=const*real(trap_2d(A,dx,dy));
end
end

I think that there is a problem here.

 

[I like Serena]

Your constant appears to be off.
I haven't checked your 2D Gauss function yet. Can I assume it is normalized?

I suspect your inverse transform should have a constant 1/2pi, as you can see here:
http://en.wikipedia.org/wiki/Fourier_transform#Two-dimensional_functions

 

[hunt_mat]

It appears to be fine, I am taking the 2d Fourier transform of:
<br /> f(x,y)=e^{-(x^{2}+y^{2})}<br />
This has a Fourier transform of:
<br /> \hat{f}(k,l)=\pi e^{-(k^{2}+l^{2})/4}<br />
So I don't think my constants off.