楊錄峰,馬 寧,趙雙鎖,3

(1.北方民族大學(xué)信息與計(jì)算科學(xué)學(xué)院,寧夏銀川 750021;2.吳忠市氣象局,寧夏吳忠 751300;3.蘭州大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,甘肅蘭州 730001)

一種變步長和變階計(jì)算的自適應(yīng)數(shù)值積分算法

楊錄峰1,馬 寧2,趙雙鎖1,3

(1.北方民族大學(xué)信息與計(jì)算科學(xué)學(xué)院,寧夏銀川 750021;2.吳忠市氣象局,寧夏吳忠 751300;3.蘭州大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,甘肅蘭州 730001)

將自適應(yīng) Simpson算法和 Romberg外推算法相結(jié)合,提出一種新型的自適應(yīng) S-R(S impson-Romberg)算法,它兼有變步長計(jì)算和逐步提高數(shù)值積分法收斂階的優(yōu)點(diǎn).若干數(shù)值比較算例表明,當(dāng)被積函數(shù)在積分區(qū)間上變化性態(tài)急劇多變時,與自適應(yīng) Simpson算法和 Romberg外推算法相比,它具有明顯優(yōu)勢.

數(shù)值積分法;自適應(yīng) Simpson算法;Romberg外推算法;自適應(yīng) S-R算法

數(shù)值積分

以自適應(yīng) Simpson算法[1-3]為優(yōu)秀代表的一類自適應(yīng)算法的主要優(yōu)點(diǎn)是能根據(jù) f(x)在積分區(qū)間[a,b]內(nèi)變化的不同性態(tài)自動選擇積分步長,因而在區(qū)間[a,b]上,該算法是變步長計(jì)算的.但這種算法在[a,b]上始終采用同一求積公式,故收斂階是固定不變的,不管 f(x)多么光滑也是如此.數(shù)值積分(1)的另一類數(shù)值方法是外推算法,它在積分過程中可以不斷使用由Richardson外推技術(shù)[4-5]得到的收斂階更高的數(shù)值積分公式,只要 f(x)滿足外推公式要求的光滑性條件和端點(diǎn)條件且控制條件許可,充分利用f(x)的光滑性利用Richardson外推過程以盡可能提高收斂階,從而加快收斂速度,是這一類方法的優(yōu)點(diǎn)[6].但這一算法在[a,b]上卻是等步長的,它不能根據(jù) f(x)在[a,b]上不同子區(qū)間的不同變化性態(tài),選取最合理的變積分步長,除非事先把[a,b]劃分成若干個小區(qū)間,然后分別求其數(shù)值積分,最后求它們的和.然而對[a,b]事先完成這一剖分卻是困難的,因此,當(dāng)被積函數(shù) fx在[a,b]上變化性態(tài)是急劇多變時,該算法的有效性將受到限制,這一算法的優(yōu)秀代表是由梯形公式不斷外推的 Romberg算法[4].利用 Romberg求解數(shù)值積分的思想,MladenMestrovic與 Eva Ocvirk[7]將 Romberg算法應(yīng)用于第 2類Volterra積分方程的數(shù)值求解問題.其他類數(shù)值積分方法,如蒙特卡洛方法[8],基于有限元的積分方法[9]等均有較快的發(fā)展.

本文將自適應(yīng) Simpson算法和 Romberg外推算法相結(jié)合構(gòu)造了一種新型自適應(yīng)算法,稱之為自適應(yīng)S-R(Simpson-Romberg)算法.該算法既能根據(jù) f(x)在[a,b]上變化性態(tài)的多變性,對[a,b]自動剖分,又能充分利用 f(x)在[a,b]上的光滑性以盡可能提高收斂階.可見,S-R算法兼有變步長計(jì)算和逐步提高數(shù)值積分收斂階的優(yōu)點(diǎn),因而更具有有效性.對于變化性態(tài)改變不大的被積函數(shù) f(x)的數(shù)值積分,數(shù)值算例表明該算法與 Romberg算法同樣有效.

1 自適應(yīng) S-R算法

1.1 算法梗概

首先,在“低精度”ε1的控制下,用自適應(yīng) Simpson算法完成對積分區(qū)間[a,b]的剖分,使得 f(x)在剖分后所得到的每個小區(qū)間內(nèi)每一點(diǎn)上的變化性態(tài)是相近的;然后在控制精度ε的要求下,利用Romberg外推算法完成每個小區(qū)間上的積分;最后把每一個小區(qū)間上所得到的數(shù)值結(jié)果求和,便得到 (1)的高精度近似積分結(jié)果.

