呂毅斌,王 堅(jiān),王櫻子,吳 爽
(1.昆明理工大學(xué)理學(xué)院,云南昆明 650500)
(2.昆明理工大學(xué)計(jì)算中心,云南昆明 650500)
保角變換是復(fù)變函數(shù)中非常重要的理論之一,廣泛應(yīng)用于物理學(xué)和工學(xué)領(lǐng)域.特別是在電磁理論、膜和板的振動(dòng)、彈性理論、熱傳輸、流體力學(xué)等方面有很多應(yīng)用[1?5].通常將保角變換的求解方法分為解析法和數(shù)值法.解析法只能在一些特殊區(qū)域給出變換函數(shù)表達(dá)式,對(duì)于復(fù)雜區(qū)域問題沒有解決方法.因此,對(duì)于很多實(shí)際中的復(fù)雜問題必須采用數(shù)值法求解保角變換函數(shù).數(shù)值保角變換計(jì)算法的主要有:積分方程式法[6]、正交多項(xiàng)式法[7?8]和有限差分法[9?10]等.德國人Steinbigler[11]首次提出用若干個(gè)虛設(shè)電荷來模擬電極表面上電荷分布的電場(chǎng)計(jì)算法,形成了模擬電荷法的基本思想;日本的天野要等數(shù)學(xué)學(xué)者從20世紀(jì)80年代開始對(duì)模擬電荷法和數(shù)值保角變換作了大量研究工作,并提出了基于模擬電荷法的數(shù)值保角變換計(jì)算法(天野法)[12?16].天野法適用于單連通區(qū)域及多連通區(qū)域的數(shù)值保角變換問題[15,17,18,19,20].
單連通區(qū)域的數(shù)值保角變換,分為內(nèi)部數(shù)值保角變換和外部數(shù)值保角變換.本文通過對(duì)基于模擬電荷法的雙方向的內(nèi)部數(shù)值保角變換計(jì)算法[13,15,16]的研究,提出了外部數(shù)值保角逆變換計(jì)算法.該方法的原理是基于模擬電荷法來求Laplace方程的Dirichlet問題的解,并通過預(yù)先建立的邊界對(duì)應(yīng)關(guān)系構(gòu)造從標(biāo)準(zhǔn)區(qū)域到問題區(qū)域的近似保角逆變換函數(shù),誤差用正則函數(shù)的最大值原理進(jìn)行評(píng)價(jià).計(jì)算數(shù)值保角逆變換最主要的就是確定模擬電荷的位置和數(shù)量,以及約束方程組的求解.
LSQR方法[21?24]是Paige和Saunders提出的一種適用于求解系數(shù)矩陣為大型、稀疏矩陣線性方程組的方法.LSQR方法求解的思路是把任意系數(shù)矩陣方程化為系數(shù)矩陣為方陣的方程,然后利用Lanczos方法,求解最小二乘解.由于在求解過程中應(yīng)用到QR分解,因此稱為LSQR(Least Square QR-factorization)方法.本文利用LSQR方法求解出了外部數(shù)值保角變換模擬電荷法中的約束方程組,得到電荷量和逆變換半徑,從而構(gòu)造出近似逆保角變換函數(shù),最后利用數(shù)值實(shí)驗(yàn)驗(yàn)證了所提計(jì)算法的有效性.
本節(jié)主要闡述了基于模擬電荷法的外部區(qū)域數(shù)值保角變換計(jì)算法[12,15,17,18].如圖1所示,C對(duì)于z平面上的任意Jordan曲線,圍繞C的外部區(qū)域?yàn)镈,通過數(shù)值保角變換將區(qū)域D映射成w平面上的單位圓的外部|w|>1,=D+C.
圖1:基于模擬電荷法的外部區(qū)域數(shù)值保角變換(+代表模擬電荷點(diǎn),?代表約束點(diǎn))
保角變換函數(shù)w=f(z),f(z)滿足正規(guī)化條件f(∞)=∞,f(∞)>0時(shí),可以表示如下
g(z)是Dirichlet型場(chǎng)勢(shì)問題
的解,其中h(z)是g(z)的共軛調(diào)和函數(shù),且h(∞)=0.在下面的敘述中,F,G,H,Γ表示f,g,h,γ的近似值.由模擬電荷法,g(z)可以用C圍繞的區(qū)域內(nèi)部里配置的N 個(gè)電荷點(diǎn)ξj作為極的對(duì)數(shù)勢(shì)場(chǎng)的1次結(jié)合
高度近似g(z).此時(shí)h(z)的高度近似函數(shù)為
未知電荷Qj通過滿足下面邊界條件進(jìn)行求解,即
另外,由條件g(∞)=0,h(∞)=0,可得
因此,通過(2.5)式和(2.6)式能推導(dǎo)出以Qj(j=1,2,···,N)和logΓ作為未知數(shù)的N+1維線性方程組的構(gòu)成如下
式中,
通過式(2.3),(2.4)和(2.7)得到近似保角變換函數(shù)
最后利用zi,Qj,Γ,ζj計(jì)算雙連通保角變換.
根據(jù)上節(jié)內(nèi)容,本節(jié)提出基于模擬電荷法的外部區(qū)域數(shù)值逆保角變換計(jì)算法.如圖2在w平面上,單位圓圍成的外部區(qū)域|w|>1,通過數(shù)值保角逆變換將數(shù)值保角正變換映射成的單位圓的邊界及外部區(qū)域變換成z平面上封閉的Jordan曲線C及所圍成的外部區(qū)域D[16].
圖2:基于模擬電荷法的外部區(qū)域數(shù)值保角逆變換(+代表模擬電荷點(diǎn),?代表約束點(diǎn))
在不失一般性的情況下,假定映射函數(shù)z=f?(w)滿足正規(guī)化條件f?(∞)= ∞,f?(∞)>0 時(shí)是正則的,即
式中,γ?是變換半徑,g?(w)是Dirichlet型勢(shì)場(chǎng)問題
的解,其中 h?(w)是 g?(w)的共軛調(diào)和函數(shù),且 h?(∞)=0.分別用 F?,G?,H?,Γ?表示 f?,g?,h?,γ?的近似值.
根據(jù)模擬電荷法,g?(w)可以用單位圓的內(nèi)部區(qū)域配置的N?個(gè)電荷點(diǎn)ζ?j作為極的對(duì)數(shù)勢(shì)場(chǎng)的1次結(jié)合
高度近似,這里h?(w)可以用
高度近似.
同時(shí),由條件 g?(∞)=0,h?(∞)=0 可得
其中zi是外部的數(shù)值正保角變換的約束點(diǎn),wi是經(jīng)過zi數(shù)值正保角變換得到的映射結(jié)果,通過wi來確定.因此,由(3.5)和(3.6)式可得以(1≤j≤N?)和logΓ?作為未知數(shù)的(N?+1)維線性方程組如下
其中
通過(3.3),(3.4)和(3.7)式可得近似保角逆變換函數(shù)
最后利用zi,wi,,Γ?,計(jì)算外部區(qū)域數(shù)值保角逆變換.
將約束方程組(3.7)式寫成標(biāo)準(zhǔn)線性方程組的
(3)在今后的研究中可以繼續(xù)聯(lián)合實(shí)地監(jiān)測(cè)數(shù)據(jù),除植被因素外,將景觀要素和土壤要素以及周邊居民滿意度等要素,在生態(tài)重建效果評(píng)價(jià)中的重要性考慮進(jìn)去。另外下一步工作中可以進(jìn)一步結(jié)合多種評(píng)價(jià)方法,例如和層次分析法、灰色關(guān)聯(lián)度法、聚類分析法、模糊綜合評(píng)價(jià)法等做對(duì)比,對(duì)研究區(qū)的生態(tài)重建效果進(jìn)行全面評(píng)價(jià)比較和分析。
形式,其中
約束方程的系數(shù)矩陣A是非對(duì)稱的且病態(tài)的,LSQR方法[20?23]是求解系數(shù)矩陣為病態(tài)的大型稀疏矩陣線性方程組的有效算法之一.利用Lanczos雙對(duì)角化方法來求解方程的最小二乘解 minkAx?bk2. 假定 Uk=[u1,···,uk]和 Vk=[v1,···,vk]是正交陣且 Lk為如下的(k+1)×k的下雙對(duì)角陣
用下列迭代方法可以實(shí)現(xiàn)A矩陣的雙對(duì)角分解
其中αi≥0,βi≥0.使上式(3.21)可寫成
可以確定
在滿足給定精度時(shí)停止迭代.我們希望krkk2盡量小,且Uk+1理論上是正交陣,取yk使ktk+1k2最小,解最小二乘問題minkβ1e1?Lkykk2.得到LSQR算法如下:
Algorithm 1 LSQR Algorithm Input:A,b,x0,ε.Initialize β1u1=b,α1v1=ATu1,h1=v1,e?1=β1,eρ1=α1.for i=1,2,3,···while stopping criterion is not satisfied do βi+1ui+1=Avi?αiui;αi+1ui+1=ATi?βi+1Vi;ρi=eρ2i+β2i+1,ci=eρi/ρi,si=βi+1/ρi;θi+1=siαi+1,eρi+1=?ciαi+1,?i=cie?i,e?i+1=sie?i;xi=xi?1+(e?i/eρi)hi;hi+1=vi+1? (θi+1/ρi)hi;if minkAxi?bk<ε;end if end while end for Output xi.q
上述算法中,ε是給定精度.
這里給出基于LSQR法的外部區(qū)域數(shù)值保角逆變換計(jì)算法的具體步驟如下.
步驟1通過外部區(qū)域數(shù)值保角正變換(2.9)得到映射點(diǎn)F(zi),將F(zi)的位置作為數(shù)值保角逆變換的約束點(diǎn)wi的位置.
步驟2根據(jù)約束點(diǎn)wi的位置配置外部區(qū)域數(shù)值保角逆變換模擬電荷點(diǎn)的位置.
步驟3 通過LSQR 方法求解約束方程組(3.7)得到模擬電荷,,···,和逆變換半徑 logΓ?.
步驟4對(duì)單位圓的邊界及外部區(qū)域的每一個(gè)點(diǎn)通過(3.3)和(3.4)式計(jì)算得到G?(w)和H?(w)后,構(gòu)造近似保角逆變換函數(shù)(3.9),然后計(jì)算對(duì)應(yīng)的變換點(diǎn).
針對(duì)橙形為邊界的外部區(qū)域,在MATLAB 13b環(huán)境下,檢驗(yàn)雙方向的外部區(qū)域數(shù)值保角變換計(jì)算方法的有效性.基于模擬電荷法的單連通區(qū)域的外部區(qū)域數(shù)值正保角變換的誤差由Ez=max(||f(z)|?1|)確定,外部數(shù)值保角逆變換的誤差由Ew=max(|f?(w)?z|)確定.
例1橙形邊界及外部區(qū)域的雙方向數(shù)值保角變換.邊界
約束點(diǎn)的位置由
確定,其中i=0,1,···,N?1.約束點(diǎn)分布在邊界上,邊界由粗實(shí)線表示.
保角正變換模擬電荷點(diǎn)的位置由下面公式給出
保角逆變換模擬電荷點(diǎn)的位置由下面公式給出
其中a=21/16,rz=rw=3,r>0是確定模擬電荷點(diǎn)位置的參數(shù),模擬電荷點(diǎn)的分布在區(qū)域D的外部,如圖3所示為數(shù)值保角正變換的模擬電荷點(diǎn)分布.如圖4所示為數(shù)值保角逆變換的模擬電荷點(diǎn)分布.約束點(diǎn)和模擬電荷點(diǎn)一一對(duì)應(yīng),數(shù)量都為N.
圖3:橙形保角變換邊界及電荷點(diǎn)位置
圖4:橙形保角逆變換邊界及電荷點(diǎn)位置
圖5:橙形保角逆變換誤差曲線
圖6:橙形保角逆變換誤差曲線
分別用Amano和lsqr分別表示天野法和基于LSQR方法的單連通外部區(qū)域數(shù)值保角逆變換計(jì)算法.圖5給出的是當(dāng)a=21/16,rz=rw=3兩種數(shù)值保角逆變換方法的誤差曲線,圖5說明誤差隨著電荷量的增大而減小,同時(shí)可看出lsqr的誤差值一直小于Amano的誤差值,在N=61時(shí),Amano誤差為2.4243×10?5而lsqr誤差為1.3025×10?6,說明了本文采用的算法可以得到更高的誤差精度,驗(yàn)證了算法的有效性.
圖6是當(dāng)a=21/14,rz=rw=2時(shí)兩種方法的誤差曲線,由圖6可看出電荷點(diǎn)越多誤差值越小,且各個(gè)數(shù)量的電荷量上,lsqr的誤差值均比Amano的誤差值小.電荷點(diǎn)數(shù)為180時(shí),Amano誤差為1.0366×10?7,lsqr誤差為1.6639×10?8,因此數(shù)值實(shí)驗(yàn)再次驗(yàn)證了外部數(shù)值保角逆變換計(jì)算法的有效性.
圖7:橙形邊界和外部區(qū)域及等高線
圖8:圖7的保角變換
圖9:圖10的保角逆變換
圖10:圖8的邊界和其外部區(qū)域及其等高線
圖7–10中的粗實(shí)線表示邊界,細(xì)實(shí)線表示等高線.圖7是橙形的邊界及外部區(qū)域等高線,圖8是圖7通過近似保角變換函數(shù)F(z)映射后得到的結(jié)果.由圖8可知,近似保角變換函數(shù)F(z)將橙形的邊界映射成了單位圓.圖10表示的是單位圓的邊界和外部區(qū)域|w|>1及其等高線,以及模擬電荷點(diǎn)的配置位置.圖9是圖10通過近似保角逆變換函數(shù)F?(w)映射后得到的結(jié)果.由圖9和圖10可知,近似逆保角變換函數(shù)F?(w)將橙單位圓邊界映射成了橙形邊界,實(shí)現(xiàn)了數(shù)值保角逆變換.
本文利用LSQR方法求解出了基于模擬電荷法的外部數(shù)值保角逆變換中的約束方程組,提出了LSQR方法的外部數(shù)值保角逆變換計(jì)算法.并通過數(shù)值實(shí)驗(yàn)驗(yàn)證了所提計(jì)算法的有效性并用等高線模擬了外部數(shù)值保角逆變換的計(jì)算結(jié)果.在今后的研究中,本方法同樣可以應(yīng)用于多連通區(qū)域的數(shù)值保角逆變換問題.