使用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:最终结论也没手算过,(不一定对