光滑函數(shù)實(shí)根計(jì)算的漸進(jìn)顯式公式

祝平，陳小雕，*，馬維銀，姜霓裳

（1.杭州電子科技大學(xué)計(jì)算機(jī)學(xué)院，浙江杭州 310018； 2.香港城市大學(xué) 機(jī)械工程學(xué)系，中國香港； 3.杭州電子科技大學(xué)理學(xué)院，浙江杭州 310018）

0 引 言

求根在計(jì)算機(jī)輔助幾何設(shè)計(jì)［1］、計(jì)算機(jī)圖形學(xué)［2-3］、機(jī)器人技術(shù)［4］及地磁導(dǎo)航中具有廣泛應(yīng)用。理論上，利用結(jié)式等消元方法，可將多項(xiàng)式方程組轉(zhuǎn)化為一元多項(xiàng)式方程。本文只討論單變量方程的求根問題，線性方程組的其他求解方法可參見文獻(xiàn)［1，5-6］等。

計(jì)算穩(wěn)定性與計(jì)算效率是求根的2 個(gè)關(guān)鍵問題。效率因子（efficiency index，EI）常用于衡量算法的計(jì)算效率。設(shè)某算法的收斂階為p，需進(jìn)行n次函數(shù)計(jì)算（functional evaluation，F(xiàn)E），則對(duì)應(yīng)的EI 定義為p1/n。文獻(xiàn)［7］猜想：任何需要n次FE 的求根方法對(duì)應(yīng)的收斂階都不會(huì)超過最優(yōu)次數(shù)2n-1。

下文著重介紹求解光滑函數(shù)實(shí)根的3 種方法。

（1） 類牛頓法，包走 Newton，Halley 和Ostrowski 法［8-10］。從初始值開始，獲得近似于根的點(diǎn)序列，通常用第k個(gè)點(diǎn)的信息計(jì)算第(k+1)個(gè)點(diǎn)，計(jì)算性能取決于初始值的選取，當(dāng)初始值選取錯(cuò)誤時(shí)，易發(fā)散。文獻(xiàn)［9-10］以收斂階8 得到了最佳效率指數(shù)，但不易擴(kuò)展至更高收斂階［11］。理論上相應(yīng)的區(qū)間方法可提高計(jì)算穩(wěn)定性，但將耗費(fèi)更多的計(jì)算時(shí)間［12-13］?；诘芽柗▌t與Sturm 定理的多項(xiàng)式求根方法雖具有較好的計(jì)算穩(wěn)定性，但收斂速度較低。

（2）基于Bernstein-Bézier 形式的裁剪法［1，14-17］。利用Bernstein-Bézier 形式具有的系數(shù)擾動(dòng)穩(wěn)定性特點(diǎn)，通過不斷裁剪不包含實(shí)根的小區(qū)間，獲取較好的計(jì)算穩(wěn)定性［18］，收斂速度通常優(yōu)于Sturm方法［19］。

（3）漸進(jìn)法［20-21］。從包含根的初始區(qū)間開始，利用前k個(gè)點(diǎn)（而非僅第k個(gè)點(diǎn)）的信息漸進(jìn)式計(jì)算第(k+1)個(gè)點(diǎn)。其基本思想是利用分子的線性多項(xiàng)式逼近給定的光滑函數(shù)，從而逼近多項(xiàng)式對(duì)應(yīng)的實(shí)根，得到較精確的解。對(duì)于區(qū)間內(nèi)只包含一個(gè)根的情況，逼近多項(xiàng)式的實(shí)根被限制在給定區(qū)間，且近似效果漸進(jìn)式遞增，即使沒有好的初始值也可確保收斂。

本文提出了一種基于重新參數(shù)化方法（reparamaterization-based method，RBM）的非線性方程求根方法，直接尋求重新參數(shù)化函數(shù)并給出實(shí)根的漸進(jìn)式近似公式，得到了較現(xiàn)有漸進(jìn)法更高的計(jì)算效率，并兼具漸進(jìn)法的計(jì)算穩(wěn)定性。

