自适应梯形积分算法

自适应梯形积分算法

引言：笔者为学习数值分析的苦逼大二学生，为了做展示需要用到自适应梯形积分算法，但是遍历我的书本和百度也没能找到，故自己根据自适应辛普森公式推导之。

一、梯形公式和辛普森公式

对于数值积分，我们最容易想到的就是梯形公式了，也就是取该段函数的起点终点的平均值作为整段函数的估计值。就像下图所示：

书本上的公式肯定是对的（\狗头)

二、复化梯形公式

有时候，为了提高精度，而不愿意改变算法，我们可以采用复化的形式，也就是把函数分成多段，每一段用梯形公式进行计算得到该段的估计值，再进行累加。（下图为书本上的公式）

三、自适应复合梯形公式推导

但是，我们往往碰到的函数可能一段上很均匀，但是在另一段上波动很大，例如我们这里测试用的函数（下图为在定义域为$(1,3)$上的形状）。
$$F=(100/x^2)\ast (sin(10/x));$$

这时候我们的自适应梯形公式就可以派上用场了，它的实现手段就是在波动较大的地方多进行分段，而波动较小的地方则较少的分段足以达到精度要求。而如何检验波动大小就成了问题的痛点所在，下面我们重点放在如何分析现有估计的误差，并作出合理的判断是否继续细分。

1.推导

先利用梯形公式：
$$\begin{gathered} {\int }_a^{b}f(x)dx=S(a,b)-\frac{h^3}{12}{f}^{\prime \prime }(\xi ) \\ \text{\hspace{1em}(1)} \end{gathered}$$
其中，
$$\begin{gathered} S(a,b)=\frac{h}{2}[f(a)+f(b)] \\ \text{\hspace{1em}(2)} \end{gathered}$$

$$Then,applytheCompositeTrapezoidalrulewithn=3andstepsize(b-a)=h$$

$$\begin{gathered} {\int }_a^{b}f(x)dx=\frac{h}{2}[f(a)+2\times f(a+\frac{h}{2})+f(b)]-(\frac{h}{2}{)}^2\frac{b-a}{12}{f}^{\prime \prime }(\bar{\xi }) \\ =S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)-\frac{h^3}{48}{f}^{\prime \prime }(\bar{\xi }) \\ \text{\hspace{1em}(3)} \end{gathered}$$

为了计算误差，我们做出以下假设：
$$wesuppose{f}^{\prime \prime }(\bar{\xi })\approx f^{\prime \prime }(\xi )$$
于是有，
$$\begin{gathered} S(a,b)-\frac{h^3}{12}{f}^{\prime \prime }(\xi )=S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)-\frac{h^3}{48}{f}^{\prime \prime }(\xi ) \\ \text{\hspace{1em}(4)} \end{gathered}$$

$$\begin{gathered} S(a,b)-[S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)]=\frac{h^3}{12}{f}^{\prime \prime }(\xi )-\frac{h^3}{48}{f}^{\prime \prime }(\xi ) \\ \text{\hspace{1em}(5)} \end{gathered}$$

$$S(a,b)-[S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)]\approx \frac{h^3}{16}{f}^{\prime \prime }(\xi )\text{\hspace{1em}(5)}$$

因此误差可以估计为,
$$|{\int }_a^{b}f(x)dx-S(a,\frac{a+b}{2})-S(\frac{a+b}{2},b)|\approx \frac{1}{3}|S(a,b)-S(a,\frac{a+b}{2})-S(\frac{a+b}{2},b)|\text{\hspace{1em}(6)}$$

这个结果表明,我们用$S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)$进行估计积分值的大小的误差为$\frac{1}{3}|S(a,b)-S(a,\frac{a+b}{2})-S(\frac{a+b}{2},b)|$,我们就可以用函数本身的性质得到计算误差的大小了，故称之为自适应(adapted)。

2、自适应复合梯形公式的matlab 代码

定义精度要求为$1\times 10^{-3}$

format long
f1 = @(x) (100/x^2)*(sin(10/x));	% 测试函数(the test function)
% 积分上下限和精度要求，注意这里是整体的精度要求
a=1;
b=3;
e=1e-3;

p=[a,b];
p0=p;
ep=[e];
m=0;	%分点数
q=0;
I=0;	%最后的结果
while(1)
    n1=length(ep);
    n=length(p0);
    if n==1
        break;
    end
    h=p0(2)-p0(1);
    whole=h/2*(f1(p0(1))+f1(p0(1)+h));
    left=h/4*(f1(p0(1))+f1(p0(1)+h/2));
    right=h/4*(f1(p0(1)+h/2)+f1(p0(1)+h));
    if abs(whole-left-right)<=3*ep(1)	%精度满足要求，则该部分不再细分
        I=I+left+right;
        p0=p0(2:n);
        if n1>=2
            ep=ep(2:n1);
        end
        q=q+1;
    else								%精度不满足要求，继续细分
        m=m+1;	%分层数++
        p0=[p0(1),p0(1)+h/2,p0(2:n)];
        if n1==1
            ep=[ep(1)/2,ep(1)/2];	%第一次分段，将原有的精度要求提高，两段分别不超过ep/2
        else
            ep=[ep(1)/2,ep(1)/2,ep(2:n1)];%分段，将原有的精度提高，两段分别不超过ep/2
        end
        if q==0
            p=p0;
        else
            p=[p(1:q),p0];
        end
    end
end
fprintf('求得的解为：%f\n',I);
fprintf('分层数为%d\n',m);
p=p';
disp('分点为');p

3、与一般复化梯形公式的比较

clc;
clear;
f=@(x) (100/x^2)*(sin(10/x));
result = -1.426014;
 a=input('Enter lower limit a: '); % exmple a=1
 b=input('Enter upper limit b: ');  % exmple b=3
 e=input('Enter the eps: ');  % exmple eps=10^-3
 n=2;
 MAXN = 3000;	%define the maximum number of segments
 while n<MAXN
    sum=0;
    h=(b-a)/n;
    for k=1:1:n-1
      x(k)=a+k*h;
      y(k)=f(x(k));
      sum=sum+y(k);
    end
    % Formula:  (h/2)*[(y0+yn)+2*(y2+y3+..+yn-1)]
    answer=h/2*(f(a)+f(b)+2*sum);
    if (abs(answer-result)<=e)	%the precision of the result is satisfied.
        break;
    end
    n = n + 1;
 end

fprintf('\n The value of integration is %f',answer); 
fprintf('\n The number of intermediate points is %f',n);  

结论：

可以看出，自适应梯形公式计算的点数要比一般复化梯形公式少，从这个角度上看，自适应梯形公式更加适合计算该函数的数值积分。