Search code examples
matlabnumerical-integration

How to use integration limits correctly (Matlab)


I've made a script that calculates a given integral with the limits:

0 <= x <= 2 and 0 <= y <= 1 

But now I want to change the limits to:

0 <= x <= 2 and 0 <= y <= sin((pi*x)/2) 

Function:

function f = inte(x,y)

dz = 10;

f = exp(-dz*((x-1.25).^2+y.^2)).*cos(y.*(x-1.25));

end

This is my script for the earlier limits:

L = 100; M = L/2;

hx = (2)/L; hy = (1)/M; 

x=[0:L]*hx;
y=[0:M]*hy;

Fx=[];

for i = 1:L+1

    Fy=[];

    for j = 1:M+1

        f = inte(x(i),y(j));
        Fy = [Fy f];

    end

    ycor = hy*(sum(Fy) - Fy(1)/2 - Fy(end)/2);
    Fx = [Fx ycor];

end

ans = hx*(sum(Fx) - Fx(1)/2 - Fx(end)/2);

disp (ans)

I can't seem to get the right answer when I try to change the code. The correct answer should be 0.1560949...

L is amount of steps in x direction, M in y direction. hx and hy are step lengths. This is really bugging me. And no I can only use the commands integral2 or traps as reference.

Thanks in advance!


Solution

  • In your present code, the lines

    hy = (1)/M; 
    y=[0:M]*hy;
    

    refer to the y-variable. When the limits for y depend on x, these lines cannot stay outside of the loop over x: they should be moved in and use the value x(i). Like this:

    for i = 1:L+1   % as in your code
    
       hy = (sin(pi*x(i)/2))/M;     
       y = [0:M]*hy;
    
       Fy=[];  % this and the rest as in your code 
    

    I get output 0.1561, as you wanted.