鄧國強(qiáng)，唐 敏

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

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

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

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

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

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

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

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

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

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

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

圖1為函數(shù)erf(x)當(dāng)自變量在區(qū)間[-5,5]的圖像。誤差函數(shù)有對稱關(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ù)和Dawson積分是不完全Gamma函數(shù)γ(a,x)的特例：

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

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

2.1 Taylor級數(shù)法

誤差函數(shù)erf(x)的Taylor級數(shù)展開式為

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

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

2.2 Chebyshev級數(shù)法

2.2.1 Chebyshev多項式

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

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

多項式Ti(x)滿足遞推關(guān)系：

T0(x)=1；

T1(x)=x；

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

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

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

1) Chebyshev多項式有一個閉式的表達(dá)式

Ti(x)=cos(iarccos x)。

3) Ti(x)的零點

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

2.2.2 Chebyshev級數(shù)展開式

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

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

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

(1)

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

(2)

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

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

，

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

(3)

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

2.2.3 Chebyshev級數(shù)的截斷誤差

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

那么對于滿足下式的n，

，

截斷誤差滿足關(guān)系式：

3 Clenshaw算法

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

3.1 Clewshaw遞推公式

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

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

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

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

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

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

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

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

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

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

(6)

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

因為Chebyshev多項式滿足三項遞推關(guān)系，因此，可以應(yīng)用Clenshaw算法計算Chebyshev級數(shù)的部分和。

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

計算Chebyshev級數(shù)的前N項部分和

可轉(zhuǎn)化為計算

首先通過Clenshaw算法計算

，

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

yN+2=yN+1=0，

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

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

計算出T(x)。

4 數(shù)值實驗

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

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

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

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

表1為Chebyshev級數(shù)展開式和Taylor級數(shù)展開式在項數(shù)相同時的平均相對誤差。采用的計算精度為10位有效數(shù)字。

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

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

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

4.2 舍入誤差對計算部分和的影響

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

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

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

4.3 討論

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

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

5 結(jié)束語

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

[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

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

唐敏(1980-)，女(壯族)，廣西桂林人，講師，博士研究生，研究方向為高可信計算、浮點計算。E-mail:52121500011@ecnu.cn

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

O241.1

A

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