CODEWITHSUNDEEP
matlabnumerical-integration
stx1020
stx1020

Reputation: 45

Evaluate double integral using trapezoidal rule (Matlab)

I have an assignment to evaluate a double integral using trapezoidal rule. The first part was to evaluate a double integral using trapezoidal rule with limits 0 <= x <= 2, 0 <= y <= 1

I have a working script for that:

N = 100;
xh= 1.25;
x = linspace(0,2,N);
y = linspace(0,1,0.5*N);
dx = diff(x(1:2));
dy = diff(y(1:2));
[x,y] = meshgrid(x,y);
funk = exp(-10.*((x-xh).^2+y.^2)).*cos(y.*(x-xh));
funk(2:end-1,:) = funk(2:end-1,:)*2;
funk(:,2:end-1) = funk(:,2:end-1)*2;
out = sum(funk(:))*dx*dy/4;

disp(out)

Now for the second part the limits are 0 <= x <= 2, 0 <= y <= ((pi*x)/2)

How do I get y (line 4 in code) to take corresponding x values from the x-matrix to create the y-matrix? If I get that to work I shouldn't have to change anything else in the code right, or am I missing something?

Upvotes: 2

Views: 1837

Answers (1)

rayryeng
rayryeng

Reputation: 104493

Interesting problem. Imagine a grid where the x coordinates span from 0 <= x <= 2 and the y coordinates span from 0 <= y <= pi. This is because when x = 2, then y = pi*2/2 = pi.

Now imagine drawing a line with a positive slope of pi/2 from the origin up to (x,y) = (2,pi). The values of (x,y) you want to focus on are the bottom right of the grid.

To do that, simply create a meshgrid where x spans from [0,2] and y spans from [0,pi], then choose the bottom right of the grid. Assuming a grid of 100 x 100 points:

%// Generate grid of points
N = 100;
xx = linspace(0,2,N);
yy = linspace(0,pi,N);
[x,y] = meshgrid(xx,yy);

%// Obtain valid region
ind = y <= (pi/2)*x;

%// Show valid region in black and white
imagesc(ind);
axis xy;
colormap gray;
set(gca,'XTick',10:10:100);
set(gca,'YTick',10:10:100);
set(gca,'XTickLabel',xx(10:10:end));
set(gca,'YTickLabel',yy(10:10:end));

This is the figure we get:

enter image description here

The region in white is the region we are after. ind contains a logical matrix that will allow us to select which values within our meshgrid we need to select. As such, your code now is simply this:

%// Generate grid of points
N = 100;
xh = 1.25;
xx = linspace(0,2,N);
yy = linspace(0,pi,N);
[x,y] = meshgrid(xx,yy);

%// Obtain valid region
ind = y <= (pi/2)*x;

%// Perform calculations with normal grid
dx = diff(xx(1:2));
dy = diff(yy(1:2));
funk = exp(-10.*((x-xh).^2+y.^2)).*cos(y.*(x-xh));
funk(2:end-1,:) = funk(2:end-1,:)*2;
funk(:,2:end-1) = funk(:,2:end-1)*2;

%// Select out valid region coordinates
funk = funk(ind);

%// Now sum
out = sum(funk(:))*dx*dy/4;

For out, I get:

>> out

out =

         0.156821355105871

Upvotes: 2

Related Questions