复合型自适应步长的Gauss型求积

先前在做数值分析实验时，把高斯型求积公式和复合型、自适应步长的求积融合到了一起，但是后来发现题目没有这个要求。。现在就把这个思路分享一下。

上题目：

实验目的：学会Gauss型求积公式,并应用该算法于实际问题.
实验内容：求定积分 ∫−44dx1+x2\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 xp=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+1xkt=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 1if a==1&&b==1S(1,1).root=Gauss_Legendre(-4,4,n,epsilon);endif 1S(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;elseS(a,b).d=0;Binary_tree(a+1,b*2-1,n,jd);endif(b+1<=2^(a-1)&&mod(b,2)==1)Binary_tree(a,b+1,n,jd);endend
end

运行结果如下：

效果还是很不错的。