1.2 若干術(shù)語

記 a ＜ b,α ＜γ,a≤α,γ≤b,h=γ-α,xj=α+jh/4,j=0∶4,β =x2

其中 Pαγ和 Qαγ分別是對 Iαγ近似積分的 3點(diǎn)和 5點(diǎn) Simpson公式[6].所謂在 [α,γ]上做 Simpson試驗(yàn) ,就是指分別算出 Pαγ和 Qαγ并對作判斷.

易知,當(dāng)γ -α→0時,如記γ→ x＊,則σ → f′(x＊),故當(dāng)γ-α足夠小時,σ的大小在一定程度上反映了 f(x)在[α,γ]上每一點(diǎn)的變化激烈 (程)度,σ愈大,變化愈激烈,反之愈平緩.據(jù)此可以認(rèn)為,當(dāng)γ-α足夠小時,σ為函數(shù) f(x)在區(qū)間[α,γ]上每一點(diǎn)的變化激烈度的一個近似描述.

選取“低精度 ”ε0,定義

這里ε0是根據(jù)對 I近似積分的精度要求選取的.一般說來,若精度要求較低,則取ε0較大,反之較小,通常取ε0∈[10-1,10-2].考慮到,對 I近似積分精度要求實(shí)現(xiàn)的好壞,是與 f(x)的變化性態(tài)有關(guān)的事實(shí),在本文中我們規(guī)定,當(dāng)σ ＜10時取ε0=10-1,而當(dāng)σ≥10時則取ε0=10-2,由程序自動完成.

為說明如何用自適應(yīng) Simpson算法對區(qū)間[a,b]實(shí)現(xiàn)剖分,先設(shè)在某個子區(qū)間 [α,γ]上用 Pαγ和 Qαγ近似計(jì)算 Iαγ.根據(jù)復(fù)合 Simpson公式的誤差有:

其中ξ,η∈[α,γ],兩式相減,略去高階項(xiàng) O(h6)整理得

故

根據(jù)文獻(xiàn)[6],如果要求

則可常使

由 (4),(5)和 (6)知,當(dāng)α確定后,要求 (5)成立,實(shí)際上是確定γ =α+h使 (5)成立 (參見 1.3).故有下述定義:如果不等式 (5)成立,則稱在子區(qū)間 [α,γ]上以精度為ε1的 Simpson試驗(yàn),簡稱 Simpson試驗(yàn)是成功的,否則稱為是失敗的;如果是前者,我們還稱 [α,γ]是實(shí)現(xiàn) Romberg算法的一個區(qū)間,簡稱為 Romberg區(qū)間,有理由認(rèn)為 f(x)在該區(qū)間上每一點(diǎn)的變化激烈度是相近的 (當(dāng)然,所謂‘相近’是與ε0和相應(yīng)的σ有關(guān)的).易知,在 f(x)的變化激烈度相近的區(qū)間上,Romberg算法是較為有效的.現(xiàn)設(shè)已經(jīng)對區(qū)間[a,b]完成所有的 Simpson試驗(yàn),即實(shí)現(xiàn)了對積分區(qū)間[a,b]的剖分 (剖分過程詳見 1.3),得到 N個Romberg區(qū)間 [ti,ti+1],ti＜ ti+1,i=0:N-1,t0=a,tN=b.易知,諸 ti+1-ti一般并不相等;ε0愈小,諸ti+1-ti亦愈小;而且當(dāng)σ ＞1愈大時,相應(yīng)的 ti+1-ti愈小,記 Ii是 f(x)在 [ti,ti+1]上的積分,此處,Si是用 Romberg算法近似計(jì)算 Ii所得到的數(shù)值積分,滿足

顯然有

如果γ-α =(b-a)/2k,則稱 [α,γ]是 [a,b]的一個 k級子區(qū)間,簡稱為 k級子區(qū)間;顯然,[α,β]是[a,b]的 k+1級子區(qū)間,又稱 [α,β]和 [β,γ]依次是 [α,γ]上的左 k+1級子區(qū)間和右 k+1級子區(qū)間,簡稱 [α,γ]上的左、右 k+1級子區(qū)間.

1.3 算法

1)令α =a,γ =b則γ-α =b-a,稱 [α,γ]是 [a,b]的 0級子區(qū)間,簡稱為 0級子區(qū)間.如果 [α,γ]上 Simpson試驗(yàn)成功,則 [α,γ]就是 Romberg區(qū)間,否則轉(zhuǎn) 2),

