数值分析 MATLAB实现Romberg积分法(最大N=16)

使用Romberg积分法求解\int_{1}^{3}\frac{100}{x^{2}}sin\frac{10}{x}dx

 

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

积分方法为先使用复化梯形公式求解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:最终结论也没手算过,(不一定对

  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值