誤差函數(shù)Chebyshev級(jí)數(shù)的計(jì)算方法

鄧國(guó)強(qiáng)，唐 敏

(1.桂林電子科技大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣西 桂林 541004；2.桂林電子科技大學(xué) 廣西高校數(shù)據(jù)分析與計(jì)算重點(diǎn)實(shí)驗(yàn)室，廣西 桂林 541004)

誤差函數(shù)Chebyshev級(jí)數(shù)的計(jì)算方法

鄧國(guó)強(qiáng)1,2，唐 敏1,2

(1.桂林電子科技大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣西 桂林 541004；2.桂林電子科技大學(xué) 廣西高校數(shù)據(jù)分析與計(jì)算重點(diǎn)實(shí)驗(yàn)室，廣西 桂林 541004)

為了在IEEE浮點(diǎn)計(jì)算環(huán)境下對(duì)誤差函數(shù)進(jìn)行精確有效地賦值，提出了誤差函數(shù)的Chebyshev級(jí)數(shù)計(jì)算方法。采用Clenshaw算法計(jì)算級(jí)數(shù)的前N項(xiàng)部分和，減小求和的舍入誤差。實(shí)驗(yàn)結(jié)果表明，針對(duì)誤差函數(shù)的賦值問(wèn)題，Chebyshev級(jí)數(shù)比Taylor級(jí)數(shù)的收斂速度更快，即達(dá)到相同的賦值精度要求時(shí)，Chebyshev級(jí)數(shù)法需要的項(xiàng)數(shù)遠(yuǎn)少于Taylor級(jí)數(shù)法。

誤差函數(shù); Chebyshev級(jí)數(shù); Clenshaw算法; IEEE浮點(diǎn)計(jì)算標(biāo)準(zhǔn)

特殊函數(shù)是一類(lèi)重要的數(shù)學(xué)函數(shù)，在計(jì)算機(jī)科學(xué)、數(shù)學(xué)、物理、工程、化學(xué)和統(tǒng)計(jì)學(xué)上都有廣泛的應(yīng)用。通常來(lái)說(shuō)，特殊函數(shù)沒(méi)有一個(gè)形式上的定義，但必須是有廣泛應(yīng)用的、滿(mǎn)足一些特殊性質(zhì)的非初等函數(shù)，常用的特殊函數(shù)有誤差函數(shù)、Gamma函數(shù)、超幾何函數(shù)、Bessel函數(shù)、Airy積分等。因其重要性，許多學(xué)者都致力于這些函數(shù)的賦值、性質(zhì)、應(yīng)用等相關(guān)研究。美國(guó)國(guó)家標(biāo)準(zhǔn)技術(shù)研究所(NIST)在2010年出版了NIST數(shù)學(xué)函數(shù)手冊(cè)[1]，同年發(fā)布在線(xiàn)的DLMF(NIST數(shù)學(xué)函數(shù)數(shù)字圖書(shū)館)[2]，其中包含了大量的特殊函數(shù)。

誤差函數(shù)作為一類(lèi)重要的特殊函數(shù)，對(duì)其快速精確的賦值，既涉及到數(shù)值計(jì)算方法的相關(guān)理論，也涉及到浮點(diǎn)計(jì)算環(huán)境下的具體實(shí)現(xiàn)。目前為止，盡管在理論上已經(jīng)有不少方法可以處理此類(lèi)函數(shù)的計(jì)算，但是在當(dāng)前廣泛使用的IEEE浮點(diǎn)計(jì)算環(huán)境下，由于不同的數(shù)值法需要的基本運(yùn)算操作次數(shù)不同，導(dǎo)致舍入誤差的大小也相差甚遠(yuǎn)。一般來(lái)說(shuō)，要達(dá)到相同的計(jì)算精度，操作次數(shù)較少、誤差較小的算法更容易被工程上采納[3]。

