使用梯形法则计算二重积分 (Matlab)
Evaluate double integral using trapezoidal rule (Matlab)
我有一个作业要使用梯形法则计算二重积分。第一部分是使用带极限的梯形法则计算二重积分
0 <= x <= 2, 0 <= y <= 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 <= x <= 2, 0 <= y <= ((pi*x)/2)
如何让 y(代码中的第 4 行)从 x 矩阵中获取相应的 x 值来创建 y 矩阵?如果我让它工作,我就不必更改代码中的任何其他内容,或者我是否遗漏了什么?
有趣的问题。想象一个网格,其中 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
矩阵,它允许我们 select 我们 meshgrid
中的哪些值我们需要 select。因此,您的代码现在就是这样:
%// 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
我有一个作业要使用梯形法则计算二重积分。第一部分是使用带极限的梯形法则计算二重积分 0 <= x <= 2, 0 <= y <= 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 <= x <= 2, 0 <= y <= ((pi*x)/2)
如何让 y(代码中的第 4 行)从 x 矩阵中获取相应的 x 值来创建 y 矩阵?如果我让它工作,我就不必更改代码中的任何其他内容,或者我是否遗漏了什么?
有趣的问题。想象一个网格,其中 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
矩阵,它允许我们 select 我们 meshgrid
中的哪些值我们需要 select。因此,您的代码现在就是这样:
%// 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