引子

问题：如何计算一些复杂的定积分？
比如

Newton-cotes公式

  Newton-cotes(牛顿-科茨)公式是通过使用插值多项式近似代替被积函数进行积分，从而得到原积分的估计值。
  设要计算的积分为 。首先在积分域[a,b]上均匀地取N+1个点的计算出插值多项式 ，然后使用该插值多项式代替 在积分域上进行积分，得到估计值。对于取N+1个点计算出的插值多项式，其积分估计值可以直接用公式得到：


根据取点数多少，可以将Newton-cotes公式分为梯形法则、辛普森1/3法则、辛普森3/8法则等，下表给出了不同取点数下的积分估计值以及误差：


理论上取点数越多，插值多项式越接近原函数，估计值也越准确，但是当N+1≥8时，科茨系数Cn出现负数，导致积分结果不稳定。因此Newton-cotes公式的精度存在上限，对于一些曲线较为复杂的原函数，使用Newton-cotes公式计算误差较大。
Matlab代码:

% ======================================================================================
% 作者： cx
% 时间： 2025-04-18
% 实现： Newtion-cotes公式计算定积分
% ======================================================================================
function [R]=trapezoid_integral(f,a,b) 
% 梯形法则
%f：被积函数
%a：积分下限
%b：积分上限
%R：积分结果
R=(b-a)*(f(a)+f(b))/2;
end

function [R]=simpson3point_integral(f,a,b) 
% simpson1/3法则
%f：被积函数
%a：积分下限
%b：积分上限
%R：积分结果
R=(b-a)*(f(a)+4*f((a+b)/2)+f(b))/6;
end

function [R]=simpson4point_integral(f,a,b) 
% simpson3/8法则
%f：被积函数
%a：积分下限
%b：积分上限
%R：积分结果
h=(b-a)/3;
R=(b-a)*(f(a)+3*f(a+h)+3*f(a+2*h)+f(b))/8;
end

function [R]=boole_integral(f,a,b) 
% 布尔法则
%f：被积函数
%a：积分下限
%b：积分上限
%R：积分结果
h=(b-a)/4;
R=(b-a)*(7*f(a)+32*f(a+h)+12*f(a+2*h)+32*f(a+3*h)+7*f(b))/90;
end

%% 测试
clear;clc;close all;
% y=@(x) 4/3.*(4./ (((cos(x)-sin(x)/sqrt(3)).^2).*cos(x).^2)); %定义被积函数
% a=0;b=pi/12;
y=@(x) x.^2.*sin(x)+x.^3.*cos(x); %定义被积函数
a=0;b=5;                          %对于这种较大的积分区间，Newtion-cotes公式的计算精度会明显下降

R1=integral(y,a,b);
R2=trapezoid_integral(y,a,b);
R3=simpson3point_integral(y,a,b);
R4=simpson4point_integral(y,a,b);
R5=boole_integral(y,a,b);

disp("matlab内置积分函数结果为:"+R1);
disp("梯形法则积分结果为:"+R2+" ,相对误差为:"+(R1-R2)/R1);
disp("simpson1/3法则积分结果为:"+R3+" ,相对误差为:"+(R1-R3)/R1);
disp("simpson3/8法则积分结果为:"+R4+" ,相对误差为:"+(R1-R4)/R1);
disp("布尔法则积分结果为:"+R5+" ,相对误差为:"+(R1-R5)/R1);

变步长梯形方法

为解决Newton-cotes公式在计算一些复杂函数积分时无法收敛的问题，可以将积分区域[a,b]分割为多个子区间，在每个子区间中用低阶的Newton-cotes公式计算积分估计值，最后将所有子区间的估计值相加得到最终估计值。
假设使用梯形法则计算子区间的估计值，并通过迭代逐次增加子区间的数量，使最终估计值逼近真实值，这样便得到了变步长梯形方法：


Matlab代码:

% ======================================================================================
% 作者： cx
% 时间： 2025-04-18
% 实现： 变步长梯形公式计算定积分
% ======================================================================================
function [R,I,S]=Iterative_trapezoid_integral(f,a,b,eps) 
% 变步长梯形公式
%f：被积函数
%a：积分下限
%b：积分上限
%eps：要求精度
%R：积分结果
%I：迭代次数
%S：每次迭代估计值