誤差函數(shù)及其互補(bǔ)誤差函數(shù)在概率、統(tǒng)計(jì)和熱傳導(dǎo)上非常重要，它們是一類(lèi)熱方程的解。對(duì)沒(méi)有初等表達(dá)式的特殊函數(shù)，其賦值必須采用數(shù)值上的近似算法，得到近似解。常用的近似法包括級(jí)數(shù)展開(kāi)式、漸近展開(kāi)式、Pade近似、連分?jǐn)?shù)近似[4-6]等。對(duì)于數(shù)值算法的選擇，在達(dá)到要求計(jì)算精度的前提下，必須通過(guò)有限步的基本算術(shù)運(yùn)算完成。為達(dá)到較快的計(jì)算速度，需要的操作次數(shù)也盡可能地少。為了實(shí)現(xiàn)誤差函數(shù)的快速準(zhǔn)確賦值，研究了級(jí)數(shù)近似求解法，并在IEEE浮點(diǎn)標(biāo)準(zhǔn)環(huán)境下[7-8]進(jìn)行數(shù)值實(shí)驗(yàn)。通過(guò)Taylor級(jí)數(shù)法和Chebyshev級(jí)數(shù)法的比較，發(fā)現(xiàn)Chebyshev級(jí)數(shù)法在計(jì)算誤差函數(shù)時(shí)，收斂速度遠(yuǎn)優(yōu)于Taylor級(jí)數(shù)法。在達(dá)到相同的計(jì)算精度要求時(shí)，Chebyshev級(jí)數(shù)法需要的項(xiàng)數(shù)要遠(yuǎn)少于Taylor級(jí)數(shù)法。為進(jìn)一步減少舍入誤差，在求級(jí)數(shù)的部分和時(shí)采用了Clenshaw算法。

1 誤差函數(shù)的定義及性質(zhì)

誤差函數(shù)erf(x)是一個(gè)全函數(shù)，定義為

圖1為函數(shù)erf(x)當(dāng)自變量在區(qū)間[-5,5]的圖像。誤差函數(shù)有對(duì)稱(chēng)關(guān)系

erf(-x)=-erf(x),

圖1 誤差函數(shù)erf(x)Fig.1 Error function erf(x)

其性質(zhì)有

erf(x)→1, x→,。

與誤差函數(shù)erf(x)相關(guān)的Dawson積分定義為

實(shí)際上，誤差函數(shù)和Dawson積分是不完全Gamma函數(shù)γ(a,x)的特例：

因此，誤差函數(shù)erf(x)與不完全Gamma函數(shù)的關(guān)系對(duì)推導(dǎo)誤差函數(shù)的展開(kāi)式起到了關(guān)鍵作用。誤差函數(shù)的Taylor級(jí)數(shù)展開(kāi)式的推導(dǎo)過(guò)程參考文獻(xiàn)[6]。本研究關(guān)注的是Chebyshev級(jí)數(shù)在誤差函數(shù)賦值上的應(yīng)用。

2 Taylor級(jí)數(shù)法和Chebyshev級(jí)數(shù)法

2.1 Taylor級(jí)數(shù)法

誤差函數(shù)erf(x)的Taylor級(jí)數(shù)展開(kāi)式為

注意到Taylor級(jí)數(shù)展開(kāi)式使用的是xk基函數(shù)，使用x的冪次作為基函數(shù)有其弊端，因?yàn)榫哂衳2k形式的函數(shù)都具有相同的特性，另一方面，具有x2k+1形式的函數(shù)也具有相同的行為。表面來(lái)看，基函數(shù)是線(xiàn)性相關(guān)的。

不使用x的冪次作為基函數(shù)，而使用正交基構(gòu)造多項(xiàng)式，在很多數(shù)值計(jì)算方法上體現(xiàn)出了優(yōu)勢(shì)[3]。Chebyshev級(jí)數(shù)法正是使用正交基作為基函數(shù)的方法。

2.2 Chebyshev級(jí)數(shù)法

2.2.1 Chebyshev多項(xiàng)式

