單位圓到任意多邊形區(qū)域的Schwarz Christoffel變換數(shù)值解法

崔建斌，姬安召，王玉風(fēng)，于江濤，周華龍

(1.隴東學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，甘肅 慶陽 745000； 2.隴東學(xué)院 能源工程學(xué)院，甘肅 慶陽 745000)

單位圓到任意多邊形區(qū)域的Schwarz Christoffel變換數(shù)值解法

崔建斌1，姬安召2，王玉風(fēng)2，于江濤2，周華龍2

(1.隴東學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，甘肅 慶陽 745000； 2.隴東學(xué)院 能源工程學(xué)院，甘肅 慶陽 745000)

Schwarz Christoffel變換技術(shù)在處理某些工程問題時具有重要作用. 從黎曼存在定理出發(fā)，建立了單位圓到任意多邊形區(qū)域的映射函數(shù)Schwarz Christoffel變換模型，采用LevenbergMarquardt算法求解含約束條件的非線性映射函數(shù)Schwarz Christoffel變換模型參數(shù)系統(tǒng).針對映射函數(shù)中出現(xiàn)的奇異積分問題，對映射函數(shù)進行2次參數(shù)變換，將其化為高斯雅克比型積分，以積分路徑中的奇異點為界，縮短積分路徑，對子路徑采用修正高斯積分方法進行計算.通過指數(shù)變換、連乘變換和累加變換，使任意初值問題均可進行迭代計算并滿足初值的約束條件.提出以邊長絕對誤差和頂點絕對誤差為迭代計算的收斂條件，并保證了映射函數(shù)的精度.給出了11頂點多邊形區(qū)域映射函數(shù)的求解算例，4種方案的計算結(jié)果表明，Schwarz Christoffel變換數(shù)值解法操作簡單、精度高、收斂快.

修正高斯雅克比型積分；單位圓；初值變換；LevenbergMarquardt算法；Schwarz Christoffel變換

0 引 言

Schwarz Christoffel變換(以下簡稱SC變換)模型在解決實際工程問題時具有重要作用.SC變換能將二維空間上邊界復(fù)雜的幾何體映射到另一個二維空間的形狀簡單的幾何體上，從而簡化對復(fù)雜邊界工程問題的處理.SC變換在油氣地下滲流力學(xué)、巖土力學(xué)、流體力學(xué)和電磁學(xué)等研究領(lǐng)域有著廣泛的應(yīng)用.目前，國內(nèi)已有一些關(guān)于SC變換數(shù)值計算方法的研究成果，王剛等[12]采用牛頓拉夫森迭代法導(dǎo)出了從多角形區(qū)域到上半平面映射和槽型區(qū)域映射的SC變換求解方法，但未涉及多角形區(qū)域到單位圓的計算.祝江鴻等[34]采用SC洛朗級數(shù)模型導(dǎo)出了地下開挖隧道斷面到單位圓映射的計算方法，但對于復(fù)雜開挖斷面，因級數(shù)構(gòu)成項較多而求解復(fù)雜.文獻[57]均以SC級數(shù)模型建立多角形區(qū)域到單位圓映射計算，其級數(shù)模型構(gòu)成復(fù)雜、計算量大，精度難以控制.文獻[8]研究了多角形區(qū)域到單位圓的映射，但由于奇異積分計算精度較低，要得到高精度的結(jié)果則比較困難.文獻[9]研究了多角形區(qū)域到上半平面SC逆變換的計算方法，同樣，計算結(jié)果精確但耗時較長.文獻[10]主要探討了由多角形區(qū)域到上半平面SC變換的數(shù)值解法.筆者借鑒前人的研究成果，通過研究多角形區(qū)域到單位圓映射SC變換模型，經(jīng)2次參數(shù)變換，將沿單位圓復(fù)積分轉(zhuǎn)換為實積分，并建立了其與高斯雅克比型積分的關(guān)系，提出可通過搜尋奇異點縮短積分路徑，并在子路徑中采用修正高斯雅克比積分.采用上述方法解決了計算過程中出現(xiàn)的奇異積分問題，亦提高了積分精度.結(jié)合LevenbergMarquardt最優(yōu)化算法求解非線性積分方程組，尋求由單位圓到多角形區(qū)域SC變換的精確高效的數(shù)值解法.

