保角变换是复变函数的一个基本问题, 广泛应用于物理学与工学^[1-3].保角变换的主要求解方法:解析法和数值计算法.解析法指出了变换函数的存在, 但是只能给出一些特殊区域的变换函数表达式.基于实际问题的复杂性, 必须通过数值计算方法求变换函数.目前的保角变换分为两类:一类是单连通区域的保角变换, 另一类是多连通区域的保角变换.

1969年, 德国的Steinbigler提出了在计算空气中回转对称电极周围电场时可以用设置在电极内部的若干个虚设的模拟电荷来计算电极表面上电荷分布的电场的计算方法, 形成了模拟电荷法(Charge Simulation Method)的基本思想^[4].自20世纪80年代起, 天野要等对模拟电荷法、数值保角变换、数值保角变换的模拟电荷法以及数值保角变换中模拟电荷点的计算方法做了许多研究工作, 提出了基于模拟电荷的保角变换计算法(天野法)^[5-8].在复杂的数值保角变换的计算中, 相比传统的数值计算法, 基于模拟电荷法的数值保角变换计算法具有计算精度高、误差评价容易、计算时间短等优势.

广义极小残量法(Generalized Minimum Residual Methods)^[9-11]是在Krylov子空间进行迭代的计算方法, 但是GMRES法的储存量和正交化工作量会随着迭代次数的增加而增加, 因此, 一般需要利用限制Krylov子空间最大维数的算法--GMRES($m$)法求解计算问题.本文利用GMRES($m$)求解出了双连通区域模拟电荷法中的约束方程, 得到了模拟电荷和变换半径, 从而得到近似保角变换函数, 最后利用数值实验验证了算法的有效性和结果的正确性.

本节主要阐述基于模拟电荷法的双连通区域保角变换的数值计算法^[12-13].对于$z$平面上的一个有限双连通区域$D$, 它由两条封闭的Jordan曲线$C_{1}$和$C_{2}$所围成.通过保角变换把它变成了$w$平面的一个圆环$\mu<\mid\omega\mid<1$, $C_{1}$和$C_{2}$分别是外部边界和内部边界.

在不失一般性的情况下, 假定$f(0)=0$, $f(z)$满足正规化条件$f(\infty)=\infty$, $f'(\infty)>0$时, 可以表示成

$ \begin{equation} \label{eq:1} f(z)=z e^{g(z)+ih(z)}, \quad z\in \bar{D} =D+C_{1}+C_{2}, \end{equation} $

(2.1)

$g(z)$是Dirichlet型势场问题

