有趣的问题。想象一个网格,其中x 坐标跨越0 <= x <= 2,y 坐标跨越0 <= y <= pi。这是因为当x = 2,然后y = pi*2/2 = pi。
现在想象从原点到(x,y) = (2,pi) 绘制一条斜率为pi/2 的直线。您要关注的(x,y) 的值是网格的右下角。
为此,只需创建一个meshgrid,其中x 跨越[0,2],y 跨越[0,pi],然后选择网格的右下角。假设一个 100 x 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