Оцените двойной интеграл, используя трапецеидальное правило (Matlab)

Question

У меня есть задание оценить двойной интеграл, используя трапецеидальное правило. Первая часть состояла в оценке двойного интеграла с помощью правила трапеции с пределами = х < = 2, 0 < = у < = 1

У меня есть рабочий скрипт для этого:

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) 

сейчас для второй части границы являются 0 < = х < = 2, 0 < = у < = ((р * х)/2)

Как получить у (строка 4 в коде) принимать соответствующие значения х из x-матрицы для создания y-mat RIX? Если я получу это для работы, мне не нужно ничего менять в коде, или я чего-то не хватает?

Answer 1

Интересная проблема. Представьте себе сетку, где координаты x охватывают от 0 <= x <= 2 и координаты y от 0 <= y <= pi. Это происходит потому, что когда x = 2, то y = pi*2/2 = pi.

Теперь представьте рисунок линии с положительным наклоном pi/2 от источника до (x,y) = (2,pi). Значения (x,y) вы хотите сосредоточиться на внизу справа из сетки.

Чтобы сделать это, просто создать meshgrid где x пролеты от [0,2] и y пролеты от [0,pi], а затем выберите в правом нижнем углу сетки. Предполагая сетку 100 х 100 точек:

%// 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)); 

Это цифра, мы получим:

регион в белом регионе мы после. ind содержит матрицу logical, которая позволит нам выбрать значения в пределах нашего meshgrid, которые нам нужно выбрать. Таким образом, ваш код теперь просто:

%// 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; 

Для out, я получаю:

>> out 

out = 

     0.156821355105871