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

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

0 引 言

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

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

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

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

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

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

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

1 RBM

非多項式函數(shù)的實根隔離是一項艱巨的工作，超出了本文的討論范圍，當(dāng)f(t)為多項式時，可通過實根隔離方法使每個區(qū)間只含有一個實根［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)對應(yīng)的無窮級數(shù)在區(qū)間[a，b]內(nèi)收斂，即

通過漸進(jì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)：

（?。┯糜欣矶囗検絃(s)/gi(s)作為Gi(s)的形式，其中L(s)為滿足L(a)=f(a)和L(b)=f(b)的一次多項式，則s*對應(yīng)為L(s)的根，可較方便地求解相應(yīng)值。

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

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

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

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

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

1.2 算法概述

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

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

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

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

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

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

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

式（10）的求根步驟為

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

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

輸出實根t*的近似值。

步驟 1令t2，i=2。

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

步驟3由式（10），計算φ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 時，有

由式（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）等價于

即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)的實根t*。

1.4 RBM 示例

先用簡單且方便驗證的例子說明本文算法。

例 1令f1(t)=(t?0.25)(2 ?t)(t+5)2，t∈[0，1]，在區(qū)間［0，1］內(nèi)有一個單根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ù)點后100。

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

2 不同方法的定性比較

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

表2 3 類方法的符號Table 2 Symbols of three types of methods

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

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

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

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

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

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

3 實例及討論

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

3.1 M1與裁剪法的比較

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

例2測試以下4 個多項式函數(shù)：

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

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

表6 例2 中4 個多項式函數(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的每個迭代步分別需花費6 和4 次FE，漸進(jìn)法每個迭代步需花費6 次FE，總共需化費12 次FE。

例3測試3 個非多項式函數(shù)，比較漸進(jìn)法M1，M2，M3和類牛頓法N2，N3。3 個 非多項式函數(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知，可用一個初始值得到一個包含實根的小區(qū)間。表8 所列為每步的逼近誤差和收斂階。由表8 可知，對于f6（t）和f7（t），N2和N3收斂至正確結(jié)果；而對于f8（t），N2在ta=5.334 7 附近發(fā)散，N3則收斂至錯誤解tb=16.933。3 種情況下，漸進(jìn)法均能收斂至正確結(jié)果。因此，漸進(jìn)法M1，M2和M3的收斂階相差不大，且均較類牛頓法N2和N3的收斂速度快得多。表9為在不同有效位數(shù)下運行12 次FE 的總計算時間。由表9 可知，M1，N2和N3的計算效率相當(dāng)，均較M2和M3好。M1在相同時間內(nèi)獲得最小的逼近誤差，并在給定的相同誤差下，計算效率最高。

表7 4 個多項式函數(shù)在不同方法下的裁剪步數(shù)nc和計算時間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 個非多項式函數(shù)在不同方法下的逼近誤差比較Table 8 Comparisons of approximation errors of 3 non-polynomial functions in different methods

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

3.3 多根情況討論

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

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

例4用RBM 計算Wilkinson 多項式

在區(qū)間［0，25］內(nèi)的實根。該Wilkinson 多項式有20個實根Wi，i=1，2，…，20（參見文獻(xiàn)［15］中的例6）。首先，計算相應(yīng)控制多邊形的零點，即｛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 個子區(qū)間，其中有16 個子區(qū)間包含W（t）的1 個或2 個根。本文選取其中的3 個子區(qū)間Λ1=［0.27，1.55］、Λ2=［2.83，4.11］、Λ3=［16.81，18.07］用以說明更多細(xì)節(jié)。此3 個子區(qū)間分別包含1 個、2 個、2 個根?？蓪BM 直接應(yīng)用于Λ1，因為W（t）在其2 個端點的函數(shù)值異號；對于Λ2和Λ3，可選擇區(qū)間中點作為分界點，將每個區(qū)間拆分為2 個子區(qū)間，每個子區(qū)間都包含1 個根。然后，將RBM 應(yīng)用于子區(qū)間，由注1的方法獲得優(yōu)化區(qū)間。本例中，區(qū)間［0.27，1.55］、［2.83，4.11］、［16.81，18.07］對應(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 處理多根情形也具有較好的計算穩(wěn)定性。

圖1 Wilkinson 多項式的形狀Fig.1 The plots of Wilkinson polynomial

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

4 結(jié) 論

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

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

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