紀宇　何一璇　吳國群　吳敏

摘要：貝塞爾函數(shù)的數(shù)值逼近既有重要的理論意義，又在數(shù)學、物理學、工程等各個領域有著廣泛的應用.研究整數(shù)階第一類貝塞爾函數(shù)的數(shù)值逼近，基于Prony方法，采用不同三角函數(shù)（正弦、余弦）形式的Prony-like方法進行逼近.通過在符號計算軟件Maple中對函數(shù)進行數(shù)值實驗，分析不同整數(shù)階的第一類貝塞爾函數(shù)在不同自變量區(qū)間上的數(shù)值逼近，將Prony-like方法的實驗結(jié)果與基于傅里葉級數(shù)展開的方法進行對比，發(fā)現(xiàn)Prony-like方法的逼近效果遠優(yōu)于基于傅里葉級數(shù)的方法.

關鍵詞：貝塞爾函數(shù);Prony-like方法;三角函數(shù);基于傅里葉級數(shù)展開;數(shù)值逼近

中圖分類號：TP311.5

文獻標志碼：A

DOI： 10.3969/j.issn.1000-5641.2019.06.006

0 引言

特殊函數(shù)是指一些具有特定性質(zhì)的函數(shù)，它們廣泛應用于物理學、化學、工程、統(tǒng)計學以及計算機科學等領域，大多數(shù)特殊函數(shù)都沒有閉形式的表達式，其賦值只能通過數(shù)值逼近來實現(xiàn).由于特殊函數(shù)在應用中的重要性和廣泛性，因此，如何提高數(shù)值逼近的效率，具有重要的學術意義和應用價值.

有關特殊函數(shù)的性質(zhì)和數(shù)值計算方法，眾多國內(nèi)外學者作出了大量的研究.1964年，Abramowitz、Stegun和Romain在文獻[1]中詳細介紹了特殊函數(shù)的定義、性質(zhì)、數(shù)值逼近方法.2010年，美國國際標準與技術研究院（NIST）發(fā)布了NIST數(shù)學函數(shù)數(shù)字圖書館（Digital Library of Mathematical Functions，DLMF）[2]，系統(tǒng)地介紹了特殊函數(shù)的性質(zhì)和常用賦值方法.

貝塞爾函數(shù)（Bessel functions）是德國數(shù)學家貝塞爾于1817年研究三體引力系統(tǒng)的運動問題時提出的，貝塞爾函數(shù)在波的傳播、有勢場和信號處理領域有著廣泛的應用，如圓柱形波導中的電磁波傳播、調(diào)頻合成等.

貝塞爾函數(shù)是一類特殊函數(shù)，無法由初等函數(shù)直接表示.針對貝塞爾函數(shù)的數(shù)值逼近，已有的研究提出了許多數(shù)值方法[1-2].文獻[3]對整數(shù)階和半整數(shù)階第一類貝塞爾函數(shù)的冪級數(shù)、漸近級數(shù)和連分式展開進行了比較.文獻[4]對任意整數(shù)階的貝塞爾函數(shù)討論了多項式逼近.文獻[5]提出了基于泰勒級數(shù)展開、Debye漸近展開和Sommerfeld積分的算法，對貝塞爾函數(shù)進行逼近.

在之前的工作中，冪級數(shù)、漸近級數(shù)展開等數(shù)值方法對整數(shù)階第一類貝塞爾函數(shù)的逼近效率不高，且在數(shù)值上不穩(wěn)定.我們注意到貝塞爾函數(shù)表現(xiàn)出近似周期性的性質(zhì)，因此，在本文中，針對貝塞爾函數(shù)的這一性質(zhì)，考慮通過三角函數(shù)的和形式（∑墨，ai cos（cpix）或∑n：lai sin（cpix））對貝塞爾函數(shù)進行逼近.