正交基函數(shù)Ti(x)序列滿(mǎn)足：

多項(xiàng)式Ti(x)稱(chēng)為Chebyshev多項(xiàng)式，可作為多項(xiàng)式的正交基。Chebyshev級(jí)數(shù)使用Chebyshev多項(xiàng)式作為基函數(shù)，而不采用單項(xiàng)式xi作為基函數(shù)[9]。

多項(xiàng)式Ti(x)滿(mǎn)足遞推關(guān)系：

T0(x)=1；

T1(x)=x；

Ti+1(x)=2xTi(x)-Ti-1(x)， i≥1。

圖2為前6個(gè)Chebyshev多項(xiàng)式的圖像。Chebyshev多項(xiàng)式具有很多優(yōu)良的性質(zhì)：

圖2 Chebyshev多項(xiàng)式(i=0, 1,…, 5)Fig.2 Chebyshev polynomials(i=0, 1,…, 5)

1) Chebyshev多項(xiàng)式有一個(gè)閉式的表達(dá)式

Ti(x)=cos(iarccos x)。

3) Ti(x)的零點(diǎn)

每個(gè)次數(shù)為n的多項(xiàng)式都可用T0(x),T1(x),…,Tn(x)的線(xiàn)性組合表示。以下給出函數(shù)f(x)在區(qū)間[-1,1]上的Chebyshev級(jí)數(shù)展開(kāi)式計(jì)算方法，并擴(kuò)展到任意有定義的區(qū)間[a,b]。

2.2.2 Chebyshev級(jí)數(shù)展開(kāi)式

若函數(shù)f(x)在區(qū)間[-1,1]上有定義，且有連續(xù)的一階導(dǎo)數(shù)f′(x)，則函數(shù)f(x)有一個(gè)收斂的Chebyshev級(jí)數(shù)展開(kāi)式

Chebyshev級(jí)數(shù)的系數(shù)有閉式的形式，重寫(xiě)式f(x)為

系數(shù)τi可由式(1)確定，

(1)

做變量替換x=cos φ，式(1)可轉(zhuǎn)換為

(2)

當(dāng)沒(méi)有合適的微積分方法將式(2)轉(zhuǎn)換為實(shí)數(shù)時(shí)，需要采用數(shù)值近似法計(jì)算系數(shù)τi。

當(dāng)需要求解的函數(shù)f(x)是定義在區(qū)間[a,b]上時(shí)，需要做一個(gè)簡(jiǎn)單的變量替換

，

把區(qū)間[-1,1]轉(zhuǎn)換到區(qū)間[a,b]上。也就是說(shuō)，在求解系數(shù)τi時(shí)，由式(3)確定。

(3)

一旦得到f(x)在區(qū)間[-1,1]上的Chebyshev級(jí)數(shù)展開(kāi)式，那么計(jì)算f(y),y∈[a,b]時(shí)，轉(zhuǎn)換為計(jì)算

2.2.3 Chebyshev級(jí)數(shù)的截?cái)嗾`差

若定義在[-1,1]上的Chebyshev級(jí)數(shù)的系數(shù)τi滿(mǎn)足

那么對(duì)于滿(mǎn)足下式的n，

，

截?cái)嗾`差滿(mǎn)足關(guān)系式：

3 Clenshaw算法

Clenshaw算法(也稱(chēng)Clenshaw求和法)是一種計(jì)算Chebyshev多項(xiàng)式的線(xiàn)性組合的遞推方法[10]，是計(jì)算單項(xiàng)式的線(xiàn)性組合的Horner算法的推廣。Clenshaw不僅可運(yùn)用于Chebyshev多項(xiàng)式，也可應(yīng)用在任何由三項(xiàng)遞推關(guān)系定義的函數(shù)求和。

3.1 Clewshaw遞推公式