圖1 單位圓映射到任意多角形區(qū)域變換示意圖Fig.1 Transformation diagram of unit circleto the polygonal region

1 基本模型

在復(fù)平面W上有N(N≥3)邊形，其頂點和內(nèi)角分別為Wj和παj(j=1,2,…,N).將單位圓上的點映射在W平面上的SC變換公式為

(1)

其中，K為伸縮系數(shù)，C為變換中心，παj為W平面多角形內(nèi)角，zj為單位圓上的點.

但在實際工程問題研究中，一般多角形區(qū)域是已知的，需求解與多邊形頂點對應(yīng)的單位圓上的映射點.據(jù)黎曼原理，要確定式(1)，則參數(shù)zj(j=1,2,…,N)中有3個點必須選定.不妨假定zN-2=-1+0i，zN-1=0-1i，zN=1+0i，亦可根據(jù)實際工程問題選擇其他點.在上述假定條件下，式(1)中zj(j=1,2,…,N)滿足：

0…

(2)

由式(1)可得多邊形頂點Wi(i=1,2,…,N)為

(3)

由式(3)則有

(4)

用相鄰兩點之間的邊長比值可消去伸縮系數(shù)，從而減少未知量的求解，由式(4)可得

(5)

令

(6)

則式(5)可表示為

(7)

(8)

通過求解非線性系統(tǒng)積分方程組式(8)，可得未知參數(shù)zj(j=1,2,…,N-3).

由式(4)可求得伸縮系數(shù)K：

(9)

由式(3)，結(jié)合未知參數(shù)zj(j=1,2,…,N-3)與伸縮系數(shù)K的求解值，可得變換中心C：

(10)

2 Schwarz Christoffel積分

2.1 高斯雅克比型積分

根據(jù)以上分析，為了求解SC變換參數(shù)問題式(8)～(10)，必須計算式(6)的積分.分析實際問題可知，式(6)積分路徑在單位圓的圓弧上，積分的起點zi和終點zi+1為奇點.為了求解式(6)的奇異積分，對其作1次參數(shù)變換，因zi在單位圓上，令zi=eiθi，其中θi為復(fù)數(shù)zi的輻角主值，代入式(6)可得

(11)

為了計算式(11)的奇異積分，對式(11)進行2次參數(shù)變換：

令

(12)

令

(13)

其中，當j=i時，該項的模可近似表示為

則式(11)可表示為

(14)

式(14)即為高斯雅克比型積分，但需注意其積分的起點和終點都存在奇異點，這里只給出了積分起點奇異點的處理方法，為了處理積分終點的奇異點，只需將積分區(qū)間劃分為若干個子路徑，因子路徑內(nèi)不含奇異點.具體劃分方法參見2.2節(jié)積分路徑奇異點的確定.根據(jù)文獻[11]，式(14)可表示為：

(15)

(16)

第2步 尋找子區(qū)間的奇異值，由式(16)可得奇異值的條件為

(j=1,2,…,start-1,start,…,N-1;start=i).

第5步 若dist≥1，采用校正后的零點和權(quán)值進行積分.

2.3 校正高斯雅克比型積分零點與權(quán)值

結(jié)合式(12)與(14)，校正后的零點為

(17)

考慮高斯雅克比型積分式(14)，校正后的權(quán)值

(18)

3 非線性系統(tǒng)求解

對非線性系統(tǒng)式(8)的求解，還得考慮其約束條件0

由式(8)與式(2)，得到約束條件下的非線性系統(tǒng)求解問題：

(19)

3.1 約束條件的變換

在采用迭代法求解式(19)的過程中，未知量zj(j=1,2,…,N-3)在迭代過程中有可能出現(xiàn)不滿足式(2)的情況，為此對約束條件做適當?shù)淖儞Q[1415]，使其在求解過程中始終滿足約束條件.變換過程如下：

第5步 計算迭代值向量Y=[y1,y2,…,yN3]T，yj=eiθj.

通過上述變化后，只要給定任意初值向量X0，迭代值向量Y必然滿足式(2).