基于傅里葉級數(shù)的方法是基于這種思路的一種應用.傅里葉級數(shù)在理論、數(shù)學和工程領域都有著重要的研究價值.傅里葉級數(shù)方法是把類似波的函數(shù)表示成簡單正弦波的和，它把任何周期函數(shù)或周期信號分解成簡單震蕩函數(shù)（如正弦函數(shù)和余弦函數(shù)）的集合.在文獻[6]中，提出了利用傅里葉級數(shù)展開對貝塞爾函數(shù)進行逼近，然而，正如本文實驗所示，形式為ai cos（kx）和ai sin（kx）的和的逼近效果并不理想.

我們轉(zhuǎn)而考慮另一種三角函數(shù)形式的逼近方法.Prony方法作為傅里葉級數(shù)的一種拓展，最早由Prony在1795年提出[7].Prony方法通過在一個均勻采樣的樣本中提取有用的信息，并建立一系列的衰減指數(shù)或正弦函數(shù).Prony方法常作為一種線譜估計方法應用于各領域的信號處理.但Prony模型在數(shù)值上是病態(tài)的，非最優(yōu)、非嚴格的近似求解算法，對噪音的干擾非常敏感，這極大地限制了Prony模型的應用.Prony方法的復興起源于對Prony方法的修改，這顯著穩(wěn)定了原始方法，例如，ESPRIT方法，矩陣束方法或近似Prony方法[8-11].近年來，Prony方法也已經(jīng)成功應用于不同的逆問題，如，量子化學[12]或流體動力學[13]中Green函數(shù)的近似，反向散射中粒子的定位[14]等.

本文中，我們針對整數(shù)階第一類貝塞爾函數(shù)，基于擴展形式的Prony方法，采用三角函數(shù)和的形式對第一類貝塞爾函數(shù)進行逼近，并將實驗結(jié)果與基于傅里葉級數(shù)方法的實驗結(jié)果作比較.

本文第1節(jié)介紹了第一類貝塞爾函數(shù)的定義和相關性質(zhì);第2節(jié)介紹了貝塞爾函數(shù)的基于傅里葉級數(shù)展開形式;第3節(jié)和第4節(jié)分別介紹了指數(shù)形式的Prony方法和三角函數(shù)形式的Prony-like方法;第5節(jié)分別應用基于傅里葉級數(shù)方法和Prony-like方法對第一類貝塞爾函數(shù)進行數(shù)值逼近，并比較兩種方法的逼近效果，

本文中，R、Z、C分別表示實數(shù)集、整數(shù)集和復數(shù)集.

1 第一類貝塞爾函數(shù)

然而，上述Prony方法中的多項式求根問題在數(shù)值上是病態(tài)的，非最優(yōu)、非嚴格的近似求解算法對噪音的干擾非常敏感，這極大地限制了Prony方法的應用.1999年，Golub等提出Prony方法的一種重要表述，將求解Hankel系統(tǒng)和求解生成多項式的根的問題轉(zhuǎn)化為求解一個廣義特征值問題[15].廣義特征值問題一般可通過Cholesky分解和LU分解來有效求解，存在成熟的且數(shù)值穩(wěn)定的數(shù)值方法，因此，這一新表示使Prony方法的應用得以復興.

下面，我們來簡單介紹一下基于廣義特征值問題表述的Prony方法.

給定正實數(shù)b和區(qū)間|a，b|上的實函數(shù)f（x），計算如下形式的展開式：

Prony方法可以把求解ai、φi的非線性問題轉(zhuǎn)化為求解兩個線性方程組的問題，即廣義特征值問題和含范德蒙矩陣的線性方程組問題.

5 基于傅里葉級數(shù)方法與Prony-like方法的實驗結(jié)果對比

在本節(jié)中，我們用基于傅里葉級數(shù)的方法和Prony-like方法對J（y;b）進行函數(shù)逼近，并比較其逼近結(jié)果.為使計算精度盡可能高，所有實驗均在計算機代數(shù)軟件Maple中進行.我們通過設置Digits，來保證實驗結(jié)果的精度.

