玩转 matlab 之二维 gauss 数值积分公式使用及 matlab 源代码(1)-常量区间

一. 标准区间

按照一维情况类似方式，首先给出如下二重积分公式：
$$I={\int }_{-1}^{1}dx{\int }_{-1}^{1}f(x,y)dy={\Sigma }_{k=1}^m{\Sigma }_{l=1}^m{A}_{kl}f(x_k,y_l).$$
这是标准的二维gauss积分 $m\times m$ (即 x 轴和 y 轴分别取 m 个点对应$m^2$个点)公式的表达式，其中，$x_{kl}$ 称作gauss积分点，$A_{kl}$称作$x_{kl}$ 点处的权重。$m=1,2,3$ 时候草图如下所示：

下面给出部分gauss积分点和权重对应列表

gauss点个数 $m^2(m\times m)$	gauss 点 $(x_i,y_i)$	权重 $A_i$	精度(2m-1)
1 ($1\times 1$)	$x_1=y_1=0$	$A_1$=4	1
4 ($2\times 2$)	$x_1=-1/\sqrt{3},y_1=-1/\sqrt{3}$ $x_2=1/\sqrt{3},y_2=-1/\sqrt{3}$ $x_3=1/\sqrt{3},y_3=1/\sqrt{3}$ $x_4=-1/\sqrt{3},y_4=1/\sqrt{3}$	$A_1=A_2=A_3=A_4=1$	3
9 ($3\times 3$) 令$gpt=\sqrt{3/5}$	$x_1=-gpt,y_1=-gpt$ $x_2=gpt,y_2=-gpt$ $x_3=gpt,y_3=gpt$ $x_4=-gpt,y_4=gpt$ $x_5=0,y_5=-gpt$ $x_6=gpt,y_6=0$ $x_7=0,y_7=gpt$ $x_8=-gpt,y_8=0$ $x_9=0,y_9=0$	$A_1=A_2=A_3=A_4=25/81$ $A_5=A_6=A_7=A_8=40/81$ $A_9=64/81$	5

一般二重积分近似值也就是使用 $2\times 2,3\times 3$ 公式就完全足够了。所以不需要太多的点，徒增计算量。因此就可以得到gauss积分的坐标表达式为:
$$I={\int }_{-1}^{1}dx{\int }_{-1}^{1}f(x,y)dy={\Sigma }_{i=1}^{m^2}{A}_{i}f(x_i,y_i).$$

二. 一般区间

对于一般区间，先考虑区间端点为常量情况（下一节介绍区间为变量情况），即
$$I={\int }_a^{b}dx{\int }_c^{d}f(x,y)dy.$$
其中 $a,b,c,d$ 都是已知常数。与一维情况类似，只需要做变量变换，于是$[s,t]\in [-1,1]\times [-1,1]$（通俗讲就是换元法）
$$\begin{gathered} 令x=\frac{b+a+(b-a)s}{2},则dx=\frac{b-a}{2}ds, \\ y=\frac{d+c+(d-c)t}{2},则dy=\frac{d-c}{2}dt. \\ 记jac=\frac{(b-a)(d-c)}{4},则dxdy=jac\times dsdt \end{gathered}$$
于是二重积分就变成了
$$\begin{gathered} I={\int }_a^{b}dx{\int }_c^{d}f(x,y)dy=jac\times {\int }_{-1}^{1}ds{\int }_{-1}^{1}f(\frac{b+a+(b-a)s}{2},\frac{d+c+(d-c)t}{2})dt. \\ =jac\times {\Sigma }_{i=1}^{m^2}{A}_{i}f(\frac{b+a+(b-a)s_i}{2},\frac{d+c+(d-c)t_i}{2}). \end{gathered}$$
其中 $(s_i,t_i)$ 即为上表中的 gauss 节点，对应的权重因子为 $A_i$。

三. 数值实验(没有编程的计算是不完整的)

使用matlab2018a 计算结果，并且与matlab自带函数 integral2 计算的结果进行比较给出误差。

算例如下:
$$计算定积分I={\int }_a^{b}dx{\int }_c^{d}f(x,y)dy$$
其中，$a=1.4,b=2,c=1,d=1.5,f(x,y)=ln(x+2\ast y),ln$是以e为底对数函数。使用matlab的integral2 函数计算结果为$I=0.429554527548275$.

自己编程计算结果如下：

高斯点数 $m^2(m\times m)$	积分值 $I_m$	误差norm($I_m-I$)
4($2\times 2$)	0.429556088022242	1.56E-06
9($3\times 3$)	0.429554531152490	3.60E-09

四. 总结和下节预告

1. 从实验数据可以发现，二重gauss数值积分使用$2\times 2$ 4个点的情况下，结果已经很准确了，达到了1e-6的误差，所以在实际计算中，一般使用4点或者9点的gauss公式就可以满足要求了。
2. 下一节介绍当积分区间在变系数下的二重gauss公式的计算方法
3. 欢迎与我交流，数值分析，矩阵计算，PDE数值解等，QQ群 315241287

五. matlab源代码

clc;clear;
%   compute int_a^b [int_c)^d    f(x,y)]
%   (x,y) \in [a,b] X [c,d]
%%  setup the integral interval and gauss point and weight
a = 1.4;    b = 2;
c = 1;      d =1.5;
fun=@(x,y)    log(x+2*y);
fprintf('***********************************************\n')
for gauss = 2:3 %    m points rule in 2 dimensional case
    if gauss == 2
        fprintf('*******  2X2 points gauss rule result *******')
        gpt=1/sqrt(3);
        s(1) = -gpt;  t(1) = -gpt;
        s(2) =  gpt;  t(2) = -gpt;
        s(3) =  gpt;  t(3) =  gpt;
        s(4) = -gpt;  t(4) =  gpt;
        wt = [1 1 1 1];        
    elseif gauss == 3
        gpt=sqrt(0.6);
        fprintf('*******   3X3 points gauss rule     *******')
        s(1) = -gpt; t(1) = -gpt; wt(1)=25/81;
        s(2) =  gpt; t(2) = -gpt; wt(2)=25/81;
        s(3) =  gpt; t(3) =  gpt; wt(3)=25/81;
        s(4) = -gpt; t(4) =  gpt; wt(4)=25/81;
        s(5) =  0.0; t(5) = -gpt; wt(5)=40/81;
        s(6) =  gpt; t(6) =  0.0; wt(6)=40/81;
        s(7) =  0.0; t(7) =  gpt; wt(7)=40/81;
        s(8) = -gpt; t(8) =  0.0; wt(8)=40/81;
        s(9) =  0.0; t(9) =  0.0; wt(9)=64/81;        
    end
    %%   区间变换到   [-1,1] X [-1,1]
    jac = (b-a)*(d-c)/4;
    x = (b+a+(b-a)*s)/2;
    y = (d+c+(d-c)*t)/2;
    f = fun(x,y);
    comp = wt(:) .* f(:) .* jac;%无论一个向量是行还是列，写成x(:)都会变成列向量
    
    format long
    comp = sum(comp)
    exact = integral2(fun,a,b,c,d);
    
    fprintf('the error is norm(comp-exact)=%10.6e\n\n',norm(comp-exact))
    
end
fprintf('******************************************\n')
fprintf('matlab  built-in function ''integral2''\n')
exact
format short