3.2LevenbergMarquardt算法[16]的參數(shù)優(yōu)化及數(shù)值計算方案

圖2 LevenbergMarquardt算法中參數(shù)、迭代次數(shù)與誤差關(guān)系曲線Fig.2 The relationship curves of number of iterations，error and parameters in LevenbergMarquardt algorithm

在采用LevenbergMarquardt算法求解式(19)時，由ρ與σ取值與迭代次數(shù)、誤差的關(guān)系曲線(見圖(2))可以看出，在迭代初期，ρ和σ的取值對絕對誤差的影響較小，收斂速度較慢，當?shù)螖?shù)大于60時，ρ和σ的取值對絕對誤差的影響較大，收斂速度差異較大.通過反復(fù)迭代計算，得到當ρ=0.2，σ=0.5時，收斂速度最快，為此推薦使用ρ=0.2，σ=0.5.

表1 非線性積分方程組數(shù)值計算方案Table 1 Numerical scheme for nonlinear integral equations

通過上述分析，求解非線性積分方程組式(19)，初值必須滿足式(2)，對初值的處理主要采用上述約束條件的變換方法：包含連乘變換和無連乘變換，所以對初值處理有2種選擇方案.對于出現(xiàn)的奇異積分，可直接采用式(14)修正高斯雅克比積分或式(16)子路徑修正高斯雅克比積分進行計算.這樣共有4種計算方案，其計算結(jié)果見表3.

4 精度評定與算例分析

i=1,2,…,N-2,

(20)

i=1,2,…,N-2.

(21)

平面封閉區(qū)域多邊形如圖3所示，這里共有11個頂點封閉多邊形.表2是SC變換參數(shù)的計算結(jié)果，由于篇幅限制，僅列出②③、①③和②④3種計算方案的結(jié)果.從表2可以看出，前2種方案的計算結(jié)果很接近，說明初值變換選擇的2種方案均可行，但奇異積分處理方法的選擇對結(jié)果影響很大.

圖3 多邊形平面圖Fig.3 Polygon plane graph

方案1(②③)方案2(①③)方案3(②④)z1-0.9730436674+i*0.2306209473-0.9730436674+i*0.2306209474-0.989757743+i*0.142757174z2-0.9759654834+i*0.2179251598-0.9759654834+i*0.2179251599-0.990369302+i*0.138450875z3-0.9849015986+i*0.1731151093-0.9849015986+i*0.1731151093-0.992766049+i*0.120064866z4-0.9857537888+i*0.1681947321-0.9857537888+i*0.1681947321-0.993105233+i*0.117226258z5-0.9857746537+i*0.1680724015-0.9857746537+i*0.1680724015-0.99311812+i*0.117117035z6-0.9858563977+i*0.1675922523-0.9858563977+i*0.1675922523-0.993169998+i*0.116676285z7-0.9858588436+i*0.1675778641-0.9858588436+i*0.1675778641-0.99317136+i*0.116664692z8-0.9858878094+i*0.1674073695-0.9858878094+i*0.1674073695-0.993184249+i*0.11655491k-0.29320062207+i*0.48877316531-0.29320062211+i*0.48877316531-0.23382918823+i*0.43935654317c-1.5840515459+i*1.4956209022-1.5840515460+i*1.4956209021-1.6837743041+i*1.6952215685

表3 從單位圓到任意多邊形SC變換4種計算方案結(jié)果的對比Table 3 Comparison of four calculation schemes for SC transformation from unit circle to arbitrary polygon

通過Matlab軟件編寫SC變換程序，在PC機上進行計算，由表3可知，若直接采用修正高斯雅克比積分方法，無論如何選擇初值變換方案，邊長和頂點都不能達到精度要求(即方案3和4)，而且計算時間長.方案2與方案1相比，若初值不采用連乘變換，上述變換結(jié)果同樣滿足式(2)，但需要多次迭代邊長和頂點才能達到精度要求，所以初值的連乘變換能夠提高計算速度，故推薦方案1.

5 結(jié) 論

