matlab 龙贝格积分

直接上代码

clc;clear;
f=@(x)(1+(cos(x)).^2).^(1/2);
a=0;%下界
b=48;%上界
max=10;%迭代次数
R=zeros(max,max);
for row=1:max%算第一列
h=(b-a)/(2^(row-1));%步长
x1=a:h:b;
R(row,1)=(2*sum(f(x1))-f(a)-f(b))*h/2;
end
for col=2:max%递归公式
for row=col:max
    R(row,col)= R(row,col-1)+(R(row,col-1)-R(row-1,col-1))/(4^(col-1)-1);
end
end
disp(R);