Clenshaw遞推公式是一種非常巧妙且有效的求和方法。計(jì)算模式由一系列系數(shù)和其相應(yīng)的函數(shù)相乘，然后對(duì)其結(jié)果求和。這些函數(shù)必須滿(mǎn)足一定的遞推關(guān)系。

假設(shè)需要對(duì)下式求和

其中Fk滿(mǎn)足遞推關(guān)系

Fk+1(x)=α(k,x)Fk(x)+β(k,x)Fk-1(x)。

α(k,x),β(k,x)是關(guān)于k和x的函數(shù)。

通過(guò)遞推關(guān)系定義：

把ck當(dāng)未知數(shù)，求解方程(4)，那么

納入標(biāo)準(zhǔn)：（1）患者肝、腎等功能正常。（2）患者有焦慮癥狀，SAS評(píng)分＞50分。（3）有正常語(yǔ)言表達(dá)和理解能力。

注意到式(5)最后一行中β(1,x)被加了一次，又被減了一次。若檢查所有含yk的因子，則會(huì)發(fā)現(xiàn)最終式(5)只剩

f(x)=β(1,x)F0(x)y2+F1(x)y1+F0(x)c0,

(6)

式(4)、(6)就是對(duì)f(x)進(jìn)行求和的Clenshaw遞推公式。

因?yàn)镃hebyshev多項(xiàng)式滿(mǎn)足三項(xiàng)遞推關(guān)系，因此，可以應(yīng)用Clenshaw算法計(jì)算Chebyshev級(jí)數(shù)的部分和。

3.2 Clenshaw算法求Chebyshev級(jí)數(shù)部分和

計(jì)算Chebyshev級(jí)數(shù)的前N項(xiàng)部分和

可轉(zhuǎn)化為計(jì)算

首先通過(guò)Clenshaw算法計(jì)算

，

由Chebyshev多項(xiàng)式的遞推關(guān)系及式(4)，可得

yN+2=yN+1=0，

yk=2xyk+1-yk+2+τk， k=N,N-1,…,1。

對(duì)于給定的N和x及計(jì)算得到的系數(shù)τk。根據(jù)上述的遞推關(guān)系得到y(tǒng)2、y1，然后由

計(jì)算出T(x)。

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

針對(duì)誤差函數(shù)erf(x)，采用Taylor級(jí)數(shù)和Chebyshev級(jí)數(shù)2種近似方法，在IEEE浮點(diǎn)標(biāo)準(zhǔn)計(jì)算環(huán)境下進(jìn)行數(shù)值實(shí)驗(yàn)。實(shí)驗(yàn)選擇誤差函數(shù)的自變量區(qū)間范圍為[0,5]，采樣點(diǎn)為1000，均勻分布在整個(gè)[0,5]區(qū)間上。即采樣點(diǎn)坐標(biāo)為

xi=i/200, i=1,2,…,1000。

第一個(gè)實(shí)驗(yàn)的目的是測(cè)試比較2種級(jí)數(shù)表示法的近似效果。第二個(gè)實(shí)驗(yàn)的目的是考查IEEE浮點(diǎn)標(biāo)準(zhǔn)環(huán)境下舍入誤差對(duì)數(shù)值計(jì)算產(chǎn)生的影響。

4.1 Taylor級(jí)數(shù)和Chebyshev級(jí)數(shù)近似計(jì)算erf(x)

表1為Chebyshev級(jí)數(shù)展開(kāi)式和Taylor級(jí)數(shù)展開(kāi)式在項(xiàng)數(shù)相同時(shí)的平均相對(duì)誤差。采用的計(jì)算精度為10位有效數(shù)字。

表1 Chebyshev級(jí)數(shù)和Taylor級(jí)數(shù)的比較

