蔣冰峰，熊小勇，胡艷霞

(湖北民族學(xué)院 理學(xué)院， 湖北 恩施 445000)

一類振蕩無(wú)窮積分的一種數(shù)值算法

蔣冰峰，熊小勇，胡艷霞

(湖北民族學(xué)院 理學(xué)院， 湖北 恩施 445000)

振蕩函數(shù)無(wú)窮范圍積分的數(shù)值計(jì)算方法是先求出振蕩函數(shù)的各級(jí)零點(diǎn)，根據(jù)被積函數(shù)在相鄰零點(diǎn)間積分，再把積分結(jié)果求和.但是，求和的結(jié)果收斂性質(zhì)非常差.我們采用W外推算法，來(lái)加快求和結(jié)果的收斂性.我們以一個(gè)簡(jiǎn)單的例子說(shuō)明W外推算法的有效性.

振蕩函數(shù)；無(wú)窮積分；算法

在科學(xué)計(jì)算及工程技術(shù)應(yīng)用中，經(jīng)常會(huì)涉及到計(jì)算振蕩函數(shù)的無(wú)窮積分，形如[1-3]:

(1)

形如式(1)的二維數(shù)值積分非常復(fù)雜，為了簡(jiǎn)便性，也不失一般性，先討論以下一維振蕩函數(shù)f(x)的無(wú)窮積分(a≥0)：

(2)

在處理積分或數(shù)列的外推方法中，常用的有Euler 變換，ε算法和W變換[4].其中由Sidi建立的W變換[5-7]是非常有效的外推算法[8].在本文中，根據(jù)W變換編寫(xiě)外推算法程序，并以一個(gè)簡(jiǎn)單的積分為例，介紹如何應(yīng)用外推算法計(jì)算振蕩積分的無(wú)窮積分.

1 W變換

根據(jù)振蕩函數(shù)f(x)的一系列零點(diǎn)x1,x2,x3,…，構(gòu)造如下函數(shù):

(3)

(4)

(5)

(6)

(7)

其中s,p=0,1,2，….

2 ISE方法求積分流程

以下將用實(shí)例來(lái)說(shuō)明如何使用ISE方法求解振蕩函數(shù)的無(wú)窮積分.

(8)

第一步：由例子可知x0=0，然后根據(jù)被積函數(shù)求其零點(diǎn)x1,x2,x3,…本例子中，振蕩函數(shù)為三角函數(shù)，其零點(diǎn)非常容易求得.式(1)中的振蕩函數(shù)J0的零點(diǎn)的求解要復(fù)雜得多，但是根據(jù)Temme的算法[9-10]，CERNLIB已有現(xiàn)成的算法Dbzejy() 計(jì)算各種類型貝塞爾函數(shù)的零點(diǎn)[11].

第二步：求被積函數(shù)sin(πx2/2)在所有零點(diǎn)間隔[xi,xi+1]的ψ(xi)以及G(xi).出于計(jì)算的精確性，此處積分程序可采用Piessens等人[12]自適應(yīng)算法基礎(chǔ)上的子程序Dqag()或Dqng()[13-14].

從表1中可以看出，即使是積分上限取到被積函數(shù)的第1000個(gè)零點(diǎn)，計(jì)算結(jié)果0.507114061665686同標(biāo)準(zhǔn)結(jié)果相比，差距仍然是非常明顯，不可忽略的.第二，計(jì)算結(jié)果在標(biāo)準(zhǔn)值0.5附近振蕩，這也反應(yīng)了被積函數(shù)的振蕩性質(zhì).

第三步：根據(jù)ψ(xi)以及G(xi)和W變換算法，編寫(xiě)外推算法子程序.根據(jù)ψ(xi)以及G(xi)和外推程序，得到最終結(jié)果.

只需要式(8)中前10個(gè)零點(diǎn)間隔的積分結(jié)果，采用W-transformation外推算法(4)(5)(6)(7)，就可以得到精度非常高的結(jié)果，見(jiàn)表2.

表1無(wú)外推算法時(shí)991-1000個(gè)零點(diǎn)積分和結(jié)果

Tab.1 Results for partial sums for 991 to 1000 zeros of the integrand without extrapolation

NResult9910．4928537258154769920．5070426610435239930．4926609188370599940．5071354788738719950．4928680901820369960．5071283183264049970．4928752399590619980．5071211792928509990．49288239827602810000．507114061665686

表2 W-transformation所得結(jié)果

Tab.2 Results for the first ten zeros of the integrand with W-extrapolation

NResult10．43771710100162720．51612305025054130．49616345317763840．50085599440944350．49981791710354860．50003736201706270．49999253339644780．50000144939663890．499999713844729100．500000045763800

PROGRAM MAIN

double precision a,abserr,b,f,epsabs,epsrel,result,c1,c2,pi,

& psi,hi,i,x,m,n,w,resul

integer ier,neval,lst,limst,lsta,limsta,p,s,n1,l,j,pl,jl,i1,i2,i3

parameter(limst=50)

dimension hi(-1:limst),psi(-1:limst), x(-1:limst),

& m(-1:limst,0:limst)

& ,n(-1:limst,0:limst), w(-1:limst,0:limst)

external f

data pi /3.1415926535 8979323846 2643383279 50 d0/

epsabs=0.0d0

epsrel=1.0d-6

L=Limst-1

x(-1)=0.0d0

**********************************************

****** 算法程序求被積函數(shù)的前50個(gè)零點(diǎn)

**********************************************

