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

馮海琴（赤峰學(xué)院　數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，內(nèi)蒙古　赤峰　024000）

?

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

馮海琴
（赤峰學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，內(nèi)蒙古赤峰024000）

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

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

1　引言

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

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

2　預(yù)備知識

梯形公式

辛普生公式

復(fù)化梯形公式

復(fù)化辛普生公式：

復(fù)化柯特斯公式：

復(fù)化梯形公式求積余項(xiàng)：

復(fù)化辛普生公式求積余項(xiàng)：

復(fù)化柯特斯公式求積余項(xiàng)：

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

3.1辛普生公式的遞推化

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

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

顯然S1與S2有下列關(guān)系

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

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

3.2龍貝格求積公式

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

將上式移項(xiàng)整理，知

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

即有

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

即為龍貝格值.

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

3.3流程圖

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

圖1　

3.4MATLAB源程序

本文對基于辛普生公式的龍貝格求積算法利用MATLAB進(jìn)行了編程實(shí)現(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í)驗(yàn)

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

y=sin（x）/x；

再調(diào)用龍貝格算法：a=eps；

b=1；

ε=10∧（-7）；

r=romberg（a，b，ε）；

計(jì)算結(jié)果見表1：

表1　

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

4　總結(jié)

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

參考文獻(xiàn)：

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

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

〔3〕孫富玉，韓偉.MATLAB程序設(shè)計(jì)教程［M］.遠(yuǎn)方出版社，2006.

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

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

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

中圖分類號：O241

文獻(xiàn)標(biāo)識碼：A

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

收稿日期：2016-02-25