利用Matlab实现Romberg数值积分算法系统建模与仿真结课作业.docx
《利用Matlab实现Romberg数值积分算法系统建模与仿真结课作业.docx》由会员分享，可在线阅读，更多相关《利用Matlab实现Romberg数值积分算法系统建模与仿真结课作业.docx（8页珍藏版）》请在冰豆网上搜索。
利用Matlab实现Romberg数值积分算法系统建模与仿真结课作业
利用Matlab实现Romberg数值积分算法
一、内容摘要
针对于某些多项式积分，利用Newton—Leibniz积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用Matlab中的.m文件编写了复化梯形公式与Romberg的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。
二、数值积分公式
1.复化梯形公式求解数值积分的基础是将区间一等分时的Newton—Cotes求积公式:
=
其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似
在区间[a,b]上的积分值，截断误差为：
具有一次的代数精度，很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：
=
其截断误差为：
数值积分算法
使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg数值积分。
其思想主要是，根据
的近似值
加上
与
的近似误差，作为新的
的近视，反复迭代，求出满足计算精度的近似解。
用
近似
所产生的误差可用下式进行估算：
新的
的近似值：
=（012….）
Romberg数值积分算法计算顺序
i=0
（1）
i=1
（2）
（3）
i=2
（4）
（5）
（6）
i=3
（7）
（8）
（9）
（10）
i=4
（11）
（12）
（13）
（14）
…
…
…
…
其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即Romberg序列。
三、复化梯形法以及Romberg算法程序流程图
图1复化梯形法程序流程图
图2Romberg算法程序流程图
四、计算实例
依据上文所述的流程图，编写复化梯形程序以及Romberg算法程序，并且利用实例验证程序的正确性，示例如下（计算精度）：
表2计算结果
计算精度
×10^-5
×10^-7
×10^-9
复化梯形
算法
时间
近似值
3.
3.
3.
Romberg
算法
时间
近似值
3.
3.
3.
从上表中可以看出，当要求的计算精度不高时，复化梯形算法与Romberg算法计算时间相差不太大，但是Romberg算法是要快于复化梯形算法的；当要求的计算精度更高的时候，Romberg算法是明显快于复化梯形算法。
本文所编写的程序适用于多项式的数值积分，且对于积分区间内，被积函数在每一点必须有定义，在以后的学习中进一步改进。
附录：
1.复化梯形算法程序
function[]=sf（a,b,m,M,d）
tic
disp（'请输入分子多项式a，分母多项式b，积分下限m,积分上限M,以及计算精度d'）
f=poly2sym（a）/poly2sym（b）%用于给用户显示被积函数的形式
%利用梯形公式计算此数值积分
disp（'利用梯形公式计算数值积分的结果'）
kk=zeros（）;%用于存放结果
kk（1,1）=1/2*（M-m）/1*（subs（f,'x',m）+subs（f,'x',M））%先存储首项
fori=1:
1:
2^30
t=0;
forj=0:
1:
2^（i-1）-1
v=m+（2*j+1）*（M-m）/（2^i）
vv=polyval（a,v）/polyval（b,v）;
t=t+（M-m）/（2^i）*vv
end
y=1/2*kk（i,1）+t%通项公式计算各项值
kk（i+1,1）=y%存储其他项
f=i+1;%记录符合条件的值的下标
if（1/3*（kk（i+1,1）-kk（i,1））<=d）
break;
end
end
time=toc
fprintf（'Theresultis%f\n',kk（f,1））
算法程序
function[]=romberg（a,b,m,M,d）
tic
disp（'请输入分子多项式a，分母多项式b，积分下限m,积分上限M,以及计算精度d'）
f=poly2sym（a）/poly2sym（b）%用于给用户显示被积函数的形式
disp（'利用梯形公式计算数值积分的结果'）
kk=zeros（）;%用于存放结果
kk（1,1）=1/2*（M-m）/1*（subs（f,'x',m）+subs（f,'x',M））;%先存储首项
fori=1:
1:
2^40
t=0;
forj=0:
1:
2^（i-1）-1
v=m+（2*j+1）*（M-m）/（2^i）;
vv=polyval（a,v）/polyval（b,v）;
t=t+（M-m）/（2^i）*vv;
end
y=1/2*kk（i,1）+t;%通项公式计算各项值
kk（i+1,1）=y;%存储其他项
if（abs（1/3*（kk（i+1,1）-kk（i,1）））<=d）%判断梯形公式值是否达到要求
disp（'Theresultis:
'）
kk（）
kk（i+1,1）%梯形值满足要求，输出结果
break;
else
s=（4*kk（i+1,1）-kk（i,1））/（4-1）;%构造simpson各项
kk（i+1,2）=s%存储
if（i+1>=3）
if（i+1>=3&abs（1/15*（kk（i+1,2）-kk（i,2）））<=d）
kk（）
disp（'Theresultis:
'）
kk（i+1,2）%simpson值满足要求，输出结果
pan1=0;
break;
else
c=（4^2*kk（i+1,2）-kk（i,2））/（4^2-1）;%构造cotes值
kk（i+1,3）=c%存储cotes值
if（i+1>=4）
if（i+1>=4&abs（1/63*（kk（i+1,3）-kk（i,3）））<=d）
disp（'Theresultis:
'）
kk（i+1,3）
break;
else
r=（4^3*kk（i+1,3）-kk（i,3））/（4^3-1）%构造romberg值
kk（i+1,4）=r%存储romberg值
if（i+1>=5）
if（i+1>=5&abs（1/127*（kk（i+1,4）-kk（i,4）））<=d）
disp（'Theresultis:
'）
kk（i+1,4）
break;
end
end
end
end
end
end
end
end
time=toc
升级会员 