從表1可知，隨著項(xiàng)數(shù)的增加，Chebyshev級(jí)數(shù)法得到的計(jì)算結(jié)果的平均相對(duì)誤差越來(lái)越小，在展開(kāi)到19項(xiàng)時(shí)達(dá)到了8位有效數(shù)字。相比之下，Taylor級(jí)數(shù)法計(jì)算結(jié)果的平均相對(duì)誤差非常大，而且從表1實(shí)驗(yàn)的記錄來(lái)看，項(xiàng)數(shù)越多，則誤差越大，這是由于Taylor級(jí)數(shù)在近似計(jì)算誤差函數(shù)時(shí)收斂速度非常慢，且數(shù)值不穩(wěn)定。通過(guò)進(jìn)一步的實(shí)驗(yàn)可知，若要達(dá)到8位有效數(shù)字，則Taylor級(jí)數(shù)需要展開(kāi)的項(xiàng)數(shù)多達(dá)73項(xiàng)。表2為T(mén)aylor級(jí)數(shù)法在不同的展開(kāi)項(xiàng)數(shù)時(shí)平均相對(duì)誤差。

表2 Taylor級(jí)數(shù)法得到的平均相對(duì)誤差

4.2 舍入誤差對(duì)計(jì)算部分和的影響

在表2中，使用Taylor級(jí)數(shù)展開(kāi)式近似計(jì)算誤差函數(shù)，發(fā)現(xiàn)從80項(xiàng)增加到100項(xiàng)時(shí)，項(xiàng)數(shù)的增加并未導(dǎo)致相對(duì)誤差的減小，其主要原因是IEEE浮點(diǎn)計(jì)算環(huán)境下舍入誤差的影響。從級(jí)數(shù)的表達(dá)式可以看出，求部分和時(shí)需要計(jì)算的x的多個(gè)冪次x2i+1,i=1,2,…,N，項(xiàng)數(shù)越多，冪次越高，需要進(jìn)行的基本運(yùn)算次數(shù)也越多。這樣的情況下，舍入誤差是一個(gè)不能忽略的影響計(jì)算結(jié)果的因素。

表3是計(jì)算精度為10位有效數(shù)字和30位有效數(shù)字時(shí)對(duì)計(jì)算Taylor級(jí)數(shù)部分和的影響。從表3可見(jiàn)，在增加中間計(jì)算的精度后，近似效果得到大幅度地提升。

表3 舍入誤差對(duì)計(jì)算Taylor級(jí)數(shù)部分和的影響

4.3 討論

通過(guò)對(duì)Taylor級(jí)數(shù)和Chebyshev級(jí)數(shù)在近似求解誤差函數(shù)時(shí)平均相對(duì)誤差的比較，發(fā)現(xiàn)Chebyshev級(jí)數(shù)能在展開(kāi)項(xiàng)數(shù)較少的情況下，獲得較高的精度。若需要描繪誤差函數(shù)的一個(gè)粗略的圖像，則使用具有2位(及以上)有效數(shù)字的近似展開(kāi)式就可得到較好的效果。圖3為N=7時(shí)Chebyshev級(jí)數(shù)展開(kāi)式得到的區(qū)間[-1,1]的圖像和高精度計(jì)算得到的誤差函數(shù)在區(qū)間[0,5]的圖像。從圖3可以看出，Chebyshev級(jí)數(shù)法實(shí)際上是將定義域在[0,5]的函數(shù)圖像平移并壓縮到了區(qū)間[-1,1]。

圖3 Chebyshev級(jí)數(shù)法圖像與高精度erf(x)圖像Fig.3 Image by Chebyshev series and high precision image of erf(x)

5 結(jié)束語(yǔ)

