复化simpson代码matlab求二重积分

复化 Simpson 公式可以用于数值求解二重积分，一般的实现方法是将二维区域分割成多个小矩形，然后对每个小矩形应用 Simpson 公式进行数值积分，最后将所有小矩形的积分值累加起来得到二重积分的近似值。具体实现如下：

function Q = composite_simpson2(f, a, b, c, d, m, n)
% f: 被积函数
% a: 积分区域 x 方向下限
% b: 积分区域 x 方向上限
% c: 积分区域 y 方向下限
% d: 积分区域 y 方向上限
% m: x 方向上的分割数，需要为偶数
% n: y 方向上的分割数，需要为偶数

hx = (b-a)/m; % x 方向梯形宽度
hy = (d-c)/n; % y 方向梯形宽度

% 构造分割点
x = linspace(a, b, m+1);
y = linspace(c, d, n+1);

% 计算每个小矩形上的积分值
Q = 0;
for i = 1:m
    for j = 1:n
        % 计算小矩形四个顶点的函数值
        f00 = f(x(i), y(j));
        f01 = f(x(i), y(j+1));
        f10 = f(x(i+1), y(j));
        f11 = f(x(i+1), y(j+1));
        % 计算小矩形的积分值
        Qij = hx*hy/9 * (f00 + f01 + f10 + f11 + 4*(f01+f10) + 16*f11);
        Q = Q + Qij;
    end
end
end

调用方式与上面的函数类似，例如要计算在矩形区域 $[0,\pi]\times[0,\pi/2]$ 上的二重积分 $\iint_{[0,\pi]\times[0,\pi/2]} \sin(xy) \, dxdy$，可以使用以下代码：

f = @(x,y) sin(x.*y); % 被积函数
a = 0; b = pi; % x 方向积分范围
c = 0; d = pi/2; % y 方向积分范围
m = 20; n = 20; % x 方向和 y 方向上的分割数，需要为偶数
Q = composite_simpson2(f, a, b, c, d, m, n); % 计算二重积分值

在这个例子中，使用了 $20\times20$ 个小矩形进行分割，得到的结果是 $Q \approx 0.6366$，与准确值 $0.6366$ 相当接近。需要注意的是，分割数 $m$ 和 $n$ 必须为偶数，否则会报错。