引 言
自由边界问题在工程和科学问题中有广泛的应用，如渗流问题、障碍问题等［1-3］.这类问题不同于一般的边值问题，因为自由边界问题在所考虑区域内部要满足微分不等式约束条件，从而具有很强的非线性，通常不能用解析方法求解，一般只能通过适当的数值方法获得近似解.利用变分不等式理论，这类问题解的存在唯一性已得到广泛的研究，采用有限差分法、有限元法和基本解方法等可将问题转化为有限维线性互补问题、二次规划问题或向量变分不等式等优化问题［4-8］，并可用Uzawa算法等迭代法进行数值求解，但这种方法的收敛速度对参数的取值非常敏感［2，7-8］.
Uzawa块松弛算法已广泛应用于求解各种等式或不等式约束条件的非线性问题.如Glowinski等针对Bingham流问题、弹塑性扭转问题等非线性问题，通过引入一个表示梯度函数的辅助变量，导出一个用增广Lagrange函数表示的鞍点问题，提出基于交替方向的Uzawa块松弛算法，用变分不等式理论证明了算法的收敛性［2，9］.在此基础上，Koko提出了求解单侧接触问题的Uzawa块松弛算法，包括无摩擦和摩擦问题，并给出了具体的算法过程［10］.Uzawa块松弛算法的主要优点是每次迭代均把所求非线性问题分解为两个线性子问题，迭代矩阵始终保持不变.另外，算法对罚参数具有全局收敛性.然而，该方法的罚参数取值对收敛速度影响非常大，罚参数太大或太小都将大大减缓收敛速度，因此在实际应用中面临如何优化选择罚参数的问题.
本文将Uzawa块松弛算法进一步推广到单侧障碍自由边界问题，并对罚参数的选择给出具体可行的办法.首先用差分法把单侧障碍自由边界问题离散为有限维线性互补问题，通过引入一个表示待求函数的辅助变量和增广Lagrange乘子法表示的鞍点问题，再用Uzawa块松弛算法求解，从而得到一个简单的两步迭代方法.该方法的每一次迭代均是先求解一个简单的线性问题，再利用对偶理论显式计算辅助变量.为了克服罚参数对收敛速度的影响，给出基于平衡原理灵活调整罚参数的自适应法则，这种法则已成功应用于投影算法［11-14］.本文在上述Uzawa块松弛算法基础上，进一步得到了类似的法则，自动选取使得算法收敛较快的可变罚参数，从而显著提高算法性能.最后用数值算例验证了算法的有效性.
1 自由边界问题的有限差分格式
考虑经典的二维单侧障碍自由边界问题，其数学模型的一般形式为
其中 Ω 是平面 Rn中的有界区域，其边界 Γ = Ω，x=(x1，x2)，L
，f(x)，g(x)和ψ(x)是已知函数，利用变分不等式理论，便可得到该问题解的存在唯一性［1-2］.
采用五点差分格式对问题(1)中的微分算子离散化:
其中向量函数u(x)表示未知函数v(x)在网格节点上的近似值.将已知区域Ω内的函数f(x)、障碍函数ψ(x)和边界Γ上的函数g(x)代入式(2)，则上述问题可改写为标准的N维线性互补问题:
其中A为一个N阶对称正定矩阵，N维列向量q的元素依赖于已知函数f(x)，ψ(x)及g(x).
定义非空有界闭凸集:
互补问题(3)等价于向量变分不等式
或极小值问题寻找N维列向量u满足
由于该极小值问题存在唯一解u*，因此线性互补问题(3)也存在唯一解u*.
2 增广Lagrange乘子法
为了求解上述有限维问题(3)的数值解，引入区域Ω上的辅助变量p∈RN，使得u=p.首先定义Lagrange函数
和相应的增广Lagrange函数
考虑式(4)与(5)对应的如下鞍点问题:
则有如下结果［2，9］.
定理1 如果 {u，p，λ}是L在RN×RN×RN上的一个鞍点，则对任意r≥0，{u，p，λ}也是Lr的鞍点，反之也成立.而且u也是问题(3)的解，且p=u.
证明(ⅰ） 因为 {u，p，λ}是L在RN×RN×RN上的一个鞍点，L(v，q，λ)∈R，则
由式(6)的前一部分以及式(4)可得到
由μ的任意性可得到
由式(6)的后一部分以及式(7)可得到
令q=v可得到
因为p=u，所以
由式(6)和(9)可得
因为
以上证明了 {u，p，λ}也是Lr的鞍点.
下面对上述证明做一个小结，由式(6)可知 {u，p}是下式的解:
在式(11)中令q=p得
在式(11)中令v=u得
（ⅱ）当r≥0时，{u，p，λ}是Lr的一个鞍点，则
利用（ⅰ）中的方法，同理可得u=p，且u也是问题(3)的解.
由式(14)可知 {u，p}是下式的解:
在式(15)中令q=p得
因为 u∈RN，对v∈RN，有u+t(v － u) ∈RN，即t v+(1 － t)u∈RN，其中t∈(0，1］.考虑 (r/2)(‖v－p‖2－‖u－p‖2)，令v=u+t(v－u)，代入其中得(rt/2)‖v－u‖2+ρ〈u－ p，v－ u〉，当 t→0时，结合式(16)得
在式(15)中令v=u得
因为p∈RN，对q∈RN，有p+t(q－p)∈RN，即t q+(1－t)p∈RN，其中t∈(0，1］.考虑(r/2)(‖u－q‖2－‖u－p‖2)，令q=p+t(q－p)，代入其中得(rt/2)‖q－p‖2+ρ〈p－ u，q － p〉，当 t→0时，结合式(18)得
因为 u=p，所以式(17)及(19)可归纳为式(12)及(13)，故 {u，p，λ} 是式(11)的解.再由式(9)可知，{u，p，λ}也是式(6)的解.综上，此定理得证.
3 Uzawa块松弛算法
根据定理1，要求解原问题(3)的解，只需求解增广Lagrange函数Lr的鞍点 {u，p，λ}，使得
由这个问题出发可得到问题(3)的Uzawa块松弛算法［9-10］，算法过程如下:
第1步 给定初始值，{p0，λ1}∈RN×RN，ρ＞ 0，置n=1.
第2步 求得u n∈RN，使得对v∈RN有
第3步 求得p n∈RN，使得对q∈RN有
第4步 更新Lagrange乘子:
第5步 对给定的某种判定条件，若满足，则停止迭代得到数值解u n;否则，置n=n+1，返回第2步.
接下来详细讨论子问题(20)和(21)的具体求解方法.
问题(20)等价于求解以下方程:
根据式(21)，能找到p n∈RN使得
且在Ω上有
因为式(24)和(25)是约束极小值问题，故可以用鞍点理论来显式地求解该问题.问题(21)等价于以下式子:
其中γn是式(25)的Lagrange乘子，且γn＞ 0.由式(26)可得到
将式(28)代入式(27)得到
又因为γn＞0，所以
即
故得到Lagrange乘子:
将式(29)代入式(28)得到
综合以上分析，得到如下Uzawa块松弛算法.
算法1(UBR)
第1步 给定初始值，{p0，λ1}∈RN×RN，r≥0，令n=1.
第2步 求解u n∈RN，使得对v n∈RN有
第3步 计算辅助变量:
第4步 更新Lagrange乘子:
第5步 对给定的误差限ε＞0，若满足
则迭代停止，否则返回第2步.
4 算法收敛性分析
为了证明算法收敛性，在此先给出以下引理［2］.
引理1 若AN×N:RN→RN是一个连续并且强单调的映射，{p n}n≥0是RN中的一个序列，对于任意的n，存在p∈RN，使得
则
从而，可得Uzawa块松弛法的收敛性结果［2］.
定理2 当
由算法产生的序列 {u n，p n，λn}收敛于Lr在R上的鞍点 {u，p，λ}.
证明 令 δu n=u n－ u，δp n=p n－p，δλn=λn－λ.因为 {u，p，λ}是 Lr在R 上的鞍点，从而有
再通过式(20)～(22)分别定义 {u n，p n，λn+1}，有
分别在式(30)和(33)中令v=u n和v=u，相加得
分别在式(31)和(34)中令q=p n和q=p，相加得
再将式(36)和(37)相加得
结合式(32)和(35)可得
进而得到
结合式(38)和(39)可得
因为
根据式(41)可得
在式(34)中用n－1代替n得到
分别在式(34)和(43)中令 q=p n－1和 q=p n，相加得
根据式(35)可得
结合式(32)和(45)可得
根据式(44)和(46)可得
根据式(42)和(47)可得
结合式(40)和(48)可得
即
其中 〈E'(u n) － E'(u)，δu n〉=〈Au n－ Au，u n－ u〉.
当ρ=r时，显然有
则根据引理1可知
当ρ≠r时，由文献［5］可得以上结果.
λ*是 {λn}n在RN上的一个弱聚点，在式(33)和(34)中取极限并且利用E的下半连续性可得
式(49)、(50)等价于以下式子:
因为p=u，所以
根据式(51)和(52)，可知 {u，p，λ*}也是 Lr的鞍点.
5 改进Uzawa块松弛算法
在应用Uzawa块松弛法求解实际问题时，通常取固定参数ρ=r［9-10］.然而，不同罚参数对收敛速度影响较大.因此，本文提出自适应Uzawa块松弛法，利用自适应法则，得到变参数序列{ρn}代替固定参数ρ［11-14］.使用如下规则来调整参数，在算法收敛性分析中，如果ρ=r=ρn，则由算法2生成的序列 {(u n，p n，λn)}满足以下不等式:
其中β＞0，可以得到
‖λn+1－ λ‖2+ ρ2n‖p n－ p‖2≤ ‖λn－ λ‖2+ ρ2n‖p n－1－ p‖2.
显然，序列 {‖λn+1－λ‖2+ρ2n‖p n－p‖2}单调递减且有界.为了加快收敛速度，我们希望‖λn+1－ λn‖ ≈ ρn‖p n－ p n－1‖ .
这提供了一个选择ρn+1的基本思路.对于一个给定的正常数τ，如果ρn‖p n－p n－1‖ ＞(1+τ)‖λn+1－ λn‖，就在下一次迭代中减小 ρn;如果 ρn‖p n－ p n－1‖ ＜ (1/(1+ τ))‖λn+1－λn‖，就在下一次迭代中增加 ρn.令 ωn= ‖λn+1－ λn‖/(ρn‖p n－ p n－1‖)，可以通过以下方式选择参数ρn+1:
综上，得到如下自适应Uzawa块松弛算法.
算法2(MUBR)
第1步 给定初始值，{p0，λ1}∈RN×RN，参数ρ＞ 0，τ ＞ 0，令ρ1=ρ，n=1.
第2步 求解u n∈RN，使得对v n∈RN有Au n+ ρn u n=－ q+ λn+ ρn p n－1.
第3步 计算辅助变量:
第4步 更新Lagrange乘子:
第6步 对给定的误差限ε＞0，若满足
则迭代停止;否则，置n=n+1，返回第2步.
6 算例分析
为了验证算法的可靠性，用本文提出的算法得到了两个算例的数值结果.对于Uzawa块松弛算法(算法1)，取固定参数ρ=r.两种算法的迭代终止条件均取ε=10－4.
算例 1 考虑在正方形区域 Ω =(－ 1.5，1.5) × (－ 1.5，1.5) 中的自由边界问题［15］，其中 f=－2，ψ=0，Dirichlet边界条件可由解析解
确定，r=
，从而边界条件可由解析解得到.采用自适应Uzawa块松弛算法求解，图1为取步长h=1/40和τ=1的数值结果，图2给出了自由边界的数值解与精确解结果对比，它们是吻合的.
图1 数值解
Fig.1 The numerical solution
图2 自由边界的数值解与精确解
Fig.2 Numerical and exact solutions of the free boundary
对于不同的初始罚参数ρ和不同的步长h，表1对算例1的Uzawa块松弛算法和自适应Uzawa块松弛算法所需的迭代次数进行了比较，其中“－”表示迭代次数超过500次.结果表明自适应Uzawa块松弛算法不仅收敛快，而且非常稳定.
表1 两种算法的迭代次数情况(算例1)
Table 1 The number of iterations for each method(example 1)
ρ algorithm 1(UBR)algorithm 2(MUBR)h=1/10 h=1/20 h=1/40 h=1/80h=1/10 h=1/20 h=1/40 h=1/80 10－2 － － － － 32 31 29 31 10－1 282 211 201 201 28 26 27 28 1 89 90 78 76 23 21 24 25 10 27 34 28 27 22 24 23 23 102 70 69 68 68 23 24 24 26 103 －－－－27 27 27 28 104 －－－－31 29 32 32
算例2 考虑在正方形区域Ω=(－2，2)×(－2，2)中的自由边界问题［6-7，16］，其中
解析解为
其中r=
，R=2，r*=0.697 965 142 8…，r* 满足条件 r*2(1 － ln(r*/R))=1.采用自适应Uzawa块松弛算法求解，同样取h=1/40和τ=1，图3为数值解结果，图4给出了自由边界的数值解与精确解结果对比，它们同样是吻合的.
对算例2，表2列举了两种算法对不同的初始罚参数ρ和不同的步长h所需的迭代次数，通过比较结果，同样表明自适应Uzawa块松弛算法是非常有效的.而且从两个算例的数值结果看出，自适应Uzawa块松弛算法的收敛速度几乎不受初始参数ρ和步长h的影响.
图3 数值解
Fig.3 The numerical solution
图4 自由边界的数值解与精确解
Fig.4 Numerical and exact solutions of the free boundary
表2 两种算法的迭代次数情况(算例2)
Table 2 The number of iterations for each method(example 2)
ρ algorithm 1(UBR)algorithm 2(MUBR)h=1/10 h=1/20 h=1/40 h=1/80h=1/10 h=1/20 h=1/40 h=1/80 10－2 － － － － 35 29 30 33 10－1 160 214 160 180 32 26 28 30 1 67 60 55 50 32 25 28 30 10 37 43 49 55 33 29 30 34 102 291 313 361 413 36 30 33 38 103 －－－－39 34 37 39 104 －－－－43 39 40 43
7 结 论
本文提出了求解单侧障碍自由边界问题的Uzawa块松弛法及其改进算法.Uzawa块松弛法的主要优点是每次迭代计算简单.然而，Uzawa块松弛算法中的迭代次数高度依赖罚参数ρ，为了提高算法的性能，本文提出基于自适应法则的自适应Uzawa块松弛法，该方法在迭代过程中自动调整罚参数ρ，数值算例表明自适应Uzawa块松弛算法具有更快的收敛速度和稳定性.
［1］ 韩渭敏，程晓良.变分不等式简介:基本理论、数值分析及应用［M］.北京:高等教育出版社，2007.(HAN Weimin，CHENG Xiaoliang.Introduction to Variational Inequality:Element Theory，Numerical Analysis and Applications［M］.Beijing:Higher Education Press，2007.(in Chinese))
［2］ GLOWINSKI R.Numerical Methods for Nonlinear Variational Problems［M］.Berlin:Springer-Verlag，2008.
［3］ 饶玲.单调迭代结合虚拟区域法求解非线性障碍问题［J］.应用数学和力学，2018，39(4):485-492.(RAO Ling.Monotone iterations combined with fictitious domain methods for numerical solution of nonlinear obstacle problems［J］.Applied Mathematics and Mechanics，2018，39(4):485-492.(in Chinese))
［4］ LIN Y，CRYER C W.An alternating direction implicit algorithm for the solution of linear complementarity problems arising from free boundary problems［J］.Applied Mathematics and Optimization，1985，13(1):1-17.
［5］ BURMAN E，HANSBO P，LARSON M G，et al.Galerkin least squares finite element method for the obstacle problem［J］.Computer Methods in Applied Mechanics and Engineering，2017，313:362-374.
［6］ LI X，YUAN D M.Asymptotic approximation method for elliptic variational inequality of first kind［J］.Applied Mathematics and Mechanics(English Edition)，2014，35(3):381-390.
［7］ YUAN D M，CHENG X L.A meshless method for solving the free boundary problem associated with unilateral obstacle［J］.International Journal of Computer Mathematics，2012，89(1):90-97.
［8］ 王光辉，王烈衡.基于对偶混合变分形式的Uzawa型算法［J］.应用数学和力学，2002，23(7):682-688.(WANG Guanghui，WANG Lieheng.Uzawa type algorithm based on dual mixed variational formulation［J］.Applied Mathematics and Mechanics，2002，23(7):682-688.(in Chinese))
［9］ GLOWINSKI R，TALLEC P L E.Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics［M］.Philadelphia:SIAM，1989.
［10］ KOKO J.Uzawa block relaxation method for the unilateral contact problem［J］.Journal of Computational and Applied Mathematics，2011，235(8):2343-2356.
［11］ HE B S，LIAO L Z，WANG S L.Self-adaptive operator splitting methods for monotone variational inequalities［J］.Applied Numerical Mathematics，2003，94(4):715-737.
［12］ 钟艳丽，严月月，张守贵.求解单侧障碍问题的自适应投影方法［J］.重庆师范大学学报(自然科学版)，2018，35(1):70-76.(ZHONG Yanli，YAN Yueyue，ZHANG Shougui.A self-adaptive projection method for the unilateral obstacle problem［J］.Journal of Chongqing Normal University(Natural Scienes)，2018，35(1):70-76.(in Chinese))
［13］ ZHANG S G.Projection and self-adaptive projection methods for the Signorini problem with the BEM［J］.Applied Mathematics and Computation，2017，74(6):1262-1273.
［14］ ZHANG S G，LI X L.A self-adaptive projection method for contact problems with the BEM［J］.Applied Mathematical Modelling，2018，55:145-159.
［15］ MARKUS B，SCHRODER A.A posteriori error control of hp-finite elements for variational inequalities of the first and second kind［J］.Computers ＆ Mathematics With Applications，2015，70(12):2783-2802.
［16］ ZOSSO D，OSTING B，XIA M，et al.An efficient primal-dual method for the obstacle problem［J］.Journal of Scientific Computing，2017，73(1):416-437.