預處理GPBi-CG 算法的數值保角變換計算法

石允龍，呂毅斌，王櫻子，伍 康

（1.昆明理工大學 理學院；2.昆明理工大學 計算中心，云南 昆明 650500）

0 引言

保角變換是研究復變函數的重要方法之一，在物理學和工學中有著廣泛應用，如流體力學、光學、空氣動力學、圖像處理等很多領域的實際問題中都有保角變換的身影［1-2］。1952 年，Nehari［3］給出了5 種狹縫域作為多連通域保角變換的重要規(guī)范域，分別是平行狹縫域、弧形狹縫域、徑向狹縫域、帶有弧形狹縫的圓盤域及帶有弧形狹縫的圓環(huán)域，任何多連通域都可映射到這5 種規(guī)范域上；20 世紀80 年代，日本學者Amano 對模擬電荷法和數值保角變換作了大量研究，提出基于模擬電荷法的數值保角變換計算法（Amano 法）［4-6］；2007 年，Nasser 將求解保角變換問題轉化為黎曼希爾伯特問題，并用帶有廣義Neumann kernel 的邊界積分方程求解該問題［7-8］。國外研究學者們提出了許多保角變換的數值方法，但針對提高數值保角變換精度的研究較少。

目前國內許多學者針對不同類型問題域的數值保角變換，提出一些提高模擬電荷法計算精度的方法［9-11］，但針對提高有界多連通區(qū)域數值保角變換精度的研究很少。因此，本文通過基于模擬電荷法的數值保角變換法將有界多連通區(qū)域映射到有相同圓心的帶有弧形狹縫的圓環(huán)，并給出提高該問題域數值保角變換計算精度的算法。

本文首先介紹基于模擬電荷法的有界多連通區(qū)域數值保角變換算法，由于該方法產生的約束方程組有很強的病態(tài)性，因此根據預處理思想［12］，采用1-范數均衡法處理約束方程組，使矩陣所有行或列有相同的模，降低矩陣的條件數，從而改善矩陣的病態(tài)性［13］；然后，通過基于雙共軛梯度的廣義乘積型算法（GPBi-CG）［14］求解預處理后的約束方程組，得到更高精度的方程組數值解，即保角變換的模擬電荷與變換半徑，進而得到近似保角變換函數；最后，通過數值實驗驗證了本文提出的數值保角變換新算法的優(yōu)越性。

1 相關研究

1.1 有界多連通區(qū)域保角變換函數

本節(jié)主要闡述基于模擬電荷法的有界多連通區(qū)域的數值保角變換計算法［15-16］。圖1 是多連通域保角變換的重要規(guī)范域之一，即將有界多連通區(qū)域D映射到帶有弧形狹縫的圓環(huán)域E上。其中，問題域D由z平面上的閉合Jordan曲線C1,C2,…,Cn構成，且閉合曲線C2,…,Cn被閉合曲線C1包圍；目標域E由w平面上的單位圓S1、同心圓S2以及與圓環(huán)有相同圓心的弧形狹縫S3,…,Sn構成，即帶有圓弧狹縫的單位圓環(huán)。

Fig.1 Numerical conformal mappings of bounded multiply connected domains by the charge simulated method（“+”denoted the position of charge points）圖1 基于模擬電荷法的有界多連通區(qū)域數值保角變換（“+”表示模擬電荷點位置）

由Riemann 映射定理可知，在滿足正則化條件f(u)=0,f(v)=1，u和v為正則化點，且點u在閉合曲線C2內部，點v是曲線C1上的點時，從多連通域D到規(guī)范狹縫域E存在唯一的保角映射函數w=f(z)，將曲線C1,C2,…,Cn分別映射成圓或圓弧狹縫S1,S2,…,Sn。

根據保角變換的幾何性質，可得到邊界約束條件：

其中，r1,r2,…,rn是圓與圓弧狹縫的半徑，且r1=1。

在不失一般性的情況下，不防令z=0 在曲線C2內部，且f(0)=0，利用調和函數g(z)及其共軛調和函數h(z)，將有界多連通區(qū)域D映射到典型狹縫域E的保角變換函數表示為：

再根據正則化條件中f(v)=1，得：

由邊界約束條件，得：

根據保角變換函數存在的唯一性，本文將求解映射函數f(z)問題轉化為求解調和函數g(z)與h(z)的問題。

1.2 基于模擬電荷法的有界多連通區(qū)域數值保角變換

模擬電荷法通過復對數函數的線性組合近似調和函數g(z)與h(z)，即：

