馮海琴（赤峰學院　數(shù)學與統(tǒng)計學院，內蒙古　赤峰　024000）

?

基于Simpson公式的龍貝格求積算法

馮海琴
（赤峰學院數(shù)學與統(tǒng)計學院，內蒙古赤峰024000）

摘要：龍貝格求積算法屬于數(shù)值積分算法的一種，該算法的特點是精度高，計算方法簡單，收斂速度快.本文對基于辛普生公式的龍貝格算法進行了研究，設計了該算法的流程圖，并編寫了MATLAB程序，最后對該算法進行了仿真實驗，實驗結果說明了該算法的有效性.

關鍵詞：數(shù)值積分；辛普生公式；龍貝格算法

1　引言

常用的數(shù)值求積公式有梯形公式、辛普生公式及柯特斯公式.但是在很多時候利用這些低階的求積公式計算出的積分值并不能滿足精度要求，所以為了改善求積公式的精度，人們研究出一種行之有效的方法，即復化求積法.使用這種方法計算積分時事先選取大小合適的步長是關鍵，也是難點.解決這個問題的方法是采用變步長的計算方案.文獻［1］中給出的基于梯形法的龍貝格算法就是采用的這種方法.

由于梯形公式和辛普生公式是并行的求積算法，而且辛普生公式求積精度高于梯形公式，所以本文將對基于辛普生公式的龍貝格求積算法進行研究，以便得到收斂速度更快的求積算法.同時由于MATLAB軟件具有很強的數(shù)值計算功能.用MATLAB編寫的程序簡單，易于調試［3］～［6］.所以本文將對新設計的算法編寫MATLAB程序，最后對該算法進行實驗.

2　預備知識

梯形公式

辛普生公式

復化梯形公式

復化辛普生公式：

復化柯特斯公式：

復化梯形公式求積余項：

復化辛普生公式求積余項：

復化柯特斯公式求積余項：

3　基于Simpson公式的龍貝格求積算法

3.1辛普生公式的遞推化

設將求積的區(qū)間［a，b］分為n等分，則一共有n+1個等分點，n.這里用Sn表示用復化辛普生公式求得的積分值，其下標n表示等分數(shù).

對于小區(qū)間［xk，xk+1］，用S1與S2分別表示在該子段上二分前后的兩個積分值，其中，則有

顯然S1與S2有下列關系

將這一關系式關于k從0到n-1累加求和，即可導出下列遞推公式

這里的步長h代表的是二分之前的步長.

3.2龍貝格求積公式

雖然辛普生法較梯形法精度有所提高，收斂速度也明顯加快，但是對很多精度要求很高的數(shù)值積分問題仍然達不到要求，所以下面要研究其加速公式.根據(jù)復化辛普生公式求積余項（6）式知，當步長變?yōu)樵瓉聿介L的二分之一后，積分值S2n的誤差大約是積分值Sn誤差的，即有

將上式移項整理，知

在（8）式中，利用積分準確值I的兩個近似值Sn和S2n估計S2n的誤差的方法，這種方法稱為誤差的事后估計法，既然知道S2n的誤差近似等于，所以只要二分前后兩個積分值Sn與S2n相當接近，就可以保證計算結果S2n的誤差很小.如果將這個誤差值作為S2n的一種補償，期望得到I的精度更高的近似值，即

即有

然后我們使用同樣的方法繼續(xù)推到，根據(jù)復化柯特斯公式求積余項公式（8）知，步長二分之后的積分值C2n的誤差是步長二分之前的積分值Cn誤差的，既有

即為龍貝格值.

使用公式（10）和（11）作為加速公式，對于積分值Sn在步長不斷二分的過程中進行加速，從而有效的提高Sn的精度，得到精度較高的龍貝格值Rn，令Rn作為積分值的近似值，這就是基于辛普生公式的龍貝格算法.

3.3流程圖

基于辛普生公式的龍貝格求積算法的流程圖如圖1所示：

圖1　

3.4MATLAB源程序

本文對基于辛普生公式的龍貝格求積算法利用MATLAB進行了編程實現(xiàn)，程序如下：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仿真實驗

解建立函數(shù)文件：function y=f（x）

y=sin（x）/x；

再調用龍貝格算法：a=eps；

b=1；

ε=10∧（-7）；

r=romberg（a，b，ε）；

計算結果見表1：

表1　

從表中數(shù)據(jù)可知，只要將積分區(qū)間［0，1］二分2次，就可得到精度很高的近似值，而文獻［1］中使用了基于梯形法的龍貝格算法，也得到了較理想的結果，只是要將積分區(qū)間進行3次二分，這里少用了一次，節(jié)省了計算量.

4　總結

本文通過利用辛普生公式來推導出龍貝格求積算法，能夠加深對龍貝格求積的基本思路.利用辛普生公式推導出的龍貝格求積算法比梯形法推導出的龍貝格求積算法要節(jié)省很多計算量，大大的提高了精度，加速過程效果比較顯著，也便于應用程序的實現(xiàn).

參考文獻：

〔1〕王能超.數(shù)值分析簡明教程（第二版）［M］.北京：高等教育出版社，2008.

〔2〕李慶揚，王能超，易大義.數(shù)值分析［M］.武漢：華中科技大學出版社，1986.148～150.

〔3〕孫富玉，韓偉.MATLAB程序設計教程［M］.遠方出版社，2006.

〔4〕韓旭里.數(shù)值分析［M］.北京：高等教育出版社，2011.

〔5〕楊杰，趙曉暉.數(shù)學軟件與數(shù)學實驗［M］.北京：清華大學出版社，2011.

〔6〕牟古芳.數(shù)學實驗［M］.北京：高等教育出版社，2012.

中圖分類號：O241

文獻標識碼：A

文章編號：1673-260X（2016）06-0003-03

收稿日期：2016-02-25