保角变换是复变函数中非常重要的理论之一, 广泛应用于物理学和工学领域.特别是在电磁理论、膜和板的振动、弹性理论、热传输、流体力学等方面有很多应用^[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方法求解出了外部数值保角变换模拟电荷法中的约束方程组, 得到电荷量和逆变换半径, 从而构造出近似逆保角变换函数, 最后利用数值实验验证了所提计算法的有效性.

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

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

$ \begin{equation} f(z) = \frac{z}{\gamma} e^{g(z)+ih(z)}, \quad z\in \bar{D}, \end{equation} $

(2.1)

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

$ \begin{equation} \left\{ \begin{aligned} & \nabla^{2}g(z) = 0, \ & z\in D, \\ & g(z) = -\log|\gamma|-\log|z|, & z\in C, \\ & g(\infty) = 0 \end{aligned} \right. \end{equation} $

(2.2)

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

$ \begin{eqnarray} G(z) = -\sum\limits_{j = 1}^N Q_j \log|z-\zeta_{j}|, \quad z\in \bar{D}, \end{eqnarray} $

(2.3)

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

$ \begin{eqnarray} H(z) = -\sum\limits_{j = 1}^N Q_j {\rm arg}|z-\zeta_{j}|, \quad z\in\bar{D}, \end{eqnarray} $

(2.4)

未知电荷$ Q_{j} $通过满足下面边界条件进行求解, 即

$ \begin{eqnarray} &&\sum\limits_{j = 1}^N Q_j \log|z_{i}-\zeta_{j}| = \log \Gamma-\log|z_{i}| \quad (i = 1, 2, \cdots, N). \end{eqnarray} $

(2.5)

另外, 由条件$ g(\infty) = 0 $, $ h(\infty) = 0 $, 可得

$ \begin{equation} \sum\limits_{j = 1}^{N}Q_j = 0. \end{equation} $

(2.6)

因此, 通过(2.5)式和(2.6)式能推导出以$ Q_{j}\; (j = 1, 2, \cdots, N) $和$ \log \Gamma $作为未知数的$ N+1 $维线性方程组的构成如下

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

(2.7)

式中,

$ \begin{equation} a_{ij} = \log|z_{j} - \zeta_{i}|, \end{equation} $

(2.8)

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

$ \begin{eqnarray} F(z) = \frac{z}{\Gamma}\exp(G(z)+iH(z)), & z\in \bar D. \end{eqnarray} $

(2.9)

最后利用$ z_{i} $, $ Q_{j} $, $ \Gamma $, $ \zeta_{j} $计算双连通保角变换.

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

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

$ \begin{eqnarray} f^*(w) = \frac{w}{\gamma^*}\exp(g^*(w)+ih^*(w)), & |w|>1, \end{eqnarray} $

(3.1)

式中, $ \gamma^* $是变换半径, $ g^*(w) $是Dirichlet型势场问题

$ \begin{equation} \left\{ \begin{aligned} &\nabla^2g^*(w) = 0, \quad w\in D, \\ & g^*(w) = \log|z|-\log|w|-\log\gamma^*, &|w|>1, \\ & g^*(\infty) = 0 \end{aligned} \right. \end{equation} $

(3.2)

的解, 其中$ h^*(w) $是$ g^*(w) $的共轭调和函数, 且$ h^*(\infty) = 0 $.分别用$ F^* $, $ G^* $, $ H^* $, $ \Gamma^* $表示$ f^* $, $ g^* $, $ h^* $, $ \gamma^* $的近似值.

根据模拟电荷法, $ g^*(w) $可以用单位圆的内部区域配置的$ N^{*} $个电荷点$ \zeta_j^* $作为极的对数势场的1次结合

$ \begin{eqnarray} G^{*}(w) = -\sum\limits_{j = 1}^{N^{*}}Q^{*}_{j} \log|w - \zeta^{*}_{j}| &|w|>1, \end{eqnarray} $

(3.3)

高度近似, 这里$ h^*(w) $可以用

$ \begin{eqnarray} H^{*}(w) = -\sum\limits_{j = 1}^{N^{*}}Q_{j}^{*} \arg(1-\zeta_{j}^{*}/w), & |w|>1 \end{eqnarray} $

(3.4)

高度近似.

未知电荷$ Q_j^* $可以通过满足下面的边界条件进行求解,

$ \begin{eqnarray} \sum\limits_{j = 1}^{N^{*}}Q_{j}^{*} \log|w_{i}-\zeta_{j}^{*}| = \log|w_i|-\log|z_i|-\log\Gamma^{*}, &j = 1, 2, \cdots, N^{*}, \end{eqnarray} $

(3.5)

同时, 由条件$ g^*(\infty) = 0 $, $ h^*(\infty) = 0 $可得

$ \begin{eqnarray} \sum\limits_{j = 1}^{N^{*}}Q_{j}^{*} = 0, \end{eqnarray} $

(3.6)

其中$ z_i $是外部的数值正保角变换的约束点, $ w_i $是经过$ z_i $数值正保角变换得到的映射结果, $ \zeta^*_j $通过$ w_i $来确定.因此, 由(3.5)和(3.6)式可得以$ Q^*_j\; (1 \leq j\leq N^{*}) $和$ \log\Gamma^* $作为未知数的$ (N^{*}+1) $维线性方程组如下

$ \begin{equation} \left( \begin{array}{cccccc} a_{11} & a_{12} & \cdots & a_{1N^{*}} & 1 \\ a_{21} & a_{22} & \ddots & a_{2N^{*}} & 1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{N^{*}1} & a_{N^{*}2} & \cdots & a_{N^{*}N^{*}} & 1 \\ 1 & 1 & \cdots & 1 & 0 \end{array} \right) \left( \begin{array}{c} Q_{1}^{*} \\ Q_{2}^{*} \\ \vdots \\ Q_{N^{*}}^{*} \\ \log {\Gamma^{*}} \end{array} \right) = \left( \begin{array}{c} \log|w_{1}|- \log|z_{1}| \\ \log|w_{2}|- \log|z_{2}| \\ \vdots \\ \log|w_{N^{*}}| - \log|z_{N^{*}}| \\ 0 \end{array} \right), \end{equation} $

(3.7)

其中

$ \begin{eqnarray} a_{ij} = \log|w_{i}-\zeta^{*}_{j}|. \end{eqnarray} $

(3.8)

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

$ \begin{eqnarray} F^*(w) = \frac{w}{\Gamma^*}\exp(G^*(w)+iH^*(w)), & |w|>1. \end{eqnarray} $

(3.9)

最后利用$ z_{i} $, $ w_{i} $, $ Q^{*}_{j} $, $ \Gamma^{*} $, $ \zeta^{*}_{j} $计算外部区域数值保角逆变换.

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

$ \begin{eqnarray} Ax = b \end{eqnarray} $

(4.1)

形式, 其中

$ A\in R(N+1)\times(N+1), x\in R(N+1), b\in R(N+1). $

约束方程的系数矩阵$ A $是非对称的且病态的, LSQR方法^[20-23]是求解系数矩阵为病态的大型稀疏矩阵线性方程组的有效算法之一.利用Lanczos双对角化方法来求解方程的最小二乘解$ \min{\|Ax-b\|_{2}} $.假定$ U_k = [u_1, \cdot\cdot\cdot, u_k] $和$ V_k = [v_1, \cdot\cdot\cdot, v_k] $是正交阵且$ L_k $为如下的$ (k+1)\times k $的下双对角阵

$ \begin{equation} \begin{array}{c} L_k \end{array} = \left( \begin{array}{cccccc} \alpha_{1} & 0 & \cdots & \cdots & \cdots & 0 \\ \beta_{2} & \alpha_{2} & 0 & \cdots & \cdots & 0 \\ 0 & \beta_{3} & \alpha_{3} & \cdots & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \alpha_{k} \\ 0 & \cdots & \cdots & \cdots & 0 & \beta_{k+1} \end{array} \right), \end{equation} $

(4.2)

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

$ \begin{equation} \left\{ \begin{aligned} &\beta_1 u_1 = b, \quad\alpha_1 v_1 = A^Tu_1, \\ & \beta_{i+1} = Av_i-\alpha_iu_i\quad(i = 1, 2, \cdots), \\ & \alpha_{i+1}v_{i+1} = A^T\beta_{i+1}v_i, \end{aligned} \right. \end{equation} $

(4.3)

其中$ \alpha_i\geq0 $, $ \beta_i\geq0 $.使上式(3.21)可写成

$ \begin{equation} \left\{ \begin{aligned} &U_{k+1}(\beta_1e_1) = b, \\ & AV_k = U_{k+1}L_k, \\ & A^TU_{k+1} = V_kL_k^T+\alpha_{k+1}v_{k+1}e_{k+1}^T, \end{aligned} \right. \end{equation} $

(4.4)

其中$ e_{k+1}^T $, 表示$ n $阶单位矩阵的第$ k+1 $行, 再设

$ \begin{eqnarray} x_k = V_ky_k, r_k = b-Ax_k, t_{k+1} = \beta_1e_1-L_ky_k, \end{eqnarray} $

(4.5)

可以确定

$ \begin{eqnarray} \begin{split} r_k = b-A_x = U_{k+1}(\beta_1e_1)-AV_ky_k = U_{k+1}(\beta_1e_1)-U_{k+1}L_ky_k = U_{k+1}t_{k+1} \end{split} \end{eqnarray} $

(4.6)

在满足给定精度时停止迭代.我们希望$ {\|r_k\|_{2}} $尽量小, 且$ U_{k+1} $理论上是正交阵, 取$ y_k $使$ {\|t_{k+1}\|_{2}} $最小, 解最小二乘问题$ \min{\|\beta_1e_1-L_ky_k\|_{2}} $.得到LSQR算法如下:

上述算法中, $ \varepsilon $是给定精度.

这里给出基于LSQR法的外部区域数值保角逆变换计算法的具体步骤如下.

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

步骤2  根据约束点$ w_{i} $的位置配置外部区域数值保角逆变换模拟电荷点$ \zeta^{*}_{j} $的位置.

步骤3  通过LSQR方法求解约束方程组$ (3.7) $得到模拟电荷$ Q^{*}_{1} $, $ Q^{*}_{2} $, $ \cdots, $ $ Q^{*}_{N^{*}} $和逆变换半径$ \log\Gamma^{*} $.

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

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

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

$ \begin{eqnarray} \{(x+1)^{2}+y^{2}\}\{(x-1)^{2}+y^{2}\} = a^{4}. \end{eqnarray} $

(5.1)

约束点的位置由

$ \begin{eqnarray} \left\{ \begin{aligned} &x_{i} = r\cos(\theta), \\ & y_{i} = r\sin(\theta), \\ & r = \{\cos(2\theta)+((\cos(2\theta))^{2}+a^{4}-1)^{1/2}\}^{1/2}, \\ & \theta = 2\pi(j-1)/N \end{aligned} \right. \end{eqnarray} $

(5.2)

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

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

$ \begin{eqnarray} \zeta_{j} = z_{i}+r_{z}|z_{i}-z_{i-1}|\exp^{{i}\{\arg(z_{i}-z_{i-1})+(\pi/2)\}}. \end{eqnarray} $

(5.3)

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

$ \begin{eqnarray} \zeta_{j} = z_{i}+r_{w}|z_{i}-z_{i-1}|\exp^{{i}\{\arg(z_{i}-z_{i-1})-(\pi/2)\}}, \end{eqnarray} $

(5.4)

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

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

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

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

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