二维区间上高斯数值积分的MATLAB实现

在之前发布的一维区间高斯数值积分实现中，我使用了通过syms定义的符号函数。经过大量实践发现，在一直调用syms函数的过程中，会让整个程序增加很多时间，效率十分低下。

因此，在二维区间上，我将不再使用符号函数，转而使用句柄函数来进行实现高斯数值积分。

1. 首先，我先介绍一下何谓符号函数和句柄函数。

例如问题中的函数，假如在Matlab中用符号函数去定义，那就是：

syms x y t;
f = (1-2*x)*(y-1)*cos(t);

用句柄函数去定义，是：

f = @(x,y,t) (1-2*x)*(y-1)*cos(t);

使用后者，会直接省去Matlab中自带的syms函数的调用时间。但事情都是有两面性的，句柄函数能节省时间，但它在代码的编写阶段会提升复杂程度。

例如，我们对函数进行换元，，，符号函数的代码如下：

syms x y s t;
g = (-2/3)*x*y^3+2-pi*sin(pi*x);
x_hat = s*cos(t);
y_hat = s*sin(t);
g_hat = subs(g,{x,y},{x_hat,y_hat});

句柄函数的代码如下：

g = @(x,y) (-2/3)*x*y^3+2-pi*sin(pi*x);
x_hat = @(s,t) s*cos(t);
y_hat = @(s,t) s*sin(t);
g_hat = @(s,t) g(x_hat(s,t),y_hat(s,t));

可以看到，句柄函数的换元，我们还需要在程序中考虑换元后函数的自变量是什么，增加到代码当中；符号函数则不需要，直接换元就可以了。

2. 现在我们进入正题，用句柄函数来实现二维区间上的高斯数值积分。

首先我们给出高斯点和高斯权重。对标准区间的三级形，它的高斯点及权重为：

高斯点	x	y	权重
1	1/3	1/3	-27/96
2	3/5	1/5	25/96
3	1/5	1/5	25/96
4	1/5	3/5	25/96

代码如下：

function [r]=gauss_int_2D(f)
T = [1/3,3/5,1/5,1/5;1/3,1/5,1/5,3/5];
Ak = [-27/96,25/96,25/96,25/96];
r = 0;
for i=1:4
    Xk = T(1,i);
    Yk = T(2,i);
    f1 = f(Xk,Yk);
    r = r+Ak(i)*f1;
end
end

3. 另外，利用句柄函数，我们还可以将这个程序进行向量化处理，去掉for循环。

首先，我们需要先考虑乘法。两个向量对应元素相乘（除），需要用点乘（除）。

如上述代码中定义的向量T，第一行乘以第二行，是

T(1,:).*T(2,:);

因此，定义句柄函数时，假如涉及到自变量的乘法和除法，我们都需要用点乘和点除。

4. 最后，我们给出对函数进行在标准区间上的高斯数值积分实现程序。

% 主程序
h = @(x,y) -(1-2*x).*(y-1)+((1-x).*(y-1)-x.*(y-1));
r = gauss_int_2D(h)

%向量化处理后的高斯数值积分子程序
function [r]=gauss_int_2D(f)
T = [1/3,3/5,1/5,1/5;1/3,1/5,1/5,3/5];
Ak = [-27/96,25/96,25/96,25/96];
r = sum(f(T(1,:),T(2,:)).*Ak)
end

本次的分享就这些了~