2)將 [a,b]分成 2個區(qū)間 [α,β]和 [β,γ],β -α =γ -β =(b-a)/2,則 [α,β]和 [β,γ]依次是 0級子區(qū)間 [a,b]的左、右 2個一級子區(qū)間.一般而言,設(shè) [α,γ]是 [a,b]的 k級子區(qū)間 (它是 [a,b]的 k-1級子區(qū)間 (記為 [α,δ])的左 k級子區(qū)間,此處δ=α+(b-a)/2(k-1)),并設(shè) [α,γ]上的 Simpson試驗(yàn)是失敗的.現(xiàn)按“先左后右”的原則,在 [α,γ]上的左、右 2個 k+1級子區(qū)間 [α,β]和 [β,γ]上做 Simpson試驗(yàn),共有 3種可能:

(i)如果分別在 [α,β]和 [β,γ]上的 Simpson試驗(yàn)均成功,則依次得到 2個 Romberg區(qū)間 [α,β]和[β,γ];下一步轉(zhuǎn)向 k-1級子區(qū)間 [α,δ]上的右 k級子區(qū)間 [γ,δ]上的 Simpson試驗(yàn),此處δ-γ =γ-α =(b-a)/2k.

(ii) 如果 [α,β]上的 Simpson試驗(yàn)成功 ,但在 [β,γ]上的試驗(yàn)失敗 ,則 [α,β]是 Romberg區(qū)間 ,[β,γ]不是;下一步需做 [β,γ]上的 2個 k+2級子區(qū)間 [β,ξ]和 [ξ,γ]的 Simpson試驗(yàn) ,這里ξ =(β +γ)/2.

(iii)如果 [α,β]上的 Simpson試驗(yàn)失敗,則下一步要考察 [α,β]上的 2個 k+2級子區(qū)間 [α,η]和[η,β]的 Simpson試驗(yàn),這里η =(α +β)/2.

現(xiàn)設(shè)[a,b]上的 Simpson試驗(yàn)已全部完成,這表明已用自適應(yīng) Simpson算法對區(qū)間[a,b]完成了剖分.可以看出,不管 f(x)在[a,b]上的變化性態(tài)如何,這種剖分總可完成.

3)從左至右,對每個Romberg子區(qū)間均按用戶要求的控制精度ε使用Romberg算法,并完成所有子區(qū)間上的 Romberg結(jié)果求和,即得到 (1)的近似積分.

實(shí)際計(jì)算時,只要得到 1個 Romberg區(qū)間 (一定是從左到右依次得到的),就在該子區(qū)間上使用 Romberg算法已得到該子區(qū)間上的滿足控制精度ε的 Romberg積分近似值,并與前面得到的近似結(jié)果累加,所得到的結(jié)果即為 (1)式的積分值.

說明在計(jì)算過程中,凡在 Simpson試驗(yàn)中所計(jì)算出的函數(shù)值均需要存儲,在Romberg計(jì)算過程中需要時不再重新計(jì)算.

2 數(shù)值實(shí)驗(yàn)

下面所述的自適應(yīng) S-R外推算法相對于某一算法L的節(jié)省率μ(簡稱為節(jié)省率)是指:(算法L的函數(shù)計(jì)值個數(shù) -自適應(yīng) S-R算法的函數(shù)計(jì)值個數(shù))/自適應(yīng) S-R外推算法的函數(shù)計(jì)值個數(shù).把“函數(shù)計(jì)值個數(shù)”簡記為 f數(shù),并記 H為最大積分步長,h為最小積分步長.

例 1計(jì)算積分[6]

分別取ε =1.0×10-4,ε =1.0×10-7(精確解為:I=-1.426 024 756 346 27).利用自適應(yīng) Simpson算法、Romberg外推算法與自適應(yīng) S-R算法分別計(jì)算積分,計(jì)算結(jié)果見表1,Romberg區(qū)間端點(diǎn)分布 (計(jì)算函數(shù)點(diǎn)每隔 3點(diǎn)取 1點(diǎn))見圖 1.

表1 在 2種控制精度下的計(jì)算結(jié)果

例 2計(jì)算積分

控制精度為 =1.0×10 , =1.0×10 精確解為:I=29.326 213 804 391 15.分別利用自適應(yīng)Simpson算法、Romberg算法與自適應(yīng) S-R算法計(jì)算積分,計(jì)算結(jié)果如表2,計(jì)算函數(shù)點(diǎn)的Romberg區(qū)間端點(diǎn)分布見圖2.