截?cái)嗾`差和舍入誤差是計(jì)算一個(gè)復(fù)雜算術(shù)表達(dá)式時(shí)出現(xiàn)的2種誤差來(lái)源，對(duì)序列的收斂速度進(jìn)行分析可計(jì)算截?cái)嗾`差，舍入誤差本質(zhì)上是因?yàn)橛?jì)算機(jī)表示的實(shí)數(shù)受限于尾數(shù)的固定精度，而且計(jì)算機(jī)硬件只支持有限位機(jī)器數(shù)的運(yùn)算，導(dǎo)致舍入誤差在計(jì)算機(jī)中被引入和傳播，舍入誤差比截?cái)嗾`差更難于分析和控制，計(jì)算時(shí)必須同時(shí)考慮這兩種誤差。值得一提的是，Chebyshev級(jí)數(shù)展開(kāi)式的系數(shù)是需要計(jì)算積分式的，當(dāng)項(xiàng)數(shù)較多時(shí)，會(huì)花費(fèi)相當(dāng)多的時(shí)間求解系數(shù)，從而導(dǎo)致計(jì)算速度迅速下降。而Taylor級(jí)數(shù)只需要基本算術(shù)運(yùn)算就能完成，顯示了其具有的優(yōu)勢(shì)。

[1] OLVER F W J,LOZIER D W,BOISVERT R F,et al.NIST Handbook of Mathematical Functions[M].New York: Cambridge University Press,2010:160-170.

[2] OLVER F W J,LOZIER D W,BOISVERT R F,et al.NIST digital library of mathematical functions[EB/OL].[2016-09-05].http://dlmf.nist.gov/.

[3] HIGHAM N J.Accuracy and Stability of Numerical Algorithms[M]. Philadelphia:Society for Industrial and Applied Mathematics(SIAM),2002:678-698.

[4] OLVER F W J.Asymptotics and Special Functions[M].Wellesley:A K Peters,1997:145-163.

[5] GIL A,SEGURA J,TEMME N M.Numerical Methods for Special Functions[M].Philadelphia:Society for Industrial and Applied Mathematics(SIAM),2007:51-83.

[6] CUYT A,PETERSEN V B,VERDONK B,et al.Handbook of Continued Fractions for Special Functions[M].Birkhauser Boston: Springer,2008:253-259.

[7] Floating-Point Working Group. IEEE standard for binary floating-point arithmetic[J].SIGPLAN,1987,22:9-25.

[8] MULLER J M,BRISEBARRE N,DINECHIN F,et al.Handbook of Floating-Point Arithmetic[M].Birkhauser Boston:Springer,2010:55-58.

[9] MATHEWS J H,FINK K D.Numerical Methods Using Matlab[M].New York:Pearson,2004:365-389.[10] PRESS W H,FLANNERY B P.Numerical Recipes in FORTRAN 77-The Art of Scientific Computing[M].England:Cambridge University,1992:34-45.

編輯：梁王歡

Chebyshev series method for the error function

DENG Guoqiang1,2, TANG Min1,2

(1.School of Mathematics and Computing Science, Guilin University of Electrical Technology, Guilin 541004, China;2.Guangxi Colleges and Universities Key Laboratory of Data Analysis and Computation,Guilin University of Electrical Technology, Guilin 541004, China)

In order to evaluate error function efficiently and accurately in the IEEE floating-point arithmetic environment. This paper investigates the Chebyshev series expansion method for the evaluation of the error function. In addition, using Clenshaw algorithm to compute the partial sum of the firstNterms leads to the rounding errors reduction. The experimental results show that Chebyshev series method has faster convergence speed than Taylor series method. In other words, to attain the same computation precision Chebyshev method needs fewer terms than Taylor method.

Error function; Chebyshev series; Clenshaw algorithm; IEEE floating-point arithmetic standard

2016-03-05

國(guó)家自然科學(xué)基金(11561015)；廣西自然科學(xué)基金(2016GXNSFFA380009)

唐敏(1980-)，女(壯族)，廣西桂林人，講師，博士研究生，研究方向?yàn)楦呖尚庞?jì)算、浮點(diǎn)計(jì)算。E-mail:52121500011@ecnu.cn

鄧國(guó)強(qiáng)，唐敏.誤差函數(shù)Chebyshev級(jí)數(shù)的計(jì)算方法[J].桂林電子科技大學(xué)學(xué)報(bào)，2016,36(6)：508-512.

O241.1

A

1673-808X(2016)06-0508-05