复化梯形、复化辛普森、复化科特斯、龙贝格MATLAB实现

%下面除了龙贝格外，其他均以此fun函数作为被积函数
%梯形、辛普森、科特斯：都是已知积分上下限和分段数，求解积分近似值。
function a = fun(x)
    a = cos(x*x);
end

一、复化梯形：

function Tn = ComplexTrap(xL, xR, n)
% xL,xR：积分区间左右端点，n：分段数
% Tn ：复化梯形积分结果
h = (xR-xL)/n;%求步长
Tn = 0;
x = xL:h:xR;
for k = 1:n
    Tn = Tn + fun(x(k)) + fun(x(k+1));
end 
Tn = Tn*h/2;
Tn = vpa(Tn,6);
end 

二、复化辛普森：

function Sn = ComplexSimpson(xL, xR, n)
% xL,xR：积分区间左右端点，n：分段数
% Sn ：复化辛普森积分结果
h = (xR-xL)/n;%求步长
Sn = 0;
x = xL:h:xR;
for k = 2:n+1
    Sn = Sn + 4*fun((x(k)-0.5)) + 2*fun(x(k)); 
end 
Sn = Sn + fun(xL) - fun(xR);
Sn = vpa(Sn,7);
Sn = Sn*h/6;
Sn = vpa(Sn,6);
end 

三、复化科斯特：

function Cn = ComplexCotes(xL, xR,n)
% 判断积分区间个数是否是2的倍数，满足则进行计算，否则打印提示
if mod(n,4) ~= 0
        print('等分区间数错误！');return
else
    % 获取步长h
    h = (xR-xL)/n;
    %累积计算
    Cn = 0;
    for i = 1:n-3
        Cn = Cn + 7*fun(xL+h*(i-1))...
                 +32*fun(xL+h*i)...
                 +12*fun(xL+h*(i+1))...
                 +32*fun(xL+h*(i+2))...
                 +7*fun(xL+h*(i+3));
    end % 循环结束
    Cn = Cn*h/90;
    Cn = vpa(Cn,6);
end % if判断结束
end % 函数结束

四、龙贝格：

function I = romberg(fun,a,b,e)
% 使用龙贝格(Romberg数值求解公式)
% 输入例如：I=romberg(@(x)x^(3/2),0,1,0.000001)
% 判断输入参数是否足够
if nargin~=4
    error('请输入需要求积分的f、上界和下界以及误差e')
end
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a; %上下界间距
T(1,1)=h/2*(fun(a)+fun(b));
d=b-a; %误差间距
while e<=d
    k=k+1;
    h=h/2;
    sum=0;
    % 计算第一列T
    for i=1:n
        sum=sum+fun(a+(2*i-1)*h);
    end
    T(k+1,1)=T(k)/2+h*sum;
    % 迭代
    for j=1:k
        T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
    end
    n=n*2;
    d=abs(T(k+1,k+1)-T(k,k));
end
T 
I=T(k+1,k+1);