数值分析 MATLAB实现Romberg积分法（最大N=16）

使用Romberg积分法求解

 

其他积分函数可以按注释对代码进行修改。

积分方法为先使用复化梯形公式求解T序列，再使用Romberg法求解复化simpson，即S序列，再使用Romberg法求解复化cotes，即C序列，最后同理求R序列。

代码支持N=1，2，4，8，16;再添加需要按模块手动添加。

function Romberg_my
x1=1:2;
x2=1:3;
x3=1:5;
x4=1:9;
x5=1:17;
x0=1;%x0，xn为积分上下限
xn=3;

x1=[x0,xn];
x2(1)=x0;
x2(3)=xn;
x3(1)=x0;
x3(5)=xn;
x4(1)=x0;
x4(9)=xn;
x5(1)=x0;
x5(17)=xn;

%N=1
for i=1:2
    y1(i)=(100/(x1(i)*x1(i)))*sin(10/x1(i));%被积函数
end
h1=(xn-x0)/2^0;
t11=h1/2*(y1(1)+y1(2));


%N=2
for i=2
    x2(i)=0.5*(x1(i-1)+x1(i));
end
for i=1:3
    y2(i)=(100/(x2(i)*x2(i)))*sin(10/x2(i));%被积函数
end
t20=y2(1)+y2(3);
for i=2
    t20=t20+2*y2(i);
end
h2=(xn-x0)/2^1;
t21=h2/2*t20;

%N=4
for i=1:2
    x3(2*i)=0.5*(x2(i)+x2(i+1));
    x3(2*i-1)=x2(i);
end
for i=1:5
    y3(i)=(100/(x3(i)*x3(i)))*sin(10/x3(i));%被积函数
end
t30=y3(1)+y3(5);
for i=2:4
    t30=t30+2*y3(i);
end
h3=(xn-x0)/2^2;
t31=h3/2*t30;

%N=8
for i=1:4
    x4(2*i)=0.5*(x3(i)+x3(i+1));
    x4(2*i-1)=x3(i);
end
for i=1:9
    y4(i)=(100/(x4(i)*x4(i)))*sin(10/x4(i));%被积函数
end
t40=y4(1)+y4(9);
for i=2:8
    t40=t40+2*y4(i);
end
h4=(xn-x0)/2^3;
t41=h4/2*t40;

%N=16
for i=1:8
    x5(2*i)=0.5*(x4(i)+x4(i+1));
    x5(2*i-1)=x4(i);
end
for i=1:17
    y5(i)=(100/(x5(i)*x5(i)))*sin(10/x5(i));%被积函数
end
t50=y5(1)+y5(17);
for i=2:16
    t50=t50+2*y5(i);
end
h5=(xn-x0)/2^4;
t51=h5/2*t50;


t=[t11,t21,t31,t41,t51];
for i=1:4
    s(i)=(4^1*t(i+1)-t(i))/(4^1-1);
end
for i=1:3
    c(i)=(4^2*s(i+1)-s(i))/(4^2-1);
end
for i=1:2
    r(i)=(4^3*c(i+1)-c(i))/(4^3-1);
end

t
s
c
r

end

PS：最终结论也没手算过，（不一定对