其中，Q0為復常數，Qli為未知電荷。如圖1 所示，模擬電荷點ζli在區(qū)域D外部，即N1個模擬電荷點分布在邊界C1外部，剩余的N2,…,Nn個模擬電荷點分別分布在邊界C2,…,Cn內 部。

因為G(z)、H(z)分別是g(z)與h(z)的近似映射函數，所以由式（3）、式（4）可知：

其中，Ri是ri的近似，zmk是邊界曲線Cm上的約束點。通過式（7）計算出Q0并代入式（5），得：

聯立式（6）與式（8），有：

由于保角變換函數具有單值性與伸縮率不變性［17-18］，近似映射函數也應滿足相應性質，從而有：

通過解式（9）、（10）與（11）聯立所得的N1+N2+…+Nn+n階大型線性方程組，得到電荷Qli,l=1,…,n,i=1,…,Nl和變換半徑Rm,m=1,…,n的值，從而求得保角變換近似函數：

2 本文方法

2.1 1-范數均衡法

將由式（9）、（10）與（11）聯立所得的N(=N1+N2+…+Nn+n)階大型線性方程組寫成Ax=b的形式。因為該方程為非對稱病態(tài)方程，所以本文首先基于1-范數均衡法降低方程組的條件數，再采用GPBi-CG 迭代算法求解預處理后的線性方程組。

對線性方程組Ax=b兩邊同乘以矩陣U，且令x=Bx*，得：

式中，U、B為對角矩陣，且矩陣U的對角元素為uii=,i=1,2,…,N，矩陣B的對角元素為bjj=,j=1,2,…,N，其中S與T為任意常數。矩陣A左乘U后可使每行的1-范數相同，矩陣UA再右乘B后每列的1-范數相同。行列的1-范數均衡可連續(xù)重復下去，一般經過幾次重復后，系數矩陣的條件數則不再改變，方程的病態(tài)性也不再得到改善，此時可停止重復，得到最后經過預處理的線性方程：

2.2 基于預處理的GPBi-CG 算法

將式（14）改寫為：

其中，A*=Uk…U2U1AB1B2…Bk,b*=Uk…U2U1b。

利用GPBi-CG 算法求解得到x*，再由x=B1B2…Bk x*得到原方程組的解x，該方法稱為基于預處理的GPBi-CG算法。GPBi-CG 算法不僅避免了雙共軛梯度法（Bi-CG）中使用轉置矩陣生成Krylov 子空間，從而需要進行的轉置矩陣與向量乘法運算，而且提高了收斂速度。

根據文獻［19］、［20］，可得到求解線性方程組的Bi-CG算法。算法具體步驟如下：

根據文獻［14］可知，GPBi-CG 算法是在Bi-CG 算法基礎上，通過使用Bi-CG 算法的殘差多項式Rn與其它具有標準三項遞推關系的多項式Hn乘積，定義新算法的殘差多項式，并且由Rn的三項遞推關系：

給出Hn的三項遞推關系：

其中，參數ηn與參數ζn可由殘差的最小二范數來確定，即：

由這兩組多項式的乘積得到相應的多項式遞推關系，從而得到新算法，即GPBi-CG 算法的遞推公式。

最后本文參考文獻［21］的思想，將1-范數均衡預處理與GPBi-CG 算法相結合，建立求解模擬電荷與變換半徑的預處理GPBi-CG 算法。算法具體步驟如下：

3 數值實驗

本文通過以下兩個基于模擬電荷法的有界多連通區(qū)域保角變換數值實驗檢驗本文方法，即基于預處理的GP?Bi-CG 算法相比Gauss-Seidle 法［22-23］與Amano 法的優(yōu)越性，并使用MATLAB 軟件編寫程序。

基于模擬電荷法的數值保角變換根據拉普拉斯方程最大值原理在邊界上進行誤差評價，誤差計算公式為：

例1：多連通區(qū)域D以正方形C1為外邊界，以正方形C2、C3為內邊界，其邊界曲線方程為：

約束點為：

模擬電荷點為：

剩余約束點與模擬電荷點按照軸對稱配置即可。

圖2 為當Nl=56(l=1,2,3)時問題域邊界及模擬電荷點位置圖，圖3 是圖2 基于本文方法的保角變換圖。