表2 在 2種控制精度下的計(jì)算結(jié)果

例 3計(jì)算積分,控制精度分別為ε=1.0×10-5,ε=1.0×10-8(I＊=0.959 572 318 005)分別利用自適應(yīng) Simpson算法、Romberg算法與自適應(yīng) S-R算法計(jì)算積分,計(jì)算結(jié)果見表3.

表3 在 2種控制精度下的計(jì)算結(jié)果

3 結(jié)論

大量數(shù)值實(shí)驗(yàn)表明,針對數(shù)值求積,本文構(gòu)造的既能變步長計(jì)算又能逐步分段提高收斂階的自適應(yīng)S-R算法,與自適應(yīng) Simpson算法和 Romberg外推算法相比較,當(dāng)被積函數(shù)在積分區(qū)間上變化性態(tài)急劇多變時,具有明顯的效率優(yōu)勢,而且,變化性態(tài)愈急劇多變,效率優(yōu)勢愈明顯,而在相反的情形,其效率也不低于Romberg外推算法.而且可以看出積分精度要求愈高,節(jié)省率愈高或可以獲得更高的精度.

[1]LYNESS J N.Notes on the adaptive Simpson quad-rature routine[J].J Assoc Comput,1969,16:83-495.

[2]GANDER W,GAUTSCH IW.Adaptive quadrature-revisited[J].BitNumericalMathematics,2000,40(1):84-101.

[3]ESPEL ID TO.Doubly adaptive quadrature routines based on Newton-Cotes rules[J].BitNumericalMathematics,2003,43(2):319-337.

[4]JOYCE D C.Survey of extrapolation processes in numerical analysis[J].S IAM Review,1971,13(4):435-490.

[5]楊錄峰,金云超.一類含奇點(diǎn)函數(shù)的數(shù)值積分方法[J].云南民族大學(xué)學(xué)報(bào):自然科學(xué)版,2010,19(1):16-19.

[6]R ICHARD L,BURDEN J.DOUGLAS F.Numerical analysis[M].7th ed.北京:高等教育出版社,2001.

[7]MESTROV ICM,OCV I RK E.An application of Romberg extrapolation on quadrature method for solving linearVolterra integral equations of the second kind[J].AppliedMathematics and Computation,2007,194(2):389-393.

[8]劉長虹,關(guān)永亮,壽卓佳,等.蒙特卡洛法在數(shù)值積分上的應(yīng)用[J].上海工程技術(shù)大學(xué)學(xué)報(bào),2010,24(1):43-46.

[9]石東洋,石東偉.一個非常規(guī)高效的數(shù)值積分方法[J].河南師范大學(xué)學(xué)報(bào):自然科學(xué)版,2008,36(1):6-8.

(責(zé)任編輯萬志瓊)

An Adaptive Numerical Quadrature Algorithem of Comput ing with Both Changeable Step and Changeable Order

YANG Lu-feng1,MA Ning2,ZHAO Shuang-suo1,3
(1.Institute of Info rmation and Computation Science,BeifangUniversity ofNationalties,Yinchan 750021,China;2.WuzhongMeteorologicalBureau,Wuzhong 751300,China;3.Institute ofMathematics and Statistics,Lanzhou University,Lanzhou 730001,China)

Adaptive S-R(s impson-romberg)Algorithm,a new type of adaptive numerical quadrature algorithm,is presented,which combines the adaptive S impson algorithm with Romberg extrapolation algorithm.It possesses both the merits of the step-changeable computation and those of the gradual convergence order of numerical quadrature.The comparison of several numerical examples has indicated that the adaptive S-R algorithm has obvious superiority compared with the adaptive Simpson algorithm and the Romberg extrapolation algorithm when the integrand has some irregular property on the integration interval and the variation property over the interval of integration changes significantly.

numerical integration;adaptive simpson algorithm;romberg extrapolation algorithm;adaptive S-Ralgorithm

O 241.4

A

1672-8513(2011)01-0032-05

10.3969/j.issn.1672-8513.2011.01.008

2010-09-29.

國家自然科學(xué)基金(1096100);國家民委科研項(xiàng)目(09BF07);北方民族大學(xué)科研項(xiàng)目(2008Y032).

楊錄峰 (1980-),男,碩士,助教.主要研究方向:數(shù)值計(jì)算、計(jì)算流體力學(xué).