1 RBM

非多項(xiàng)式函數(shù)的實(shí)根隔離是一項(xiàng)艱巨的工作，超出了本文的討論范圍，當(dāng)f(t)為多項(xiàng)式時(shí)，可通過實(shí)根隔離方法使每個(gè)區(qū)間只含有一個(gè)實(shí)根［1，13，22-23］。本文假設(shè)函數(shù)f(t)在給定區(qū)間[a，b]內(nèi)具有唯一單根t*。

1.1 思路及算法原理

首先 ，設(shè)函數(shù)Ri(t) 插值f(t) 于t=t0，t1，…，ti∈[a，b]，即

其中，t0=a和t1=b。

則Ri(t)和f(t)對(duì)應(yīng)的無窮級(jí)數(shù)在區(qū)間[a，b]內(nèi)收斂，即

通過漸進(jìn)式插值更多的點(diǎn)可得到更小的逼近誤差。

設(shè)φi(s)為重新參數(shù)化函數(shù)，使得φi(sj)=tj且

若s*∈[a，b]滿足Gi(s*)=0，則類似于式（1），可得

由式（2），φi(s*)可作為t*的近似值。

然后，做以下技術(shù)處理，以方便高效地獲取相應(yīng)的s*和重新參數(shù)化函數(shù)φi(s)：

（ⅰ）用有理多項(xiàng)式L(s)/gi(s)作為Gi(s)的形式，其中L(s)為滿足L(a)=f(a)和L(b)=f(b)的一次多項(xiàng)式，則s*對(duì)應(yīng)為L(s)的根，可較方便地求解相應(yīng)值。

（ⅱ）用重新參數(shù)化函數(shù)φi(s)作為次數(shù)為i的多項(xiàng) 式，并 用(s/gi(s)，L(s)/gi(s))插 值(t，f(t))。于是，每 給 定 一 個(gè)tj，由(sj/gi(sj)，L(sj)/gi(sj))=(tj，f(tj))，可得

式（3）是關(guān)于sj的一次多項(xiàng)式，可方便地求解sj對(duì)應(yīng)的值。給定t0，t1，…，ti，由式（3），可得到對(duì)應(yīng)的s0，s1，…，si。進(jìn)而由i+1 個(gè)線性方程：

確定相應(yīng)的多項(xiàng)式φi(s)。

最后，由式（2），得φi(s*)的近似實(shí)根t*。

本文方法無需計(jì)算gi(s)的表達(dá)式，計(jì)算效率比現(xiàn)有的漸進(jìn)法更高，且兼具漸進(jìn)法良好的計(jì)算穩(wěn)定性。

1.2 算法概述

引入文獻(xiàn)［24］中的定理3.5.1，即

定理1令w0，w1，…，wr為區(qū)間［a，b］內(nèi)的r+1個(gè)不同點(diǎn)，并且n0，w1，…，wr為r+1 個(gè)大于0 的整數(shù)，令N=n0+n1+…+nr，假設(shè)g(t)是次數(shù)為N的多項(xiàng)式，使得g(i)(wj)=f(i)(wj)，i=0，1，…，nj?1，j=0，1，…，r。存在ξ0(t)∈[a，b]，使得

設(shè)t0=a，t1=b，可驗(yàn)證

求解式（6），得到ti+1。

其中，s0=t0，s1=t1，sj為線性方程

的實(shí)根。令φi(s)為i次多項(xiàng)式，i=1，2…，使得

由式（7）和式（9），得到漸近式

式（10）的求根步驟為

算法1求解區(qū)間[a，b]內(nèi)f(t)的實(shí)根t*。

輸入函數(shù)f(t)、區(qū)間[a，b]和誤差精度ε。

輸出實(shí)根t*的近似值。

步驟 1令t2，i=2。

步驟2若轉(zhuǎn)步驟4。否則，轉(zhuǎn)步驟3。