do 10 i1=0, limst

x(i1)=sqrt(2.0d0*(i1+1))

10 continue

IF(lst.le.(-1)) hi(lst)=0.0D0

IF(lst.lt.(-1)) psi(lst)=0.0D0

**********************************************

****** 在相鄰零點(diǎn)間被積函數(shù)的積分及結(jié)果的和

**********************************************

do 20 lst=-1,l

c1=x(lst)

c2=x(lst+1)

call dqng(f,c1,c2,epsabs,epsrel,result,abserr,neval,ier)

Psi(lst)=result

hi(lst+1)=hi(lst)+psi(lst)

write(*,*)lst, psi(lst), hi(lst+1)

20continue

**********************************************

****** 應(yīng)用外推算法，得到最終結(jié)果resul

**********************************************

call mWextra(l,psi,hi,x,m,n,w,resul,epsrel)

write(*,*)resul

stop

end

**********************************************

****** 根據(jù)W算法編寫(xiě)的外推算法子程序

**********************************************

subroutine mWextra(L,psi,hi,x,m,n,w,resul,epsrel)

integer l,p,s,j,jl

double precision psi,hi,x,m,n,w,resul,epsrel

dimension psi(-1:l+1),hi(-1:l+1),x(-1:l+1),m(-1:l+1,0:l+1),

& n(-1:l+1,0:l+1),w(-1:l+1,0:l+1)

do 30 s=0,L

m(-1,s)=hi(s)/psi(s)

n(-1,s)=1.0d0/psi(s)

m(0,s)=(m(-1,s)-m(-1,s+1))/(x(s)**(-1)-x(s+1)**(-1))

n(0,s)=(n(-1,s)-n(-1,s+1))/(x(s)**(-1)-x(s+1)**(-1))

write(*,*)s,m(-1,s),n(-1,s)

write(*,*)s,m(0,s),n(0,s)

30 continue

do 40 p=1,l

jl=l-p

do 50 j=0,jl

m(p,j)=(m(p-1,j)-m(p-1,j+1))/(x(j)**(-1)-x(j+p+1)**(-1))

n(p,j)=(n(p-1,j)-n(p-1,j+1))/(x(j)**(-1)-x(j+p+1)**(-1))

50 continue

w(p,0)=m(p,0)/n(p,0)

write(*,*) p, w(p,0)

if(abs(w(p,0)-w(p-1,0)).le.epsrel) then

resul=w(p,0)

write(*,*)"The final result"

write(*,*) p, resul

go to 60

end if

40 continue

if(abs(w(l,0)-w(l-1,0)).gt.epsrel) then

write(*,*)"please enhance the L"

end if

60 return

end

**********************************************

****** 定義被積函數(shù)

**********************************************

double precision function f(x)

double precision x

data pi /3.1415926535 8979323846 2643383279 50 d0/

f=sin(pi/2*x**2)

return

end

[1] Jiang B F,Li J R.The polarized charge distribution induced by a fast parton in the viscous quark-gluon plasma[J].J Phys G:Nucl Part Phys,2012,39:1-10.

[2] Jiang B F,Li J R.The wake potential in the viscous quark-gluon plasma[J]. Nucl Phys A, 2011, 856:121-133.

[3] Jiang B F,Li J R Non-abelian effects on wake potential in quark-gluon plasma[J]. Nucl Phys A, 2010,832:100-111.

[4] Davis P J,Rabinowitz P. Methods of Numerical Interaction[M]. London: Acadamic Press,1984.

[5] Sidi A. The numerical evaluation of very oscillatory infinite integrals by extrapolation[J].Math Comput,1982,38:517-529.

[6] Sidi A. A user-friendly extrapolation method for oscillatory infinite integrals[J]. Math Comput,1988,51:249-266.

[7] Sidi A. Practical extrapolation methods[M]. Cambridge: Cambridge Univ Press,2003.

[8] Lucas S K, Stone H A.Evaluating infinite integrals involving Bessel functions of arbitrary order[J].J Comput Appl Math, 1995, 64:217-231.

[9] Temme N M. On the numerical evaluation of the ordinary Bessel function of the second kind[J].J Comput Phys, 1976, 21:343-350.

[10] Temme N M.An algorithm with Algol 60 program for the computation of the zeros of ordinary Bessel functions and those of their derivatives[J].J Comput Phys,1979,32:270-279.

[11] Piessens R， Doncker-Kapenga E De， Uberhuber C W,et al. A subroutine package for automaticintergration[M].Berlin：Spinger，1983.

TheNumericalAlgorithmofInfiniteIntegralsInvolvingaKindofOscillatoryFunction

JIANG Bing-feng,XIONG Xiao-yong,HU Yan-xia

(School of Science,Hubei University for Nationalities,Enshi 445000,China)

The method of evaluating infinite integrals involving oscillatory function is divided into three steps: finding zeros of the oscillatory function, doing integral at its zeros and forming a sequence of partial sums, using extrapolation to accelerate convergence to obtain the final result. In this paper, W-transformation is used to accelerate the convergence of the infinite integral of oscillatory function.

Oscillatory function; infinite integral;aglorithm

2012-11-13.

國(guó)家自然科學(xué)基金項(xiàng)目(11147012)；湖北民族學(xué)院博士啟動(dòng)基金項(xiàng)目(MY2011B002).

蔣冰峰(1977- ),男，博士，講師，主要從事夸克物質(zhì)理論研究.

O241.4

A

1008-8423(2013)02-0163-04