一類新的求解雙邊障礙問題的數(shù)值方法

冉清華 程曉良

(1.浙江大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,浙江杭州 310000;2.貴州大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院,貴州貴陽 550000)

§1 引言

障礙問題是典型的第一類變分不等式,它出現(xiàn)在物理,工程,金融和管理科學(xué)等領(lǐng)域并有著廣泛的應(yīng)用,比如潤滑理論中的彈塑性扭轉(zhuǎn)問題,彈性理論中的薄膜問題,尾流問題以及期權(quán)定價等.

假設(shè)Ω是R2上的一個具有Lipschitz邊界?Ω的有界域.給定函數(shù)f ∈L2(Ω),φ,ψ ∈H1(Ω)滿足φ ≤ψ,φ|?Ω ≤0及ψ|?Ω ≥0.令在Ω中幾乎處處成立},雙邊障礙問題描述為:求u ∈K,使得

對該問題的物理解釋為:位于區(qū)域Ω且邊界固定在?Ω的彈性膜,在垂直力τf(τ為彈性膜的張量系數(shù))的作用下的平衡位置u(膜的垂直位移分量)滿足式(1),其中φ,ψ分別為下障礙和上障礙的高度.問題(1)存在唯一解.根據(jù)能量最小原理,解u可以用變分不等式形式來刻畫[2]:

進(jìn)一步,假設(shè)f ∈C(Ω),ψ ∈C(Ω),φ ∈C(Ω)和u ∈C2(Ω)∩C()，障礙問題又可表述為自由邊界的偏微分方程的形式:

其中,Ω0={x ∈Ω|φ ＜u ＜ψ}表示膜與障礙不接觸的區(qū)域,Ω1={x ∈Ω|u=φ}和Ω2={x ∈Ω|u=ψ}分別表示膜與下障礙,上障礙接觸的部分,滿足Ω=Ω0∪Ω1∪Ω2.由于Ωi(i=0,1,2)是未知的,(3)屬于非線性問題.

作為一類經(jīng)典的問題,障礙問題受到了許多研究者的關(guān)注,下面簡單回顧一些已有算法.投影松弛法[3-4]是最著名的方法.它們具有收斂,易于實(shí)現(xiàn)的特點(diǎn).然而,其收斂速度依賴于松弛參數(shù)的選擇和與網(wǎng)格尺寸相關(guān)的自由度的數(shù)目[5].有效集策略[5-6]是一種求解離散障礙問題的有效方法,其基本迭代包括兩個步驟.首先,根據(jù)某個指定準(zhǔn)則將離散區(qū)域分解為活動部分和非活動部分.然后,在非活動集部分求解簡化的線性系統(tǒng).多重網(wǎng)格算法[7-9]可以使用迭代算法構(gòu)造合理的域劃分,其收斂速度不依賴于網(wǎng)格尺寸,但需要在某些子域上采用非線性求解算法.“活動障礙”法[10],分片線性迭代算法[11-12],罰方法[13],區(qū)域分解法[14]等也相繼被提出,具有較好的效果.基于離散障礙問題的互補(bǔ)形式,本文通過引入一個參變量將原不易求解的互補(bǔ)問題等價改寫為易于數(shù)值實(shí)現(xiàn)的等式方程組,并利用牛頓迭代法進(jìn)行數(shù)值求解.

論文余下部分的結(jié)構(gòu)如下:§2給出了雙邊障礙問題的等價描述.§3將所得的等價系統(tǒng)正則化并給出了迭代算法.§4數(shù)值實(shí)驗(yàn)說明了本文算法的有效性.最后,§5對全文進(jìn)行了總結(jié).

§2 雙邊障礙問題的等價描述

本文考慮障礙問題的數(shù)值求解,故首先將問題(1)離散,得到有限維障礙問題:

其中

和fff=(f1,f2,……,fm)T分別對應(yīng)φ(x),ψ(x) 和f(x)的離散形式,(·,·)表示歐式空間Rm中的內(nèi)積.剛度矩陣A=(aij)∈Rm×m是對稱正定的.記I={1,2,……,m}為所有離散節(jié)點(diǎn)構(gòu)成的集合.類似于連續(xù)情況,問題(4)具有唯一解,并且可以得到如下等價的離散障礙問題[3-6,12].

引理2.1下面三個條件等價.