步驟3由式（10），計(jì)算φi(s)和并令i=i+1，轉(zhuǎn)步驟2。

步驟4輸出即t*的近似值。

1.3 收斂階

由式（8）和式（11），可得

定理2設(shè)則有

其中，h=b?a，為區(qū)間[a，b]的長度。

證明由式（13），有

由式（15），當(dāng)i=2 時(shí)，有

由式（7），可得

由式（15）和式（17），可得

由式（16）和式（18），可得式（14）。

證畢。

定 理3設(shè)κ1，κ2∈[a，b]，若2|κ2?t*|＜|κ1?t*|，則2κ2?κ1和κ1將包住t*。

證明不失一般性，不妨設(shè)κ1＜κ2，則有

即

由κ1＜κ2和式（19），可得

式（20）等價(jià)于

即2κ2?κ1和κ1包住了t*。

證畢。

注1在式（10）中，所有ti和si（i=2，3，…）應(yīng)落在給定區(qū)間[a，b]內(nèi)，但在某些情況下，ti或si可能落在區(qū)間[a，b]外。由定理3 知，通過獲得更小的子區(qū)間，并用算法1 得到相應(yīng)的實(shí)根t*。

1.4 RBM 示例

先用簡單且方便驗(yàn)證的例子說明本文算法。

例 1令f1(t)=(t?0.25)(2 ?t)(t+5)2，t∈[0，1]，在區(qū)間［0，1］內(nèi)有一個(gè)單根t*=0.25。令可得線性函數(shù)L(s)=?12.5+39.5s及s*=由式（3），可得由式（6）和式（3），可分別得到t3=φ2(s*)≈0.040 274，s3≈0.041 84。

更多ei=|ti?t*|的值（i=3，4，…，9)，見表1。本例中，精度設(shè)置為小數(shù)點(diǎn)后100。

表1 例1 中不同i 值下逼近誤差ei=|ti?t*|的對(duì)應(yīng)結(jié)果Table 1 The corresponding results of the approximation error ei=|ti?t*|in example 1 when different i

2 不同方法的定性比較

表2 給出了3 類方法的符號(hào)，本節(jié)將對(duì)此3 類方法進(jìn)行定性比較。所有示例均在具有12 GB 內(nèi)存和4 核2.2 G CPU 的PC 上 用Mathematica 軟 件 進(jìn) 行 測(cè)試，時(shí)間單位為ms。

表2 3 類方法的符號(hào)Table 2 Symbols of three types of methods

首先，比較不同方法的漸進(jìn)效率因子（AEI）。AEI 與EI 類似，通常用于測(cè)量算法的計(jì)算效率，其含義為每增加一次FE 所對(duì)應(yīng)的收斂階，AEI 越大效率越高。表3 為不同方法的計(jì)算效率比較，其中nfe、CR、AEI 分別表示FE 次數(shù)、收斂階和漸近效率因子。表3 中，從第2 步或第3 步開始，漸進(jìn)法M1，M2和M3每增加一次FE 可獲取的收斂階為2。結(jié)果表明，漸進(jìn)法較類牛頓法和裁剪法具有更高的計(jì)算效率。

表3 3 類方法的計(jì)算效率比較Table 3 Comparisons of computational efficiency of three types of methods

其次，表4 中，nfe為不同方法達(dá)到收斂階為16 時(shí)所對(duì)應(yīng)的FE 次數(shù)。理論上，C1，C2和C3分別需花費(fèi)4n，7n和5n次FE 才能獲得收斂階4n，7n，5n，且每步裁剪將分別增加4，7 和5 次FE；而本文方法M1花費(fèi)n次FE 可 達(dá) 到 的 收 斂 階 為3 ?2n?2(n=3，4，…)，每增加一次FE 可達(dá)到的收斂階為2。因此，M1較Ci（i=1，2，3）更靈活和高效。

