蔣冰峰，熊小勇，胡艷霞

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

一類振蕩無窮積分的一種數值算法

蔣冰峰，熊小勇，胡艷霞

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

振蕩函數無窮范圍積分的數值計算方法是先求出振蕩函數的各級零點，根據被積函數在相鄰零點間積分，再把積分結果求和.但是，求和的結果收斂性質非常差.我們采用W外推算法，來加快求和結果的收斂性.我們以一個簡單的例子說明W外推算法的有效性.

振蕩函數；無窮積分；算法

在科學計算及工程技術應用中，經常會涉及到計算振蕩函數的無窮積分，形如[1-3]:

(1)

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

(2)

在處理積分或數列的外推方法中，常用的有Euler 變換，ε算法和W變換[4].其中由Sidi建立的W變換[5-7]是非常有效的外推算法[8].在本文中，根據W變換編寫外推算法程序，并以一個簡單的積分為例，介紹如何應用外推算法計算振蕩積分的無窮積分.

1 W變換

根據振蕩函數f(x)的一系列零點x1,x2,x3,…，構造如下函數:

(3)

(4)

(5)

(6)

(7)

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

2 ISE方法求積分流程

以下將用實例來說明如何使用ISE方法求解振蕩函數的無窮積分.

(8)

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

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

從表1中可以看出，即使是積分上限取到被積函數的第1000個零點，計算結果0.507114061665686同標準結果相比，差距仍然是非常明顯，不可忽略的.第二，計算結果在標準值0.5附近振蕩，這也反應了被積函數的振蕩性質.

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

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

表1無外推算法時991-1000個零點積分和結果

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所得結果

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

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

****** 算法程序求被積函數的前50個零點

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

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

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

****** 在相鄰零點間被積函數的積分及結果的和

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

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

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

****** 應用外推算法，得到最終結果resul

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

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

write(*,*)resul

stop

end

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

****** 根據W算法編寫的外推算法子程序

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

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

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

****** 定義被積函數

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

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.

國家自然科學基金項目(11147012)；湖北民族學院博士啟動基金項目(MY2011B002).

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

O241.4

A

1008-8423(2013)02-0163-04