if nargin < 4
    eps = 1e-4; % 默认精度
end

i=1;        %迭代计数
N=2^(i-1);  %划分区间数
h=(b-a)/N;  %步长
err=1;
S(i)=h/2*(f(a)+f(b));
while( err > eps )
    i=i+1;
    N=2^(i-1);
    h=(b-a)/N;
    S(i)=S(i-1)/2;
    for n=1:2:N-1
      S(i)=S(i)+h*f(a+n*h);
    end
    err=abs((S(i)-S(i-1))/S(i));
end
R=S(i);
I=i;
end

%% 测试
clear;clc;close all;
% y=@(x) 4/3.*(4./ (((cos(x)-sin(x)/sqrt(3)).^2).*cos(x).^2)); %定义被积函数
% a=0;b=pi/12;
y=@(x) x.^2.*sin(x)+x.^3.*cos(x); %定义被积函数
a=0;b=5;

R1=integral(y,a,b);
[R2,I,S]=Iterative_trapezoid_integral(y,a,b,1e-4);

disp("matlab内置积分函数结果为:"+R1);
disp("变步长梯形法则积分结果为:"+R2+",划分区间数:"+2^(I-1)+",迭代次数:"+I+" ,相对误差为:"+(R1-R2)/R1);

Romberg方法

变步长梯形方法的精度可以满足要求，但是迭代收敛速度还不够快，在其基础上增加外推加速算法以加快收敛速度，就得到了Romberg(龙贝格)方法。
外推加速算法介绍：在变步长梯形方法的第i次和第i+1次迭代中，可以将积分真实值 与估计值 的关系表示如下：

于是可以构造一个新的估计值：

可以看到新估计值的误差小于原来的两个估计值，这样就完成了一次外推加速，即使用低精度的估计值得到高精度的估计值。接着可以继续外推:
对于第i次迭代，总共可以外推i-1次，并将误差量级从 减少到 ，从而大大提高收敛速度，其中第j次外推的计算公式为:

综上，Romberg方法的计算步骤如下：


Matlab代码:

% ======================================================================================
% 作者： cx
% 时间： 2025-04-18
% 实现： Romberg方法计算定积分
% ======================================================================================
function [R,I,ST]=Romberg_integral(f,a,b,eps) 
% 变步长梯形公式
%f：被积函数
%a：积分下限
%b：积分上限
%eps：要求精度
%R：积分结果
%I：迭代次数
%ST：每次迭代估计值、外推估计值

if nargin < 4
    eps = 1e-4; % 默认精度
end

i=1;        %迭代计数
N=2^(i-1);  %划分区间数
h=(b-a)/N;  %步长
err=1;
S=zeros(100,1);
T=zeros(100,100);
S(1)=h/2*(f(a)+f(b));
while( err > eps || i>100)
    i=i+1;
    N=2^(i-1);
    h=(b-a)/N;
    S(i)=S(i-1)/2;
    for n=1:2:N-1
       S(i)=S(i)+h*f(a+n*h);  %变步长梯形公式
    end
    T(i,1)=(4*S(i)-S(i-1))/3;
    if i>=3
        for j=2:1:i-1
          T(i,j)=(4^j*T(i,j-1)-T(i-1,j-1))/(4^j-1);  %外推加速
        end
        err=abs((T(i,i-1)-T(i,i-2))/T(i,i-2));
    end
end
R=T(i,i-1);
I=i;
ST=[S T];
end

%% 测试
clear;clc;close all;
% y=@(x) 4/3.*(4./ (((cos(x)-sin(x)/sqrt(3)).^2).*cos(x).^2)); %定义被积函数
% a=0;b=pi/12;
y=@(x) x.^2.*sin(x)+x.^3.*cos(x); %定义被积函数
a=0;b=5;

R1=integral(y,a,b);
[R2,I,ST]=Romberg_integral(y,a,b,1e-4);

disp("matlab内置积分函数结果为:"+R1);
disp("Romberg方法积分结果为:"+R2+",划分区间数:"+2^(I-1)+",迭代次数:"+I+" ,相对误差为:"+(R1-R2)/R1);