最后，表5 為AEI 為2 時(shí)不同漸進(jìn)法Mi(i=1，2，3)在第n步的計(jì)算成本比較，其中nfe，na，nm分別表示FE 的次數(shù)、加（減）法運(yùn)算的次數(shù)和乘（除）法運(yùn)算的次數(shù)。由表5 可知，M1比M2和M3的計(jì)算效率更高，且當(dāng)n較大時(shí)越為有利。

表4 收斂階為16 時(shí)對(duì)應(yīng)的FE 次數(shù)的比較Table 4 Comparison results on FE times under convergence rate 16

表5 AEI 為2 時(shí)Mi在第n 步的計(jì)算成本比較Table 5 Comparison of computational cost of the n-th step in Mi under AEI 2

3 實(shí)例及討論

用3 個(gè)實(shí)例，對(duì)不同方法進(jìn)行了比較?？傮w來說，M1具有與裁剪法相當(dāng)?shù)挠?jì)算穩(wěn)定性和更高的計(jì)算效率、比類牛頓法更高的計(jì)算穩(wěn)定性和更高的計(jì)算效率、與漸進(jìn)法相當(dāng)?shù)氖諗侩A和更高的計(jì)算效率。

3.1 M1與裁剪法的比較

對(duì)M1與裁剪法C2［15］和C3［16］進(jìn)行了比較。理論上，C2和C3可以計(jì)算多項(xiàng)式在給定區(qū)間內(nèi)的所有實(shí)根。計(jì)算多項(xiàng)式在某個(gè)區(qū)間內(nèi)的唯一單根，M1的效果更好，可將其作為裁剪法的補(bǔ)充。C2［15］和C3［16］在每個(gè)裁剪步驟中分別需要7 和5 次FE，并且第i個(gè)裁剪步驟子區(qū)間的長度為ei(i=1，2，…)。表6 中，M1每步需花費(fèi)5 次FE，第i個(gè)誤差ei映射至|t5i?1?t*|，使得M1中有5i次FE。平均計(jì)算時(shí)間與精度有關(guān)，在本文中，若無特別說明，精度為小數(shù)點(diǎn)后20 位。為更準(zhǔn)確地測(cè)試對(duì)應(yīng)的收斂階，將最大精度設(shè)置為小數(shù)點(diǎn)后3 000 位。

例2測(cè)試以下4 個(gè)多項(xiàng)式函數(shù)：

區(qū)間［0，1］內(nèi)分別有1 個(gè)單根1/4，1/5，1/3 和2/3。表6 為4 個(gè)多項(xiàng)式函數(shù)的近似誤差和計(jì)算時(shí)間，其中，計(jì)算時(shí)間和CR 分別表示每個(gè)迭代步驟的平均計(jì)算時(shí)間（精度為小數(shù)點(diǎn)后100 位以下，單位為ms）和平均收斂階。M1、C2和C3的收斂階CR 分別為32，7，5。M1每步的計(jì)算效率是C2、C3的10～50 倍。對(duì)于f2（t），M1的e3可達(dá)5.1×10?2011，對(duì)應(yīng)的有效精度位數(shù)為［20，3 000］。結(jié)果表明，M1的性能較C2和C3更佳。

給定誤差為10?16，測(cè)試了例2，結(jié)果見表7。其中，nc表示所需的（裁剪）步數(shù)，Tc表示總的計(jì)算時(shí)間（ms）。表明M1，C2和C3可在2 個(gè)（裁剪）步驟內(nèi)達(dá)到給定的誤差（10-16）要求，并且M1的計(jì)算效率更高，約為C2和C3的7～30 倍。由于M1為漸進(jìn)法，在獲取tk時(shí)，可隨時(shí)在任意k處停止，所需步驟較少，因此計(jì)算效率較高。

表6 例2 中4 個(gè)多項(xiàng)式函數(shù)在不同方法下的誤差和計(jì)算時(shí)間比較Table 6 Comparisons of the error and computational time of 4 polynomial functions in different methods of example 2

3.2 漸進(jìn)法與類牛頓法的比較

