先前在做数值分析实验时，把高斯型求积公式和复合型、自适应步长的求积融合到了一起，但是后来发现题目没有这个要求。。现在就把这个思路分享一下。
实验目的：学会Gauss型求积公式,并应用该算法于实际问题.
实验内容：求定积分
∫
−
4
4
d
x
1
+
x
2
\int_{{\rm{ - }}4}^4 {\frac{{dx}}{{1 + {x^2}}}}
∫−441+x2dx
实验要求：
（1）把Gauss点的表格存入计算机,以Gauss-Legendre求积公式作为本实验的例子,要求程序可以根据不同的阶数,自动地用阶Gauss-Legendre求积公式计算上述定积分的近似值.体会Gauss型求积公式是具有尽可能高的代数精度的数值求积公式。
（2）可用MATLAB中的内部函数int求得此定积分的准确值与Gauss型求积公式求得的值进行比较。
首先是写出高斯型求积公式的代码，根据清华大学出版社出版的第五版数值分析上给出的定义：
简单来说就是求积公式具有2n+1次的代数精度，那就是高斯型求积公式。
接下来看高斯公式的代码实现（MATLAB）：
这里的输入参数分别为积分的上下限，n为节点个数，epsilon为求解精度。
function S=Gauss_Legendre(a,b,n,epsilon)
syms x
p=sym2poly(diff((x^2-1)^(n+1),n+1))/(2^(n+1)*factorial(n+1));
tk=roots(p); % 求积节点
% 计算求积系数
Ak=zeros(n+1,1);
for i=1:n+1
xkt=tk;
xkt(i)=[];
pn=poly(xkt);
fp=@(x)polyval(pn,x)/polyval(pn,tk(i)); % Lag基底函数
Ak(i)=quad(fp,-1,1,epsilon); % 计算求积系数
end
% 积分变量代换，将[a,b]变换到[-1,1]
xk=(b-a)/2*tk+(b+a)/2;
% 检验积分函数fun有效性
% fun1=fcnchk(fun1,'vectorize');
% 计算变量代换之后积分函数的值
fx=fun(xk)*(b-a)/2;
% 计算积分值
S=sum(Ak.*fx);
接下来的复合型高斯求积公式思路很简单，就是将所求区间分割，将多个区间套入高斯型求积公式进行求解，最后将结果相加。
我的思路是先将区间二分，再进行高斯型求积公式。很自然的，我们可以将二分出来的两个区间看做父区间生成的子区间，子区间再生成子区间，即可以当做二叉树进行处理。那么有些枝生成后，对应的区间函数变化幅度小，就可以不再进一步生成子区间，我们就可以把枝减去。那么进行如此迭代后即可产生我们任意指定精度的求积公式。
代码如下：
以下是主函数
%main.m
clc;clear;
a=-4;b=4;n=5;epsilon=1e-10;jd=1e-8;
syms x;global sum
p=1;q=1;sum=0;
Binary_tree(1,1,n,jd);
sum
S1=double(int(fun(x),a,b))
diff=abs(sum-S1)
以下是所求的积分函数
function y=fun(x)
y=1./(1+x.^2);
% y=sqrt(x);
return
这是生成二叉树的迭代函数
function Binary_tree(a,b,n,jd) %a,b为数值位于二叉树的位置
global S sum;
epsilon=1e-8;
h=8/(2^(a));
if 1
if a==1&&b==1
S(1,1).root=Gauss_Legendre(-4,4,n,epsilon);
end
if 1
S(a+1,2*b-1).root=Gauss_Legendre(-4+h*(2*b-2),-4+h*(2*b-1),n,epsilon);
S(a+1,2*b).root=Gauss_Legendre(-4+h*(2*b-1),-4+h*(2*b),n,epsilon);
S(a,b).L=S(a+1,2*b-1).root;S(a,b).R=S(a+1,2*b).root;
diff=abs(S(a,b).root-S(a,b).L-S(a,b).R)*2^(a-1);
if(diff<jd)
S(a,b).d=1;sum=sum+S(a,b).R+S(a,b).L;
else
S(a,b).d=0;
Binary_tree(a+1,b*2-1,n,jd);
end
if(b+1<=2^(a-1)&&mod(b,2)==1)
Binary_tree(a,b+1,n,jd);
end
end
end
运行结果如下：
效果还是很不错的。