顯然,離散障礙問題的上述三種等價形式分別對應(yīng)于能量極小化問題(1),變分不等式(2),自由邊界的偏微分方程(3)的離散.關(guān)于上述離散結(jié)果的收斂性可類似于[3]中單邊障礙問題的討論進(jìn)行分析,這里不給出細(xì)節(jié).本文著重介紹離散系統(tǒng)(6)的求解.對于離散障礙問題(6)的數(shù)值求解,已有很多算法被提出.例如,活動集策略是求解離散障礙問題的有效方法,但是為了使迭代過程有效,需要應(yīng)用適用于非正則幾何區(qū)域的快速線性系統(tǒng)求解器,在大規(guī)模計算中收斂較慢.[12]首先將離散障礙問題轉(zhuǎn)化為分片線性系統(tǒng)

然后應(yīng)用皮卡迭代算法求解得到uuu=min(ψ,max(zzz,φ)).雖然該算法能在有限步收斂到離散障礙問題的解,但是其收斂速度較慢,而且計算過程復(fù)雜.基于此,本文提出如下計算方法.

首先介紹兩個輔助函數(shù).對于i=1,2,……,m,定義函數(shù)gi:θ ∈R→gi(θ)∈R:

和pi:θ ∈R→pi(θ)∈R:

為了更直觀的表達(dá),圖1給出了gi(θ)和pi(θ)的圖像.

圖1 gi(θ)和pi(θ)的圖像

現(xiàn)在給出本文的主要結(jié)果.引入變量Θ=(θ1,θ2,……,θm)T∈Rm,構(gòu)造如下等式方程組:

其中,ggg:Θ ∈Rm →ggg(Θ)=(g1(θ1),g2(θ2),……,gm(θm))T∈Rm,ppp:Θ ∈Rm →ppp(Θ)=(p1(θ1),p2(θ2),……,pm(θm))T∈Rm.下面證明所構(gòu)造的系統(tǒng)(9)與離散障礙問題(6)是等價的.

定理2.2存在Θ=(θ1,θ2,……,θm)T∈Rm使得系統(tǒng)(9)與問題(6)等價.

證 (?) 假設(shè)uuu是(9)的解.根據(jù)(9)(b)和ppp(Θ)的定義可得φi ≤ui ≤ψi,i=1,2,……,m.下面進(jìn)行分類討論.

若φi ＜ui ＜ψi即,則由(9)(b)可得φi ＜pi(θi)＜ψi.根據(jù)pi(·)的定義可知φi ＜θi ＜ψi.從而gi(θi)=0.因此,(Auuu)i-fi=0.

若ui=φi即,通過(9)(b)和(8)可以得到θi ≤φi.因此

若ui=ψi即,則由(9)(b)可得pi(θi)=ψi.故θi ≥ψi.從而

綜上可得,uuu為問題(6)的解.

(?) 假設(shè)uuu是(6)的解,下證存在Θ=(θ1,θ2,……,θm)T∈Rm使得uuu滿足系統(tǒng)(9).

若,則(Auuu)i-fi=0且φi ＜ui ＜ψi.取θi=ui,則成立

若,則(Auuu)i-fi ≥0且ui=φi.根據(jù)gi(·),pi(·)的定義,顯然存在θi ≤φi使得

而且,由gi(·)的單調(diào)性可知θi是唯一存在的.

若,(Auuu)i-fi ≤0且ui=ψi.類似地,存在唯一的θi ≥ψi使得

綜上所述,存在Θ=(θ1,θ2,……,θm)T∈Rm使得(uuu,Θ)滿足(9).

顯然從定理的證明過程可以得到如下結(jié)論.

定理2.3系統(tǒng)(9)存在唯一解(uuu,Θ).

本節(jié)通過引入一個參數(shù)Θ,將原非線性的不等式問題(6)等價轉(zhuǎn)化為等式方程組(9),使得問題更容易計算.由于(9)是一個非線性系統(tǒng),可以考慮用牛頓迭代法進(jìn)行求解.

§3 正則化及迭代算法

眾所周知,牛頓迭代法的最大優(yōu)點(diǎn)是收斂速度快,其在單根附近具有二階收斂速度.牛頓迭代法也因此而被廣泛使用.下面,本文利用該方法計算系統(tǒng)(9).