將漸進(jìn)法（M1，M2［20］和M3［21］）與類牛頓法（N2［11］和N3［10］）進(jìn) 行 了 比 較。N2和N3的每個(gè)迭代步分別需花費(fèi)6 和4 次FE，漸進(jìn)法每個(gè)迭代步需花費(fèi)6 次FE，總共需化費(fèi)12 次FE。

例3測(cè)試3 個(gè)非多項(xiàng)式函數(shù)，比較漸進(jìn)法M1，M2，M3和類牛頓法N2，N3。3 個(gè) 非多項(xiàng)式函數(shù)分別為f6（t）=（t?1/2）［esin（10（t?π））+4（t?π）?1］，t∈［3，3.3］，f7（t）=tanh（2t?π/25），t∈［?0.5，0.5］，f8（t）=?1/t+sin（t）+1，t∈［0.01，1.3］。其單根分別為由定理3知，可用一個(gè)初始值得到一個(gè)包含實(shí)根的小區(qū)間。表8 所列為每步的逼近誤差和收斂階。由表8 可知，對(duì)于f6（t）和f7（t），N2和N3收斂至正確結(jié)果；而對(duì)于f8（t），N2在ta=5.334 7 附近發(fā)散，N3則收斂至錯(cuò)誤解tb=16.933。3 種情況下，漸進(jìn)法均能收斂至正確結(jié)果。因此，漸進(jìn)法M1，M2和M3的收斂階相差不大，且均較類牛頓法N2和N3的收斂速度快得多。表9為在不同有效位數(shù)下運(yùn)行12 次FE 的總計(jì)算時(shí)間。由表9 可知，M1，N2和N3的計(jì)算效率相當(dāng)，均較M2和M3好。M1在相同時(shí)間內(nèi)獲得最小的逼近誤差，并在給定的相同誤差下，計(jì)算效率最高。

表7 4 個(gè)多項(xiàng)式函數(shù)在不同方法下的裁剪步數(shù)nc和計(jì)算時(shí)間Tc比較（誤差精度為10?16）Table 7 Comparisons of clipping steps nc and computational time Tc of 4 polynomial functions under different methods（error tolerance is 10?16）

上述結(jié)果表明，在方法M1，M2，M3，N2和N3中，M1性能最好。

表8 3 個(gè)非多項(xiàng)式函數(shù)在不同方法下的逼近誤差比較Table 8 Comparisons of approximation errors of 3 non-polynomial functions in different methods

表9 3 個(gè)非多項(xiàng)式函數(shù)在不同方法下運(yùn)行12 次FE 的時(shí)間比較Table 9 Comparisons of the computational time of three functions with cost of twelve functional evaluations under different methods 單位：ms

3.3 多根情況討論

（1）考慮如何獲取包含一個(gè)根的區(qū)間。設(shè)在初始值x0下，x1通過牛頓法或其他迭代方法得到。若2|x1?t*|＜|x0?t*|，t*是f（t）的根，由定理3，t*被x0和2x1?x0所包圍。因此，通過選取適當(dāng)?shù)某跏贾祒0可得到一個(gè)包含實(shí)根的小區(qū)間。

（2）考慮區(qū)間內(nèi)包含重根的情況，k為實(shí)根的重?cái)?shù)，k≥2。若k為偶數(shù)，利用Ai（s）=（s/gi（s），L（s）/gi（s））逼近C1（t）=(f(t)，f′(t))；否則，利用Ai（s）逼近C2（t）=(f′(t)，f(t))。理論上，估算k值且利用Ai（s）逼 近（f（k?2）（t），f（k?1）（t））是最好的選擇，有利于提高收斂階。

（3）考慮區(qū)間內(nèi)可能包含2 個(gè)或多個(gè)根的情況。若給定多項(xiàng)式函數(shù)，則一種可行的方法是將其轉(zhuǎn)換為Bernstein-Bézier 形式，并利用其控制多邊形的零點(diǎn)將給定的區(qū)間分為若干個(gè)子區(qū)間；對(duì)每個(gè)子區(qū)間，用與（1）相同的方法進(jìn)行實(shí)根隔離。對(duì)于非多項(xiàng)式情況，可將區(qū)間采樣為多個(gè)較小的子區(qū)間，并用與（1）相同的方法將其迭代拆分為若干個(gè)包含一個(gè)根的子區(qū)間。

