# Numerical computation of double integral

1. Apr 22, 2012

### hunt_mat

Hi,

I wrote a piece of MATLAB code to compute a double integral of the form:
$$\int_{a}^{b}\int_{c}^{d}f(x,y)dxdy$$
I went about it using the trapezium rule, so what I did was apply the rule to the $x$ variable first to obtain:
$$\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$$
I then went on to apply the trapezium rule in the $y$ variable to obtain:
$$\begin{array}{rcl} \int_{a}^{b}\int_{c}^{d}f(x,y)dxdy & = & \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 \\ & - & \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 \\ & - & \frac{1}{2}(f(a,c)+f(b,c)+f(a,d)+f(b,d))\delta x\delta y \end{array}$$
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.

Last edited: Apr 23, 2012
2. Apr 22, 2012

### hunt_mat

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

Thoughts?

#### Attached Files:

• ###### gauss_2d.png
File size:
33.4 KB
Views:
159
3. Apr 23, 2012

### 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.

4. Apr 23, 2012

### 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.

5. Apr 23, 2012

### 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.

#### Attached Files:

• ###### gauss_2d.jpg
File size:
9.7 KB
Views:
141
6. Apr 23, 2012

### I like Serena

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

7. Apr 23, 2012

### 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
$$\pi e^{(k^{2}+l^{2})/4}$$
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.

8. Apr 23, 2012

### I like Serena

9. Apr 23, 2012

### hunt_mat

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