由于系統(tǒng)(9)中的gi(θ)和pi(θ)不可導(dǎo),首先對其進(jìn)行正則化處理.分別定義gi(θ)及pi(θ)的C1近似gic(θ),pic(θ):

其中常數(shù)c ＞0.gic(θ),pic(θ) 的函數(shù)圖像如圖2所示.不難證明,當(dāng)c →∞時,

圖2 gic(θ)和pic(θ)的圖像

因此得到系統(tǒng)(9)的正則化結(jié)果:

其中g(shù)ggc(Θ)=(g1c(θ1),g2c(θ2),……,gmc(θm))T,pppc(Θ)=(p1c(θ1),p2c(θ2),……,pmc(θm))T.利用牛頓迭代法的思想,得到關(guān)于(uuu,Θ)的如下線性系統(tǒng):

其中

從而將原問題轉(zhuǎn)化為如下代數(shù)系統(tǒng):

其中,E是Rm×m中的單位矩陣,

基于上述分析,給出求解障礙問題的如下算法步驟.

算法

第1步 輸入?yún)?shù):正則化參數(shù)c,網(wǎng)格尺寸h,誤差精度∈0,初始值uuu0,Θ0,n ←0.

第2步 根據(jù)(16)求解uuun+1,Θn+1.

第3步 如果‖Θn+1-Θn‖＜∈0,停止迭代并輸出uuun+1;否則n ←n+1,轉(zhuǎn)第2步.

注通過系統(tǒng)(13)不難發(fā)現(xiàn),當(dāng)Θn+1=Θn時,(uuun+1,Θn+1)是問題(12)的精確解.故算法中將‖Θn+1-Θn‖設(shè)置為停止準(zhǔn)則是合理的.

§4 數(shù)值模擬

圖3 例1的數(shù)值結(jié)果

圖4 例2的數(shù)值結(jié)果

圖5 例3的數(shù)值結(jié)果

在本節(jié)中,通過以下四個例子來說明本文所提出的算法的有效性.所有實(shí)驗(yàn)均在系統(tǒng)配置為8G內(nèi)存,2.3GHz Intel Core i5處理器,Matlab版本為R2018b的個人計算機(jī)上進(jìn)行.以下各例中,Ω=[0,1]×[0,1],并采用絕對誤差en+1=‖Θn+1-Θn ‖2＜∈0=10-6作為迭代停止準(zhǔn)則,其中‖ · ‖2為歐式范數(shù).除非特別說明,所得實(shí)驗(yàn)結(jié)果均在c=108,h=1/64,Θ0=uuu0=Afff的情況下獲得.

例1f(x,y)=-8,φ(x,y)=-dist((x,y),?Ω),ψ(x,y)=dist((x,y),?Ω).[13]

例2f(x,y)=11(x+y-1),φ(x,y)=-dist((x,y),?Ω),ψ(x,y)=dist((x,y),?Ω).[13]

例3[5,15]φ(x,y)=-dist((x,y),?Ω),ψ(x,y)=0.2,

例4f(x,y)=45(x-x2)(sin(11πx)+1),ψ(x,y)=dist((x,y),?Ω),φ(x,y)=-∞.[5]

以上各例的數(shù)值結(jié)果如圖3-圖6所示.圖中上接觸集和下接觸集的白色區(qū)域分別表示數(shù)值解與上障礙和下障礙接觸的區(qū)域.例1的數(shù)值解與上障礙沒有接觸.例4的數(shù)值解與下障礙沒有接觸.各例的迭代次數(shù)n及誤差en隨網(wǎng)格尺寸h變化的情況如表1所示.圖3-6以及表1的數(shù)據(jù)說明,本文算法能有效,快速地計算出數(shù)值解,而且收斂速度不依賴于網(wǎng)格尺寸.

圖6 例4的數(shù)值結(jié)果

表1 迭代次數(shù)及誤差隨網(wǎng)格尺寸變化的情況

§5 結(jié)語

本文通過引入一個參數(shù),將離散雙邊障礙問題的互補(bǔ)形式等價地改寫為易于求解的方程組,并將其正則化處理之后,運(yùn)用牛頓迭代法進(jìn)行數(shù)值求解.數(shù)值結(jié)果顯示該方法能快速,準(zhǔn)確地計算出數(shù)值解及接觸集,其收斂速度不依賴于網(wǎng)格尺寸的選取.