冯海琴（赤峰学院　数学与统计学院，内蒙古　赤峰　024000）



基于Simpson公式的龙贝格求积算法

冯海琴
（赤峰学院数学与统计学院，内蒙古赤峰024000）

摘要：龙贝格求积算法属于数值积分算法的一种，该算法的特点是精度高，计算方法简单，收敛速度快.本文对基于辛普生公式的龙贝格算法进行了研究，设计了该算法的流程图，并编写了MATLAB程序，最后对该算法进行了仿真实验，实验结果说明了该算法的有效性.

关键词：数值积分；辛普生公式；龙贝格算法

1　引言

常用的数值求积公式有梯形公式、辛普生公式及柯特斯公式.但是在很多时候利用这些低阶的求积公式计算出的积分值并不能满足精度要求，所以为了改善求积公式的精度，人们研究出一种行之有效的方法，即复化求积法.使用这种方法计算积分时事先选取大小合适的步长是关键，也是难点.解决这个问题的方法是采用变步长的计算方案.文献［1］中给出的基于梯形法的龙贝格算法就是采用的这种方法.

由于梯形公式和辛普生公式是并行的求积算法，而且辛普生公式求积精度高于梯形公式，所以本文将对基于辛普生公式的龙贝格求积算法进行研究，以便得到收敛速度更快的求积算法.同时由于MATLAB软件具有很强的数值计算功能.用MATLAB编写的程序简单，易于调试［3］～［6］.所以本文将对新设计的算法编写MATLAB程序，最后对该算法进行实验.

2　预备知识

梯形公式

辛普生公式

复化梯形公式

复化辛普生公式：

复化柯特斯公式：

复化梯形公式求积余项：

复化辛普生公式求积余项：

复化柯特斯公式求积余项：

3　基于Simpson公式的龙贝格求积算法

3.1辛普生公式的递推化

设将求积的区间［a，b］分为n等分，则一共有n+1个等分点，n.这里用Sn表示用复化辛普生公式求得的积分值，其下标n表示等分数.

对于小区间［xk，xk+1］，用S1与S2分别表示在该子段上二分前后的两个积分值，其中，则有

显然S1与S2有下列关系

将这一关系式关于k从0到n-1累加求和，即可导出下列递推公式

这里的步长h代表的是二分之前的步长.

3.2龙贝格求积公式

虽然辛普生法较梯形法精度有所提高，收敛速度也明显加快，但是对很多精度要求很高的数值积分问题仍然达不到要求，所以下面要研究其加速公式.根据复化辛普生公式求积余项（6）式知，当步长变为原来步长的二分之一后，积分值S2n的误差大约是积分值Sn误差的，即有

将上式移项整理，知

在（8）式中，利用积分准确值I的两个近似值Sn和S2n估计S2n的误差的方法，这种方法称为误差的事后估计法，既然知道S2n的误差近似等于，所以只要二分前后两个积分值Sn与S2n相当接近，就可以保证计算结果S2n的误差很小.如果将这个误差值作为S2n的一种补偿，期望得到I的精度更高的近似值，即

即有

然后我们使用同样的方法继续推到，根据复化柯特斯公式求积余项公式（8）知，步长二分之后的积分值C2n的误差是步长二分之前的积分值Cn误差的，既有

即为龙贝格值.

使用公式（10）和（11）作为加速公式，对于积分值Sn在步长不断二分的过程中进行加速，从而有效的提高Sn的精度，得到精度较高的龙贝格值Rn，令Rn作为积分值的近似值，这就是基于辛普生公式的龙贝格算法.

3.3流程图

基于辛普生公式的龙贝格求积算法的流程图如图1所示：

图1　

3.4MATLAB源程序

本文对基于辛普生公式的龙贝格求积算法利用MATLAB进行了编程实现，程序如下：function r=romberg（a，b，ε）

h=b-a；

c=1/2*（a+b）；

s1=h/6*（f（a）+4*f（c）+f（b）；

x=a；

s=0；

while x＜b

x1=x+h/4；

x2=x1+h/4；

x3=x2+h/4；

s=s+2*f（x1）-f（x2）+2*f（x3）；

x=x+h；

end

s2=1/2*s1+h/6*s；

c1=s2+1/15*（s2-s1）；

h=h/2；

s1=s2；

x=a；

s=0；

while x＜b

x1=x+h/4；

x2=x1+h/4；

x3=x2+h/4；

s=s+2*f（x1）-f（x2）+2*f（x3）；

x=x+h；

end

s2=1/2*s1+h/6*s；

c2=s2+1/15*（s2-s1）；

r1=c2+1/63*（c2-c1）；

h=h/2；

s1=s2；

c1=c2；

x=a；

s=0；

while x＜b

x1=x+h/4；

x2=x1+h/4；

x3=x2+h/4；

s=s+2*f（x1）-f（x2）+2*f（x3）；

x=x+h；

end

s2=1/2*s1+h/6*s；

c2=s2+1/15*（s2-s1）；

r2=c2+1/63*（c2-c1）；

while abs（r2-r1）＞=

h=h/2；

s1=s2；

c1=c2；

r1=r2；

x=a；

s=0；

while x＜b

x1=x+h/4；

x2=x1+h/4；

x3=x2+h/4；

s=s+2*f（x1）-f（x2）+2*f（x3）；

x=x+h；

end

s2=1/2*s1+h/6*s；

c2=s2+1/15*（s2-s1）；

r2=c2+1/63*（c2-c1）；

end

r=r2；

3.5仿真实验

解建立函数文件：function y=f（x）

y=sin（x）/x；

再调用龙贝格算法：a=eps；

b=1；

ε=10∧（-7）；

r=romberg（a，b，ε）；

计算结果见表1：

表1　

从表中数据可知，只要将积分区间［0，1］二分2次，就可得到精度很高的近似值，而文献［1］中使用了基于梯形法的龙贝格算法，也得到了较理想的结果，只是要将积分区间进行3次二分，这里少用了一次，节省了计算量.

4　总结

本文通过利用辛普生公式来推导出龙贝格求积算法，能够加深对龙贝格求积的基本思路.利用辛普生公式推导出的龙贝格求积算法比梯形法推导出的龙贝格求积算法要节省很多计算量，大大的提高了精度，加速过程效果比较显著，也便于应用程序的实现.

参考文献：

〔1〕王能超.数值分析简明教程（第二版）［M］.北京：高等教育出版社，2008.

〔2〕李庆扬，王能超，易大义.数值分析［M］.武汉：华中科技大学出版社，1986.148～150.

〔3〕孙富玉，韩伟.MATLAB程序设计教程［M］.远方出版社，2006.

〔4〕韩旭里.数值分析［M］.北京：高等教育出版社，2011.

〔5〕杨杰，赵晓晖.数学软件与数学实验［M］.北京：清华大学出版社，2011.

〔6〕牟古芳.数学实验［M］.北京：高等教育出版社，2012.

中图分类号：O241

文献标识码：A

文章编号：1673-260X（2016）06-0003-03

收稿日期：2016-02-25