基于LSQR法的外部数值保角逆变换计算法

吕毅斌,王 坚,王樱子,吴 爽

(1.昆明理工大学理学院,云南昆明 650500)

(2.昆明理工大学计算中心,云南昆明 650500)

1 引言

保角变换是复变函数中非常重要的理论之一,广泛应用于物理学和工学领域.特别是在电磁理论、膜和板的振动、弹性理论、热传输、流体力学等方面有很多应用[1−5].通常将保角变换的求解方法分为解析法和数值法.解析法只能在一些特殊区域给出变换函数表达式,对于复杂区域问题没有解决方法.因此,对于很多实际中的复杂问题必须采用数值法求解保角变换函数.数值保角变换计算法的主要有:积分方程式法[6]、正交多项式法[7−8]和有限差分法[9−10]等.德国人Steinbigler[11]首次提出用若干个虚设电荷来模拟电极表面上电荷分布的电场计算法,形成了模拟电荷法的基本思想;日本的天野要等数学学者从20世纪80年代开始对模拟电荷法和数值保角变换作了大量研究工作,并提出了基于模拟电荷法的数值保角变换计算法(天野法)[12−16].天野法适用于单连通区域及多连通区域的数值保角变换问题[15,17,18,19,20].

单连通区域的数值保角变换,分为内部数值保角变换和外部数值保角变换.本文通过对基于模拟电荷法的双方向的内部数值保角变换计算法[13,15,16]的研究,提出了外部数值保角逆变换计算法.该方法的原理是基于模拟电荷法来求Laplace方程的Dirichlet问题的解,并通过预先建立的边界对应关系构造从标准区域到问题区域的近似保角逆变换函数,误差用正则函数的最大值原理进行评价.计算数值保角逆变换最主要的就是确定模拟电荷的位置和数量,以及约束方程组的求解.

LSQR方法[21−24]是Paige和Saunders提出的一种适用于求解系数矩阵为大型、稀疏矩阵线性方程组的方法.LSQR方法求解的思路是把任意系数矩阵方程化为系数矩阵为方阵的方程,然后利用Lanczos方法,求解最小二乘解.由于在求解过程中应用到QR分解,因此称为LSQR(Least Square QR-factorization)方法.本文利用LSQR方法求解出了外部数值保角变换模拟电荷法中的约束方程组,得到电荷量和逆变换半径,从而构造出近似逆保角变换函数,最后利用数值实验验证了所提计算法的有效性.

2 基于模拟电荷法的外部数值保角变换计算法

本节主要阐述了基于模拟电荷法的外部区域数值保角变换计算法[12,15,17,18].如图1所示,C对于z平面上的任意Jordan曲线,围绕C的外部区域为D,通过数值保角变换将区域D映射成w平面上的单位圆的外部|w|>1,=D+C.

图1:基于模拟电荷法的外部区域数值保角变换（+代表模拟电荷点,•代表约束点）

保角变换函数w=f(z),f(z)满足正规化条件f(∞)=∞,f(∞)>0时,可以表示如下

g(z)是Dirichlet型场势问题

的解,其中h(z)是g(z)的共轭调和函数,且h(∞)=0.在下面的叙述中,F,G,H,Γ表示f,g,h,γ的近似值.由模拟电荷法,g(z)可以用C围绕的区域内部里配置的N 个电荷点ξj作为极的对数势场的1次结合

高度近似g(z).此时h(z)的高度近似函数为

未知电荷Qj通过满足下面边界条件进行求解,即

另外,由条件g(∞)=0,h(∞)=0,可得

因此,通过(2.5)式和(2.6)式能推导出以Qj(j=1,2,···,N)和logΓ作为未知数的N+1维线性方程组的构成如下

式中,

通过式(2.3),(2.4)和(2.7)得到近似保角变换函数

最后利用zi,Qj,Γ,ζj计算双连通保角变换.

3 外部区域数值保角逆变换计算法

根据上节内容,本节提出基于模拟电荷法的外部区域数值逆保角变换计算法.如图2在w平面上,单位圆围成的外部区域|w|>1,通过数值保角逆变换将数值保角正变换映射成的单位圆的边界及外部区域变换成z平面上封闭的Jordan曲线C及所围成的外部区域D[16].

图2:基于模拟电荷法的外部区域数值保角逆变换（+代表模拟电荷点,•代表约束点）

在不失一般性的情况下,假定映射函数z=f∗(w)满足正规化条件f∗(∞)= ∞,f∗(∞)>0 时是正则的,即

式中,γ∗是变换半径,g∗(w)是Dirichlet型势场问题

的解,其中 h∗(w)是 g∗(w)的共轭调和函数,且 h∗(∞)=0.分别用 F∗,G∗,H∗,Γ∗表示 f∗,g∗,h∗,γ∗的近似值.

根据模拟电荷法,g∗(w)可以用单位圆的内部区域配置的N∗个电荷点ζ∗j作为极的对数势场的1次结合

高度近似,这里h∗(w)可以用

高度近似.

同时,由条件 g∗(∞)=0,h∗(∞)=0 可得

其中zi是外部的数值正保角变换的约束点,wi是经过zi数值正保角变换得到的映射结果,通过wi来确定.因此,由(3.5)和(3.6)式可得以(1≤j≤N∗)和logΓ∗作为未知数的(N∗+1)维线性方程组如下