$ \left\{ \begin{aligned} &\nabla^{2}g(z) =0, &z\in D, \\ &g(z) =-\log|z|,&z\in C_{1}, \\ &g(z) =\log\mu - \log|z|,&z\in C_{2} \end{aligned} \right. $

(2.2)

的解, $h(z)$是$g(z)$的共轭调和函数^[14]由模拟电荷法可知$g(z)$可以用$C_{1}$和$C_{2}$所围成的区域外部配置$N$个电荷点作为极的对数势场的一次结合

$ G(z) = - \sum\limits_{i = 1}^N {{Q_i}} \log |z - {\zeta _i}|,z \in \bar D $

(2.3)

高度近似.因此$g(z)$在区域$D$内的共轭调和函数$h(z)$可以用

$ H(z)=-\sum\limits_{i=1}^N Q_i \arg(z-\zeta_{i}), z\in\bar{D} $

(2.4)

近似.另外, $\mu$由$M$近似. $\zeta_{i}(i=1, 2, \cdots, N)$为电荷点, 分布在给定区域的外部, 即$N/2$个点分布在边界$C_{1}$的外部, 另外$N/2$个点分布在边界$C_{2}$的内部.根据(2.2) 式, 双连通区域的数值保角变换是单连通区域的内部保角变换与外部保角变换的组合, 因此边界条件可以看成内部保角变换边界条件与外部保角变换边界条件的组合.未知电荷$Q_{i}$可以用边界上选择的约束点$z_{j}$在满足Dirichlet边界条件时进行求解, 即

$ \sum\limits_{i=1}^N Q_i \log|z_{j}-\zeta_{i}|=\log|z_{j}|, \quad z_{j}\in C_{1}, $

(2.5)

$ \sum\limits_{i=1}^N Q_i \log|z_{j}-\zeta_{i}|+\log M =\log|z_{j}|, z_{j}\in C_{2}. $

(2.6)

又由于$g(\infty)=0$, $h(\infty)=0$, 根据(2.3) 和(2.4) 式可推导出

$ \sum\limits_{i=N/2+1}^{N}Q_i =0. $

(2.7)

上述(2.5) 式中$z_{j}~(j=1, 2, \cdots, N/2)$取在边界$C_{1}$上.类似的, (2.6) 式中$z_{j}~(j=N/2+1, N/2+2, \cdots, N)$取在边界$C_{2}$上.

联立(2.5)-(2.7) 式, 化为$N+1$维约束方程

$ \begin{equation} \left( \begin{array}{cccccc} a_{11}& \cdots& a_{1, N/2+1} & \cdots& a_{1, N} & 0 \\ \vdots& \ddots& \vdots & \ddots& \vdots & \vdots \\ a_{N/2+1, 1}& \cdots& a_{N/2+1, N/2+1} & \cdots& a_{N/2+1, N} & 1 \\ \vdots& \vdots& \vdots & \ddots& \vdots & \vdots \\ a_{N1}& \cdots& a_{N, N/2+1} & \cdots& a_{NN} & 1 \\ 0 & \cdots& 1 & \cdots& 1 & 0 \end{array} \right) \left( \begin{array}{c} Q_{1} \\ \vdots \\ Q_{N/2+1} \\ \vdots \\ Q_{N} \\ \log M \end{array} \right) = \left( \begin{array}{c} \log|z_{1}| \\ \vdots \\ \log|z_{N/2+1}| \\ \vdots \\ \log|z_N| \\ 0 \end{array} \right), \end{equation} $

(2.8)

其中$a_{ij}=\log|z_{j} - \zeta_{i}|$.最后, 利用$z_{i}, \zeta_{i}, Q_{i}, M$计算双连通保角变换.

将约束方程(2.8) 式写成标准线性方程组$Ax=b$的形式, 其中

$ A\in R^{(N+1)\times (N+1)}, ~~ x\in R^{N+1}, ~~ b\in R^{N+1}. $

约束方程的系数矩阵$A$非对称的, 广义极小残量法(GMRES)是当前求解大型非对称线性方程非常有效的算法之一.它是在Krylov子空间

$ K_{m}(A, r_{0})= \mbox{span}\{r_{0}, Ar_{0}, Ar^{2}_{0}, \cdots, Ar^{m-1}_{0}\} $

的基础上发展起来的算法, 通过寻找近似解$x_{m}\in x_{0}+K_{m}(A, r_{0})$, 来逼近约束方程$(2.8)$的准确解.但是GMRES法在求解约束方程时, 随着电荷点数的增加, 计算量和储存量都会随着迭代次数的增加而大大增加.为了避免这种情况, 每$m$步必须重新启动, 这就是GMRES($m$)法^[15]. GMRES($m$)算法是GMRES算法的进一步改进, 根据文献[16], 可以得到求解模拟电荷的GMRES($m$)法, 其算法步骤如下:

上述算法中, $x_{0}$是给定的初始近似解, 取为零向量, $r_{0}=b-Ax_{0}$为初始向量, $\varepsilon$是给定精度.

在MATLAB 7.0环境下, 以椭圆为边界的双连通区域$C_{1}:x^{2}/a_{1}^{2} + y^{2}/b_{1}^{2}=1$, $C_{2}:x^{2}/a_{2}^{2} + y^{2}/b_{2}^{2}=1$为例, $C_{1}$和$C_{2}$为边界.利用模拟电荷法对双连通区域的保角变换进行数值实验, 保角变换误差由

$ {\rm{error}} = \max \left( {\mathop {\max }\limits_{_{{c_1}}} \left| {|f(z)| - 1} \right|,\mathop {\max }\limits_{_{{c_2}}} ||f(z)| - M|} \right) $

确定, 检验GMRES($m$)法的有效性, 具体步骤如下:

步骤1  设定电荷点$\zeta_{i}(i=1, 2, \cdots, N)$, 约束点$z_{i}(i=1, 2, \cdots, N)$以及各个参数;

步骤2  根据边界条件确定关于模拟电荷$Q_{i}(i=1, 2, \cdots, N)$和变换半径$M$的约束方程；

步骤3  通过GMRES($m$)法求解约束方程, 得到模拟电荷$Q_{i}(i=1, 2, \cdots, N)$和变换半径$M$, 从而构造近似保角变换函数$f(z)$.

例1  $a_{1}=7$, $b_{1}=5$, $a_{2}=5$, $b_{2}=1$时.分别用Method 1和Method 2表示天野法和基于GMRES($m$)法的双连通区域保角变换的计算法. 图 2给出的是这两种方法数值保角变换误差.由图 2可知电荷点数增加时, 误差会随着电荷点数的增加而减小.当电荷点数多于90时, Method 2的保角变换误差要小于Method 1.电荷点数$N=200$时, Method 2的保角变换误差为$4.1729\times 10^{-6}$, 而Method 1的保角变换误差为$8.7557\times 10^{-6}$, 说明本文所采用的算法可以得到更高精度的模拟电荷$Q_{i}~(i=1, 2, \cdots, N)$和变换半径$M$. 图 3给出的模拟电荷点数$N$为200时, 模拟电荷点的分布情况.

为进一步验证本文算法的有效性, 取以

$ C_{1}:x^{2}/7^{2} + y^{2}/5^{2}=1, ~~ C_{2}:x^{2}/5^{2} + y^{2}=1 $

为边界的内部区域的等高线. 图 4中, 粗实线表示的边界, 细线表示的等高线.电荷点数取$200$, 利用GMRES($m$)法求解约束方程, 从而构造变换函数. 图 5为区域及其等高线的映射结果, 可以看出, 对于$C_{1}$和$C_{2}$的所围成的区域内的任意部分, 经过保角变换, 对应的仍然是变换后所围成区域的内部.且边界经过保角变换对应的同心圆也是变换后图像的边界, 从而验证了基于GMRES($m$)法的双连通区域保角变换的计算法的有效性.

例2  $a_{1}=8$, $b_{1}=4$, $a_{2}=4$, $b_{2}=2$时.类似于图 2和图 3, 图 6和图 7分别给出了两种方法的数值保角变换误差和模拟电荷点的分布.从图 6中可以看出, 电荷点数越多, 两种方法的保角变换误差越小.但是, 当电荷点数大于180时, Method 2的保角变换误差明显要小于Method 1.电荷点数为300, Method 2的保角变换误差已经达到$3.7814\times10^{-13}$. 图 8中细线表示以

$ C_{1}:x^{2}/8^{2} + y^{2}/4^{2}=1, ~~C_{2}:x^{2}/4^{2} + y^{2}/2^{2}=1 $

为边界的内部区域的等高线, 粗线表示边界.从图 9可以看出, Method 2可以得到了很好的保角变换结果.

本文利用GMRES($m$)法求解出了双连通区域保角变换模拟电荷法中的约束方程, 进而提出了基于GMRES($m$)法的双连通区域数值保角变换的计算法.并通过数值实验验证了新算法的有效性, 且用等高线模拟了保角变换的计算结果.本算法同样可以运用解决其它的多连通区域数值保角变换问题.