使用梯形规则的 MATLAB 双积分

Question

我必须用梯形复合法则来积分：

%f.m
function [ab]=f(x,y)
ab = sqrt((exp(sin(x.*pi - y.*pi)))+(exp(cos(x.*y.*pi))))

还有我的梯形函数：

%trap.m
function [T] = trap(a,b,N)
h=(b-a)./N;
x=a+[0:N]*h;
y=f(x);
T=(y(1)+y(N+1)+2*sum(y(2:N)))*h/2

我对如何进行双重积分感到困惑，但我想它会得到类似的东西：

%main.m
clear;
Tr = zeros(6);
a=0;
b=3 ; 
N = 2.^([1:10]);

for i=1:2:10,
    [Tr(i)]=trap(a, b, N(i), %first integration); %second integration...
end

我的陷阱功能需要第一次集成？

Answer 1

在一维梯形规则中，函数值乘以 h/2, h, h, ..., h, h/2，其中 h 是步长。这些是该积分规则的权重。当积分规则应用于二维时，权重会相乘。一般来说，您在二维上会有不同的步长，例如 h1 和 h2，权重将为 h1*h2 乘以以下数字：

1/4 1/2 1/2 1/2 1/4
1/2  1   1   1  1/2
1/2  1   1   1  1/2
1/4 1/2 1/2 1/2 1/4

实现此目的的一种快速方法是将所有值相加，减去每边的 1/2 值（这会减去角点两次），然后添加 1/4 角点值。不需要

for

循环：您可以使用

meshgrid

设置评估网格，使用

sum

进行双重求和。示例：

f = @(x,y) sqrt((exp(sin(x.*pi - y.*pi)))+(exp(cos(x.*y.*pi))));
a1 = 0; b1 = 3; N1 = 20;
h1 = (b1-a1)/N1;
a2 = -1; b2 = 4; N2 = 30;
h2 = (b2-a2)/N2;
[X,Y] = meshgrid(a1:h1:b1, a2:h2:b2);
Z = f(X,Y); 
T = h1*h2*(sum(Z(:)) - 0.5*(sum(Z(1,:)+Z(end,:)) + sum(Z(:,1)+Z(:,end))) + 0.25*(Z(1,1)+Z(1,end)+Z(end,1)+Z(end,end)))

最后一行与我上面所说的完全一致：添加所有值，减去 1/2 边，添加 1/4 角。