5.1 單位圓到多邊形區(qū)域共形映射函數(shù)SC變換模型求解的難點，在于計算過程中出現(xiàn)的奇異積分方程的處理和初值必須滿足一定的約束條件.從常規(guī)高斯雅克比型積分出發(fā)，對SC變換模型進行2次參數(shù)變換，通過變換積分方程的形式，以積分路徑中的奇異點為界，縮短積分路徑，在子路徑中采用修正高斯積分方法可有效克服奇異積分計算的困難.對于初值的約束條件，通過指數(shù)變換、連乘變換和累加變換巧妙處理迭代值，使其對任意初值均可進行迭代并滿足初值的約束條件.

5.2 通過對4種方案的計算分析知，初值和積分方式的選擇決定了SC變換模型系數(shù)迭代計算的快速性和有效性，提出的子路徑修正高斯雅克比積分和初值變換方法可有效解決該問題；以邊長的絕對誤差和頂點絕對誤差作為迭代計算的收斂條件，保證了計算映射函數(shù)的映射精度.

5.3 在求解SC變換模型時，初值的連乘變換能夠較大幅度減少非線性系統(tǒng)求解過程中的迭代次數(shù).

5.4 LevenbergMarquardt算法對映射函數(shù)SC變換模型非線性系統(tǒng)是可行的，當參數(shù)ρ與σ取值適當時，LevenbergMarquardt算法收斂較快.

5.5 11個頂點多邊形映射函數(shù)的求解算例表明，本文給出的從單位圓到多邊形區(qū)域共形映射函數(shù)SC變換模型求解算法具有一定的可操作性.

[1] 王剛，許漢珍，顧王明，等.數(shù)值許瓦爾茲克力斯托夫變換與數(shù)值高斯雅可比型積分[J].海軍工程學(xué)院學(xué)報，1994，1(2)：2534. WANG G，XU H Z，GU W M，et al. Numerical schwarzchristoffel transformation and numerical gaussjacobi quadrature[J]. Journal of Naval Academy of Engineering，1994，1(2)：2534.

[2] 王剛，陸小剛，顧王明.槽形內(nèi)域中的數(shù)值許瓦爾茲克力斯托夫保角變換[J].海軍工程學(xué)院學(xué)報，1995，1(4)：1623. WANG G，LU X G，GU W M. Numerical schwarzchristoffel conformal mapping in channel region[J].Journal of Naval Academy of Engineering，1995，1(4)：1623.

[3] 祝江鴻.隧洞圍巖應(yīng)力復(fù)變函數(shù)分析法中的解析函數(shù)求解[J].應(yīng)用數(shù)學(xué)和力學(xué)，2013，34(4)：345354. ZHU J H. Analytic functions in stress analysis of the surrounding rock for caverns with the complex variable theory[J]. Applied Mathematics and Mechanics，2013，34(4)：345354.

[4] 祝江鴻，楊建輝，施高萍，等.單位圓外域到任意開挖斷面隧洞外域共形映射的計算方法[J].巖土力學(xué)，2014，35(1)：175183. ZHU J H，YANG J H，SHI G P，et al. Calculating method for conformal mapping from exterior of unit circle to exterior of cavern with arbitrary excavation crosssection[J]. Rock and Soil Mechanics，2014，35(1)：175183.

[5] 皇甫鵬鵬，伍法權(quán)，郭松峰.基于邊界點搜索的洞室外域映射函數(shù)求解法[J].巖石力學(xué)，2011，32(5)：14181424. HUANGFU P P，WU F Q，GUO S F. A new method for calculating mapping function of external area of cavern with arbitrary shape based on searching points on boundary[J]. Rock and Soil Mechanics，2011，32(5)：14181424.

[6] 朱大勇，錢七虎，周早生.復(fù)雜形狀洞室映射函數(shù)的新解法[J].巖石力學(xué)與工程學(xué)報，1999，18(3)：279282. ZHU D Y，QIAN Q H，ZHOU Z S. New method for calculating mapping function of opening with complex shape[J]. Chinese Journal of Rock Mechanics and Engineering，1999，18(3)：279282.

[7] 王潤富.一種保角映射法及其微機實現(xiàn)[J].河海大學(xué)學(xué)報，1991，19(1)：8689. WANG R F. A method of conformal mapping and its computer implementation[J]. Journal of Hohai University，1991，19(1):8689.