其中

通过(3.3),(3.4)和(3.7)式可得近似保角逆变换函数

最后利用zi,wi,,Γ∗,计算外部区域数值保角逆变换.

4 基于LSQR方法的数值保角逆变换模拟电荷求解

将约束方程组(3.7)式写成标准线性方程组的

(3)在今后的研究中可以继续联合实地监测数据，除植被因素外，将景观要素和土壤要素以及周边居民满意度等要素，在生态重建效果评价中的重要性考虑进去。另外下一步工作中可以进一步结合多种评价方法，例如和层次分析法、灰色关联度法、聚类分析法、模糊综合评价法等做对比，对研究区的生态重建效果进行全面评价比较和分析。

形式,其中

约束方程的系数矩阵A是非对称的且病态的,LSQR方法[20−23]是求解系数矩阵为病态的大型稀疏矩阵线性方程组的有效算法之一.利用Lanczos双对角化方法来求解方程的最小二乘解 minkAx−bk2. 假定 Uk=[u1,···,uk]和 Vk=[v1,···,vk]是正交阵且 Lk为如下的(k+1)×k的下双对角阵

用下列迭代方法可以实现A矩阵的双对角分解

其中αi≥0,βi≥0.使上式(3.21)可写成

可以确定

在满足给定精度时停止迭代.我们希望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法的外部区域数值保角逆变换计算法的具体步骤如下.

步骤1通过外部区域数值保角正变换(2.9)得到映射点F(zi),将F(zi)的位置作为数值保角逆变换的约束点wi的位置.

步骤2根据约束点wi的位置配置外部区域数值保角逆变换模拟电荷点的位置.

步骤3 通过LSQR 方法求解约束方程组(3.7)得到模拟电荷,,···,和逆变换半径 logΓ∗.

步骤4对单位圆的边界及外部区域的每一个点通过(3.3)和(3.4)式计算得到G∗(w)和H∗(w)后,构造近似保角逆变换函数(3.9),然后计算对应的变换点.

5 数值实验

针对橙形为边界的外部区域,在MATLAB 13b环境下,检验双方向的外部区域数值保角变换计算方法的有效性.基于模拟电荷法的单连通区域的外部区域数值正保角变换的误差由Ez=max(||f(z)|−1|)确定,外部数值保角逆变换的误差由Ew=max(|f∗(w)−z|)确定.

例1橙形边界及外部区域的双方向数值保角变换.边界

约束点的位置由

确定,其中i=0,1,···,N−1.约束点分布在边界上,边界由粗实线表示.

保角正变换模拟电荷点的位置由下面公式给出

保角逆变换模拟电荷点的位置由下面公式给出

其中a=21/16,rz=rw=3,r>0是确定模拟电荷点位置的参数,模拟电荷点的分布在区域D的外部,如图3所示为数值保角正变换的模拟电荷点分布.如图4所示为数值保角逆变换的模拟电荷点分布.约束点和模拟电荷点一一对应,数量都为N.

图3:橙形保角变换边界及电荷点位置

图4:橙形保角逆变换边界及电荷点位置

图5:橙形保角逆变换误差曲线

图6:橙形保角逆变换误差曲线

分别用Amano和lsqr分别表示天野法和基于LSQR方法的单连通外部区域数值保角逆变换计算法.图5给出的是当a=21/16,rz=rw=3两种数值保角逆变换方法的误差曲线,图5说明误差随着电荷量的增大而减小,同时可看出lsqr的误差值一直小于Amano的误差值,在N=61时,Amano误差为2.4243×10−5而lsqr误差为1.3025×10−6,说明了本文采用的算法可以得到更高的误差精度,验证了算法的有效性.

图6是当a=21/14,rz=rw=2时两种方法的误差曲线,由图6可看出电荷点越多误差值越小,且各个数量的电荷量上,lsqr的误差值均比Amano的误差值小.电荷点数为180时,Amano误差为1.0366×10−7,lsqr误差为1.6639×10−8,因此数值实验再次验证了外部数值保角逆变换计算法的有效性.

图7:橙形边界和外部区域及等高线

图8:图7的保角变换

图9:图10的保角逆变换

图10:图8的边界和其外部区域及其等高线

图7–10中的粗实线表示边界,细实线表示等高线.图7是橙形的边界及外部区域等高线,图8是图7通过近似保角变换函数F(z)映射后得到的结果.由图8可知,近似保角变换函数F(z)将橙形的边界映射成了单位圆.图10表示的是单位圆的边界和外部区域|w|>1及其等高线,以及模拟电荷点的配置位置.图9是图10通过近似保角逆变换函数F∗(w)映射后得到的结果.由图9和图10可知,近似逆保角变换函数F∗(w)将橙单位圆边界映射成了橙形边界,实现了数值保角逆变换.

6 结束语

本文利用LSQR方法求解出了基于模拟电荷法的外部数值保角逆变换中的约束方程组,提出了LSQR方法的外部数值保角逆变换计算法.并通过数值实验验证了所提计算法的有效性并用等高线模拟了外部数值保角逆变换的计算结果.在今后的研究中,本方法同样可以应用于多连通区域的数值保角逆变换问题.