Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Numerical computation of double integral

  1. Apr 22, 2012 #1

    hunt_mat

    User Avatar
    Homework Helper

    Hi,

    I wrote a piece of MATLAB code to compute a double integral of the form:
    [tex]
    \int_{a}^{b}\int_{c}^{d}f(x,y)dxdy
    [/tex]
    I went about it using the trapezium rule, so what I did was apply the rule to the [itex]x[/itex] variable first to obtain:
    [tex]
    \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
    [/tex]
    I then went on to apply the trapezium rule in the [itex]y[/itex] variable to obtain:
    [tex]
    \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}
    [/tex]
    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. jcsd
  3. Apr 22, 2012 #2

    hunt_mat

    User Avatar
    Homework Helper

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

    Thoughts?
     

    Attached Files:

  4. Apr 23, 2012 #3

    I like Serena

    User Avatar
    Homework Helper

    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.
     
  5. Apr 23, 2012 #4

    hunt_mat

    User Avatar
    Homework Helper

    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.
     
  6. Apr 23, 2012 #5

    hunt_mat

    User Avatar
    Homework Helper

    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:

  7. Apr 23, 2012 #6

    I like Serena

    User Avatar
    Homework Helper

    Sorry, but I've got no clue how to give a meaningful response to your last 2 replies.
     
  8. Apr 23, 2012 #7

    hunt_mat

    User Avatar
    Homework Helper

    The integral algorithm wasn't the problem as it turns out, but another problem with the definition, I defined the 2D Gauss function as
    [tex]
    \pi e^{(k^{2}+l^{2})/4}
    [/tex]
    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.
     
  9. Apr 23, 2012 #8

    I like Serena

    User Avatar
    Homework Helper

  10. Apr 23, 2012 #9

    hunt_mat

    User Avatar
    Homework Helper

    It appears to be fine, I am taking the 2d Fourier transform of:
    [tex]
    f(x,y)=e^{-(x^{2}+y^{2})}
    [/tex]
    This has a Fourier transform of:
    [tex]
    \hat{f}(k,l)=\pi e^{-(k^{2}+l^{2})/4}
    [/tex]
    So I don't think my constants off.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Numerical computation of double integral
Loading...