表1 給出了取不同數量（32、56、80、104、120）的模擬電荷點時，本文方法與Gauss-Seidel 法及Amano 法的對比，包括誤差、迭代次數及條件數情況。從表1 可以看出，本文方法經過預處理后，方程組的條件數降低到傳統(tǒng)方法條件數的左右。對于內、外邊界Cl，本文方法的誤差明顯小于Gauss-Seidel 法，且迭代次數遠少于Gauss-Seidel 法；對于外邊界C1，本文方法的誤差明顯小于Amano 法，充分驗證了本文方法的優(yōu)越性。

Fig.2 Boundary and position of charge points（example 1）圖2 邊界及模擬電荷點位置（例1）

Fig.3 Fig.2 Transformation diagram based on the proposed method（example 1）圖3 圖2 基于本文方法的保角變換圖（例1）

Table 1 Numerical conformal mapping error，iterations and condi?tion number（example 1）表1 數值保角變換誤差、迭代次數與條件數（例1）

圖4 給出了3 種方法更直觀的誤差曲線圖。結合表1與圖4 可以看出，當各邊界模擬電荷點數N=104 時，邊界曲線C1基于本文方法的誤差為E1=2.90 × 10-4，迭代次數為1 694 次；Gauss-Seidel 法的誤差為E1=1.38 × 10-2，迭代次數為2 963 次；Amano 法的誤差為E1=6.40 × 10-3，且本文算法隨著模擬電荷點的增加，誤差穩(wěn)定下降。

Fig.4 Error curve of conformal mapping（example 1）圖4 保角變換誤差曲線（例1）

為了進一步驗證本文算法的有效性，下面以菱形為外邊界的有界多連通區(qū)域的保角變換情況為例進行說明。

例2：多連通區(qū)域D以菱形C1為外邊界，以正方形C2和C3為內邊界的邊界曲線方程為：

約束點位置為：

模擬電荷點位置為：

內邊界C2、C3剩余約束點與模擬電荷點同樣按照軸對稱配置即可。

同樣的，圖5 為當Nl=56(l=1,2,3)時問題域邊界及模擬電荷點分布圖，圖6 是圖5 基于本文方法的保角變換圖。

表2 同樣給出了取不同數量（32、56、80、104、120）的模擬電荷點時，本文方法與Gauss-Seidel 法及Amano 法的對比，包括誤差、迭代次數與方程組條件數情況。由表2 可知，本文方法經過預處理后，方程組的條件數降低到傳統(tǒng)方法條件數的左右，且對于內、外邊界Cl，本文方法的誤差明顯小于Gauss-Seidel 法，且迭代次數遠少于Gauss-Seidel 法；對于外邊界C1，本文方法的誤差明顯小于Amano法，再次驗證了本文方法的優(yōu)越性。

Fig.5 Boundary and position of charge points（example 2）圖5 邊界及模擬電荷點位置（例2）

Fig.6 Fig.5 Transformation diagram based on the proposed method（example 2）圖6 圖5 基于本文方法的保角變換圖（例2）

Table 2 Numerical conformal mapping error，iterations and condi?tion number（example 2）表2 數值保角變換誤差、迭代次數與條件數（例2）

同樣的，本文通過圖7 給出了3 種方法更直觀的誤差曲線圖。結合表2 與圖7 可以看出，當各邊界模擬電荷點數N=104 時，邊界曲線C1基于本文算法的誤差為E1=8.50×10-4，迭代次數為44 次；Gauss-Seidel 法的誤差為E1=4.50×10-2，迭代次數 為3 300 次；Amano 法的 誤差為E1=3.31×10-2，且本文算法隨著模擬電荷點的增加，誤差穩(wěn)定下降，再次驗證了本文算法的優(yōu)越性。

Fig.7 Error curve of conformal mapping（example 2）圖7 保角變換誤差曲線（例2）

4 結語

本文提出基于1-范數均衡預處理GPBi-CG 算法的有界多連通區(qū)域數值保角變換計算法，并通過多個數值實驗驗證了該方法在算法精度、穩(wěn)定性與迭代次數上都有較好表現，對于求解該類型問題域的近似保角變換函數具有明顯的優(yōu)越性。保角變換在圖像處理中有著廣泛應用，未來可進行一些具體應用的研究，比如本文研究結果可應用于將不規(guī)則或破損圖片映射形成環(huán)形，并進一步將這些圖片信息刻錄到光學媒介、CD 或DVD 上長久保存下來。但不足的是，本文算法在剩余的4 種規(guī)范域、單連通域和雙連通域等數值保角變換問題中的應用效果還有待驗證，且本文是根據經驗選取模擬電荷點，在未來研究中也可考慮給出更科學的選取方法。