目录
继续上一篇一维gauss积分的讨论，本文讨论二维gauss数值积分公式的使用，并给出数值实验。
一. 标准区间
按照一维情况类似方式，首先给出如下二重积分公式：
\\[ I = \\int_{-1}^1dx\\int_{-1}^1 f(x,y)dy=\\Sigma_{k=1}^m\\Sigma_{l=1}^mA_{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,\\quad 3\\times 3\\) 公式就完全足够了。所以不需要太多的点，徒增计算量。因此就可以得到gauss积分的坐标表达式为:
\\[ I = \\int_{-1}^1dx\\int_{-1}^1 f(x,y)dy=\\Sigma_{i=1}^{m^2}A_{i}f(x_i,y_i). \\]
二. 一般区间
对于一般区间，先考虑区间端点为常量情况（下一节介绍区间为变量情况），即
\\[ I = \\int_{a}^bdx\\int_{c}^d f(x,y)dy. \\]
其中 \\(a,b,c,d\\) 都是已知常数。与一维情况类似，只需要做变量变换，于是\\([s,t]\\in[-1,1]\\times [-1,1]\\)（通俗讲就是换元法）
\\[ 令 \\quad x =\\frac{b+a+(b-a)s}{2},\\qquad 则\\quad dx=\\frac{b-a}{2}ds,\\\\ \\quad y =\\frac{d+c+(d-c)t}{2},\\qquad 则\\quad dy=\\frac{d-c}{2}dt.\\\\ 记\\qquad jac = \\frac{(b-a)(d-c)}{4},\\qquad 则 \\quad dxdy=jac\\times dsdt \\]
于是二重积分就变成了
\\[ I = \\int_{a}^bdx\\int_{c}^d f(x,y)dy =jac\\times\\int_{-1}^1ds\\int_{-1}^1f(\\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}). \\]
其中 \\((s_i,t_i)\\) 即为上表中的 gauss 节点，对应的权重因子为 \\(A_i\\)。
三. 数值实验(没有编程的计算是不完整的)
使用matlab2018a 计算结果，并且与matlab自带函数 integral2 计算的结果进行比较给出误差。
算例如下:
\\[ 计算定积分 I = \\int_{a}^bdx\\int_{c}^d f(x,y)dy \\]
其中，\\(a=1.4,b=2,c=1,d=1.5,f(x,y)=ln(x+2*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 |
四. 总结和下节预告
- 从实验数据可以发现，二重gauss数值积分使用\\(2\\times 2\\) 4个点的情况下，结果已经很准确了，达到了1e-6的误差，所以在实际计算中，一般使用4点或者9点的gauss公式就可以满足要求了。
- 下一节介绍当积分区间在变系数下的二重gauss公式的计算方法
- 欢迎与我交流，数值分析，矩阵计算，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