[8] HU C H. A software package for computing schwarzchristoffel conformal transformation for doubly connected polygonal regions[J]. ACM Transactions on Mathematical Software，1998，24(3)：317333.

[9] COSTAMAGNA E. Numerical inversion of the schwarzchristoffel conformal transformation： Stripline case studies[J]. Microwave and Optical Technology Letters，2001，28(3)：179183.

[10] 崔建斌，姬安召，魯洪江，等.Schwarz Christoffel變換數(shù)值解法[J].山東大學(xué)學(xué)報：理學(xué)版，2016，51(4):104111. CUI J B，JI A Z，LU H J，et al. Numerical solution of Schwarz Christoffel transform[J]. Journal of Shandong University： Natural Science，2015，50(12):104111.

[11] COSTAMAGNA E. A new approach to standard schwarzchristoffel formula calculations[J]. Microwave and Optical Technology Letters，2002，32(3)：196199.

[12] DRISCOLL T A. Improvements to the schwarzchristoffel toolbox for matlab[J]. ACM Transactions on Mathematical Software，2005，31(2)：239251.

[13] HOUGH D M. Asymptotic gauss jacobi quadrature error estimation for schwarzchristoffel integrals[J]. Journal of Approximation Theory，2007，146(2)：157173.

[14] NATARAJAN S， BORDAS S P A， MAHAPATRA R D. Numerical integration over arbitrary polygonal domains based on schwarzchristoffel conformal mapping[J]. International Journal for Numerical Methods in Engineering，2009，80(1)：103134.

[15] CROWDY D. The schwarzchristoffel mapping to bounded multiply connected polygonal domains[J]. Proceedings of the Royal Society，2005，146(2061)：26532678.

[16] 劉浩.大規(guī)模非線性方程組和無約束優(yōu)化方法研究[D].南京：南京航空航天大學(xué)，2008. LIU H. Research on Methods for Largescale Nonlinear Equations and Unconstrained Optimization[D]. Nanjing： Nanjing University of Aeronautics and Astronautics，2008.

CUI Jianbin1，JI Anzhao2，WANG Yufeng2， YU Jiangtao2， ZHOU Hualong2

(1.MathematicsandStatisticsInstitute，LongdongUniversity，Qingyang745000，GusuProvince,China；2.EnergyEngineeringInstitute，LongdongUniversity，Qingyang745000，GusuProvince,China)

Schwarz Christoffel transformation technique has an important role in dealing with engineering problem. Based on the Riemann’s existence theorem, a transform model of Schwarz Christoffel from unit circle to arbitrary polygon area was established. LevenbergMarquardt algorithm was used to solve parameters of the nonlinear mapping function of Schwarz Christoffel with constraint conditions. As to the mapping function in the singular integral problem, the mapping function was twice transformed to Gauss Jacobi integral, then the singular point in the integral path was taken as the boundary; The length of integral path was narrowed, and the modified Gaussian integral was used to calculate the sub path. By the exponential transform, multiplicative transform and accumulation transform process, the arbitrary initial value could be iteratively calculated and the results satisfied the constraints of the initial value. The convergence conditions of the iterative copulation based on the absolute errors of length and absolute error of vertex were put forward to ensure the precision of the mapping function. The solution of the mapping function with 11 vertex polygon region was calculated by four calculation schemes. Results of the four schemes showed that the numerical solution of Schwarz Christoffel transform is simple with high precision and fast convergence.

modified Gauss Jacobi type integral; unit circle; initial value transformation; LevenbergMarquardt algorithm; Schwarz Christoffel transformation

20160301.

甘肅省科技計劃資助項目(1606RJZM092,1606RJYM259,1506RJYM324).

崔建斌(1972-),ORCID:http://orcid.org/0000000266933415,男,碩士,副教授,主要從事數(shù)值分析與數(shù)據(jù)挖掘研究，Email:cuijb0658@163.com.

10.3785/j.issn.10089497.2017.02.007

O 241

A

1008-9497(2017)02-161-07

Numerical solution method for Schwarz Christoffel transformation from unit circle to arbitrary polygon area. Journal of Zhejiang University(Science Edition), 2017,44(2):161167