光滑函數實根計算的漸進顯式公式

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

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

0 引 言

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

計算穩定性與計算效率是求根的2 個關鍵問題。效率因子（efficiency index，EI）常用于衡量算法的計算效率。設某算法的收斂階為p，需進行n次函數計算（functional evaluation，FE），則對應的EI 定義為p1/n。文獻［7］猜想：任何需要n次FE 的求根方法對應的收斂階都不會超過最優次數2n-1。

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

（1） 類牛頓法，包走 Newton，Halley 和Ostrowski 法［8-10］。從初始值開始，獲得近似于根的點序列，通常用第k個點的信息計算第(k+1)個點，計算性能取決于初始值的選取，當初始值選取錯誤時，易發散。文獻［9-10］以收斂階8 得到了最佳效率指數，但不易擴展至更高收斂階［11］。理論上相應的區間方法可提高計算穩定性，但將耗費更多的計算時間［12-13］。基于笛卡爾法則與Sturm 定理的多項式求根方法雖具有較好的計算穩定性，但收斂速度較低。

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

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

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

1 RBM

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

1.1 思路及算法原理

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

其中，t0=a和t1=b。

則Ri(t)和f(t)對應的無窮級數在區間[a，b]內收斂，即

通過漸進式插值更多的點可得到更小的逼近誤差。

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

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

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

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

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

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

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

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

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

本文方法無需計算gi(s)的表達式，計算效率比現有的漸進法更高，且兼具漸進法良好的計算穩定性。

1.2 算法概述

引入文獻［24］中的定理3.5.1，即

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

設t0=a，t1=b，可驗證

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

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

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

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

式（10）的求根步驟為

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

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

輸出實根t*的近似值。

步驟 1令t2，i=2。

步驟2若轉步驟4。否則，轉步驟3。

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

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

1.3 收斂階

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

定理2設則有

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

證明由式（13），有

由式（15），當i=2 時，有

由式（7），可得

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

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

證畢。

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

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

即

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

式（20）等價于

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

證畢。

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

1.4 RBM 示例

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

例 1令f1(t)=(t?0.25)(2 ?t)(t+5)2，t∈[0，1]，在區間［0，1］內有一個單根t*=0.25。令可得線性函數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。本例中，精度設置為小數點后100。

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

2 不同方法的定性比較

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

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

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

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

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

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

表4 收斂階為16 時對應的FE 次數的比較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 個實例，對不同方法進行了比較。總體來說，M1具有與裁剪法相當的計算穩定性和更高的計算效率、比類牛頓法更高的計算穩定性和更高的計算效率、與漸進法相當的收斂階和更高的計算效率。

3.1 M1與裁剪法的比較

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

例2測試以下4 個多項式函數：

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

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

表6 例2 中4 個多項式函數在不同方法下的誤差和計算時間比較Table 6 Comparisons of the error and computational time of 4 polynomial functions in different methods of example 2

3.2 漸進法與類牛頓法的比較

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

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

表7 4 個多項式函數在不同方法下的裁剪步數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）

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

表8 3 個非多項式函數在不同方法下的逼近誤差比較Table 8 Comparisons of approximation errors of 3 non-polynomial functions in different methods

表9 3 個非多項式函數在不同方法下運行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）考慮如何獲取包含一個根的區間。設在初始值x0下，x1通過牛頓法或其他迭代方法得到。若2|x1?t*|＜|x0?t*|，t*是f（t）的根，由定理3，t*被x0和2x1?x0所包圍。因此，通過選取適當的初始值x0可得到一個包含實根的小區間。

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

例4用RBM 計算Wilkinson 多項式

在區間［0，25］內的實根。該Wilkinson 多項式有20個實根Wi，i=1，2，…，20（參見文獻［15］中的例6）。首先，計算相應控制多邊形的零點，即｛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｝，并將區間［0，25］劃分為20 個子區間，其中有16 個子區間包含W（t）的1 個或2 個根。本文選取其中的3 個子區間Λ1=［0.27，1.55］、Λ2=［2.83，4.11］、Λ3=［16.81，18.07］用以說明更多細節。此3 個子區間分別包含1 個、2 個、2 個根。可將RBM 直接應用于Λ1，因為W（t）在其2 個端點的函數值異號；對于Λ2和Λ3，可選擇區間中點作為分界點，將每個區間拆分為2 個子區間，每個子區間都包含1 個根。然后，將RBM 應用于子區間，由注1的方法獲得優化區間。本例中，區間［0.27，1.55］、［2.83，4.11］、［16.81，18.07］對應的優化區間為Γ1=［0.99，1.0034］、Γ2=［2.83，3.004］、Γ3=［17.988，18.07］。如表10 所示，在誤差ei映射到重新參數化函數φi(t)的情況下，RBM 的效果較好，表明用RBM 處理多根情形也具有較好的計算穩定性。

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

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

4 結 論

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

未來的工作，包括進一步擴展RBM，選取合適的初始值，并將其應用于非線性方程或方程組的求根。

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