5.1 J0（y;b）和J1（y;b）

由J0（Y）和J1（y）的奇偶性以及知，JO （y;b）和J1（y;b）是偶函數(shù).所以應用4.1節(jié)中的cosine形式的Prony-like方法在區(qū)間[O，b]上對Jo（Y）、J1（y）進行數(shù)值逼近，并與基于傅里葉級數(shù)的方法進行比較.為比較兩種方法的逼近效果，我們在區(qū)間[0，b]上按照1/2 10的間隔取樣本點，并計算出Prony-like方法和基于傅里葉級數(shù)方法在這些樣本點的相對誤差.為簡便起見，我們?nèi)∠鄬φ`差的對數(shù)值（log10）.

（1）J0（y;b）

在n=0，b=1時，分別用基于傅里葉級數(shù)方法和Prony-like方法對J0（y;1）進行了數(shù)值逼近.

表1展示了JO （y;1）在6=l，即自變量區(qū)間為[O，1l時，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，10，25，50，75，100時的最大相對誤差（l09io）.

圖3分別展示了Jo（Y;1）在自變量區(qū)間[0，1]上，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，25，50，100時的相對誤差（log10），其中橫坐標表示自變量x的值，縱坐標表示相對誤差對數(shù)值，紅色曲線和黑色曲線分別代表基于傅里葉級數(shù)方法和Prony-like方法的相對誤差.

設6為正實數(shù)，我們分別用基于傅里葉級數(shù)的方法和Prony-like方法對JO （y;b）進行了數(shù)值逼近.

表2和表3分別展示了JO （y;b）在b=5和b=20，即在區(qū)間[0，5]和[0，20]上，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，10，25，50，75，100時的最大相對誤差（log10）.

圖4展示了JO （y;b）在區(qū)間為[0，5]時，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，25，50，100時的相對誤差（log10）.其中，橫坐標表示自變量x的值，縱坐標表示相對誤差對數(shù)值，紅色曲線和黑色曲線分別代表基于傅里葉級數(shù)方法和Prony-like方法的相對誤差.

從表1、表2、表3及圖3、圖4中可以看出，Prony-like方法的最大相對誤差遠小于基于傅里葉級數(shù)的方法，且相對誤差隨展開項數(shù)的增大而效果更好.結(jié)合JO （y;1）的圖象，發(fā)現(xiàn)相對誤差隨取值范圍[0，b]增大而增大，且相對誤差隨展開項數(shù)的增大而變好.若想獲得與JO （y;1）同樣的逼近效果，需要更大的展開項數(shù).

（2）J1（y;b）

在n=l，b=l時，分別用基于傅里葉級數(shù)方法和Prony-like方法對J1（y;1）進行數(shù)值逼近.

表4展示了Jl（y;1）在6=1，即自變量區(qū)間為[0，1]時，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，10，25，50，75，100時的最大相對誤差（l09io）.

圖5展示了J1（y;1）在自變量區(qū)間為[0，1]時，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，25，50，100時的相對誤差（log10）.其中，橫坐標表示自變量x的值，縱坐標表示相對誤差對數(shù)值，紅色曲線和黑色曲線分別代表基于傅里葉級數(shù)方法和Prony-like方法的相對誤差.

在n=l，b>l時，我們分別用基于傅里葉級數(shù)方法和Pronv-like方法對J1（y;b）進行了數(shù)值逼近，

表5、表6展示了J1（y;b）在b=5和b=20，即區(qū)間為[0，5]和[0，20]時，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，10，25，50，75，100時的最大相對誤差（log10）.

圖6給出了在區(qū)間[0，5]上，J1 （y;b）基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，25，50，100時的相對誤差（log10）.其中，橫坐標表示自變量x的值，縱坐標表示相對誤差對數(shù)值，紅色曲線和黑色曲線分別代表基于傅里葉級數(shù)方法和Prony-like方法的相對誤差.

從表4、表5、表6及圖5、圖6中可以看出，類似于JO （y;b），J1（y;b）的Prony-like方法的最大相對誤差遠小于基于傅里葉級數(shù)方法，且相對誤差隨展開項數(shù)的增大而效果更好.結(jié)合J1（y;1）的圖象，發(fā)現(xiàn)相對誤差隨展開項數(shù)的增大而變好，但同時相對誤差隨取值范圍[0，b]的增大而增大.因而，為了獲得與J1（y;1）同樣的逼近效果，需要展開更多的項數(shù).

5.2

J2 （y;b）

由于J2（x）是偶函數(shù)，因而J2（y;b）=J2（y）/ y/b是奇函數(shù)，我們應用基于傅里葉級數(shù)方法和4.2節(jié)介紹的sine形式的Prony-like方法在區(qū)間[0，b]上對J2 （y;b）進行數(shù)值逼近.

令b= 1.我們分別用基于傅里葉級數(shù)方法和Prony-like方法對J2 （y;1）進行數(shù)值逼近.

表7展示了J2 （y;1）在自變量區(qū)間為[0，1]時，基于傅里葉級數(shù)方法和Prony-like方法在展

圖7分別展示了J2（y;1）在自變量區(qū)間[O，1]上，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，25，50，100時的相對誤差（log10）.其中，橫坐標表示自變量x的值，縱坐標表示相對誤差對數(shù)值，紅色曲線和黑色曲線分別代表基于傅里葉級數(shù)方法和Prony-like方法的相對誤差.

以下，設b>l，分別應用基于傅里葉級數(shù)方法和Prony-like方法對J2 （y;b）進行了數(shù)值逼近，

表8和表9展示了J2 （y;b）在b=5和b=20，即區(qū)間為[0）5]和[O，20]時，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，10，25，50，75，100時的最大相對誤差（log10）.

圖8展示了J2 （y;b）在區(qū)間[0，5]上，基于傅里葉級數(shù)方法和Prony-like方法在展開項數(shù)m=5，25，50，100時的相對誤差（log10）.其中，橫坐標表示自變量z的值，縱坐標表示相對誤差對數(shù)值，紅色曲線和黑色曲線分別代表基于傅里葉級數(shù)方法和Prony-like方法的相對誤差.

從表7、表8、表9及圖7、圖8中可以看出，J2 （y;b）與JO （y;b）、J1（y;b）的實驗結(jié)果相似，即Prony-like方法的逼近效果遠優(yōu)于基于傅里葉級數(shù)的方法.5.3基于傅里葉級數(shù)方法與Prony-like方法計算時間對比

本小節(jié)以函數(shù)JO （y;b）、J1（y;b）和J2 （y;b）的逼近為例，比較基于傅里葉級數(shù)的方法和Prony-like方法在不同展開項數(shù)和不同區(qū)間下的計算時間.表10給出了比較結(jié)果.

從表10可以看出，基于傅里葉級數(shù)方法的計算時間比Prony-like方法的計算時間短，旦隨著展開項數(shù)的增加，兩種方法計算時間的差距逐漸拉大，兩種方法受自變量區(qū)間和階數(shù)影響均不大.然而，Prony-like方法可以在項數(shù)很小時得到非常高的精度，如J0 （y;1）在展開5項時，Prony-like方法比基于傅里葉級數(shù)方法的精度高約10，計算時間僅長0.154 s，所以Prony-like方法逼近效果更好.

6 結(jié)語

本文對整數(shù)階第一類貝塞爾函數(shù)進行數(shù)值逼近，基于前面的分析和實驗可得，cosme和slne形式Prony-like方法的數(shù)值逼近結(jié)果遠遠優(yōu)于基于傅里葉級數(shù)方法，雖然Prony-like方法的計算時間略長于基于傅里葉級數(shù)方法的計算時間，但是Prony-like方法可以在項數(shù)很小時就得到很高的精度，所以Prony-like方法的逼近效果更好.在未來的工作中，將針對復數(shù)階的貝塞爾函數(shù)應用Prony-like方法進行實驗.作為Prony-like方法的改進，下一步的研究將側(cè)重于改進通過 Hankel矩陣和廣義特征值問J題求解φi的復雜計算過程，進一步提高 -角函數(shù)和形式的計算

[參考文獻]

[1]ABRAMOWITZ M， STEGUN I A， ROMAIN J E. Handbook of Mathematical Functions with Formulas， Graphs：nd Mathematical Tables [M]. New York： Dover， 1964.

[2] LOZIER D W. NIST Digital library of mathematical functions [J]. Annals of Mathematics & Artificial Intelli-gence： 2003， 38（1/3）： 105-119.

[3] 常曉陽. JL類特V殊函數(shù)的賦值分析研究 [D] .上海：華東師范大學， 2018.

[4]MILLANE R P， EADS J L. Polynomial approximations to Bessel functions [J]. IEEE Transactions on Antennas& Propagation， 2003， 51（1）： 1398-1400.

[5] MATVIYENKO G. On the evaluation of Bessel functions [J]. Applied and Computational Harmonic Analysis，1993， 1（1）： 116-135.

[6]ANDRUSYK A. Infinite series representations for Bessel functions of the first kind of integer order [EB/OLl.[2018-10-20]. https：//arxiv.org/abs/1210.2109

[7]HILDEBRAND F B. Introduction to Numerical Analysis [M] . New Delhi： Tata McGraw-Hill Publishing Company Limited， 1956.

[8]B013MANN F， PLONKA G， PETER T， et al. Sparse deconvolution methods for ultrasonic NDT [J]. Journal ofNondestructive Evaluation， 2012， 31（3）： 225-244.

[9]ROY R， PAULRAJ A. KAILATH T. Estimation of signal parameters via rotational invariance techniques-ESPRIT [C]// Proceedings of the SPIE. International Society for Optics and Photonics， 1986： 94-101.

[10] HUA Y： SARKAR T K. Matrix pencil method for estimating parameters of exponentially damped/undampedsinusoids in noise [J]. IEEE Transactions on Acoustics， Speech， and Signal Processing， 1990， 38（5）： 814-824.

[11] POTTS D， TASCHE M. Nonlinear approximation by sums of nonincreasing exponentials [J] . Applicable Analysis：2011， 90（3/4）： 18.

[12] YANAI T， FANN G I， GAN Z. et al. Multiresolution quantum chemistry in multiwavelet bases： Analyticderivatives for Hartree-Fock and density functional theory [J] . Journal of Chemical Physics： 2004， 121（7） ： 2866-2876.

[13] BEYLKIN G， CRAMER R， FANN G， et al. Multiresolution separated representations of singular and weaklysingular operators [J] Applied & Computational Harmonic Analysis， 2007， 23（2）： 235-253.

[14]HANKE M. One shot inverse scattering via rational approximation [J]. Siam Journal on Imaging Sciences： 2012，5（1）： 465-482.

[15] GOLUB G H， MILANFAR P， VARAH J. A Stable Numerical Method for Inverting Shape from Moments [M].Society for Industrial and Applied Mathematics： 1999.

[16] GIESBRECHT M， LABAHN G， LEE W S. Symbolic-numeric sparse polynomial interpolation in Chebyshevbasis and trigonometric interpolation [C]// Proceedings of Computer Algebra in Scientific Computing （CASC2004） ， 2004： 195206.

[17]CUYT A. LEE W S. Generalized eigenvalue relations： Spread polynomials and sine [R]. 2018.

（責任編輯：林磊）

收稿日期：2018-11-30

基金項目：國家自然科學基金（11371143）

第一作者：紀宇，女，碩士研究生，研究方向為數(shù)值計算、符號計算.E-mail： 15526476747@163.com.

通信作者：吳敏，女，副教授，主要研究方向為符號計算、數(shù)值計算.E-mail： mwu@sei.ecnu.edu.cn.