例4用RBM 計(jì)算Wilkinson 多項(xiàng)式

在區(qū)間［0，25］內(nèi)的實(shí)根。該Wilkinson 多項(xiàng)式有20個(gè)實(shí)根Wi，i=1，2，…，20（參見文獻(xiàn)［15］中的例6）。首先，計(jì)算相應(yīng)控制多邊形的零點(diǎn)，即｛0.27，1.55，2.83，4.11，5.38，6.65，7.92，9.19，10.46，11.73，13.007，14.27，15.54，16.81，18.07，19.34，20.61，21.87，23.14，24.40｝，并將區(qū)間［0，25］劃分為20 個(gè)子區(qū)間，其中有16 個(gè)子區(qū)間包含W（t）的1 個(gè)或2 個(gè)根。本文選取其中的3 個(gè)子區(qū)間Λ1=［0.27，1.55］、Λ2=［2.83，4.11］、Λ3=［16.81，18.07］用以說明更多細(xì)節(jié)。此3 個(gè)子區(qū)間分別包含1 個(gè)、2 個(gè)、2 個(gè)根?？蓪BM 直接應(yīng)用于Λ1，因?yàn)閃（t）在其2 個(gè)端點(diǎn)的函數(shù)值異號(hào)；對(duì)于Λ2和Λ3，可選擇區(qū)間中點(diǎn)作為分界點(diǎn)，將每個(gè)區(qū)間拆分為2 個(gè)子區(qū)間，每個(gè)子區(qū)間都包含1 個(gè)根。然后，將RBM 應(yīng)用于子區(qū)間，由注1的方法獲得優(yōu)化區(qū)間。本例中，區(qū)間［0.27，1.55］、［2.83，4.11］、［16.81，18.07］對(duì)應(yīng)的優(yōu)化區(qū)間為Γ1=［0.99，1.0034］、Γ2=［2.83，3.004］、Γ3=［17.988，18.07］。如表10 所示，在誤差ei映射到重新參數(shù)化函數(shù)φi(t)的情況下，RBM 的效果較好，表明用RBM 處理多根情形也具有較好的計(jì)算穩(wěn)定性。

圖1 Wilkinson 多項(xiàng)式的形狀Fig.1 The plots of Wilkinson polynomial

表10 例4 中M1在3 個(gè)子區(qū)間內(nèi)的逼近誤差Table 10 The approximation errors of M1 in the three sub-intervals of example 4

4 結(jié) 論

提出了一種基于RBM 的高效求根方法，用于計(jì)算光滑函數(shù)在給定區(qū)間的實(shí)根。提供了一種能找到重新參數(shù)化函數(shù)的簡單方法，實(shí)現(xiàn)漸進(jìn)式逼近精確實(shí)根。與漸進(jìn)法M2和M3相比，本文方法M1具有與之相同的收斂階、更高的計(jì)算效率。與類牛頓法N2和N3相比，在相同的FE 下，本文方法M1的計(jì)算時(shí)間與之相近，逼近誤差更小、收斂階更高、計(jì)算穩(wěn)定性更高。結(jié)合裁剪法C2和C3的優(yōu)點(diǎn)（可以隔離給定區(qū)間內(nèi)多項(xiàng)式的所有實(shí)根），M1可用于計(jì)算給定函數(shù)在某給定區(qū)間內(nèi)的所有實(shí)根，并具有更高的計(jì)算效率和更高的收斂階。

未來的工作，包括進(jìn)一步擴(kuò)展RBM，選取合適的初始值，并將其應(yīng)用于非線性方程或方程組的求根。

感謝匿名專家對(duì)本文的審理和提供的寶貴意見！