紀(jì)宇　何一璇　吳國群　吳敏

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

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

中圖分類號(hào)：TP311.5

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

DOI： 10.3969/j.issn.1000-5641.2019.06.006

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

（1）J0（y;b）

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

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

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

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

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

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

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

（2）J1（y;b）

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

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

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

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

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

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

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

5.2

J2 （y;b）

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

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

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

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

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

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

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

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

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

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

6 結(jié)語

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

[參考文獻(xiàn)]

[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] .上海：華東師范大學(xué)， 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.

（責(zé)任編輯：林磊）

收稿日期：2018-11-30

基金項(xiàng)目：國家自然科學(xué)基金（11371143）

第一作者：紀(jì)宇，女，碩士研究生，研究方向?yàn)閿?shù)值計(jì)算、符號(hào)計(jì)算.E-mail： 15526476747@163.com.

通信作者：吳敏，女，副教授，主要研究方向?yàn)榉?hào)計(jì)算、數(shù)值計(jì)算.E-mail： mwu@sei.ecnu.edu.cn.