非線性發展方程的無網格比高精度有限元方法

石東洋,王俊俊

(1.鄭州大學數學與統計學院,河南鄭州 450001)(2.平頂山學院數學與統計學院,河南平頂山 467000)

1 前言

眾所周知,非線性發展方程的解通常無法直接用解析式寫出來,或是寫出來的表達式非常復雜,所以利用數值方法給出其近似解就顯得尤為重要.而對于有限元方法這一主流方向,我們常見的線性化BE(Backward-Euler)方法和CN(Crank-Nicolson)方法憑借可以避免在每一個時間層都要求解非線性方程的劣勢且不降低計算精度的優勢,成為了該方向的研究熱點之一.事實上,研究一個非線性發展方程的線性化有限元方法總會涉及到一個有限元解關于某種模的有界性問題,由于這些模的先驗估計不容易直接得到,通常的處理技巧就是利用逆不等式.比如在二維的情況下,考慮有限元解有界時的經典做法是

其中un是原始問題的解,Ih是某個插值算子或者投影算子.由于通常有誤差估計(m1,m2為某些正數),要使有界就不可避免的要對時間步長τ有一個限制,從而導致空間網格參數h與時間步長τ需要滿足某個比值關系(即網格比).在實際計算中,這樣的網格比經常會導致時間步長變的非常小,從而引起很長的耗時.因此,怎樣甩掉這些限制就成了備受關注的課題.最近,為了克服這一嚴重缺陷,孫偉偉、李步揚、王冀魯、高華東等學者都在此方面做出了許多有價值的工作.其主要思想(見文獻[1])是通過引入一個時間離散方程系統,并利用其解Un把誤差分裂成兩部分——時間誤差un?Un和空間誤差,利用時間誤差的結果得到關于時間離散方程解的正則性,再利用空間誤差得到有限元解的無網格比有界性事實上,由于空間誤差的分析過程中甩掉了經典誤差估計中的截斷誤差項,只要空間上的誤差能寫成(m3,m4為某些正數)的形式,網格比即可去掉.隨后,王冀魯、高華東、司志勇等又把該思想應用于非線性多孔介質流問題[2,3],非線性的Joule Heating方程[4],非線性Thermistor方程[5,6],非線性Schrdinger方程[7,8]和非線性Navier-Stokes方程[9]等.以上研究都考慮了這些非線性發展方程在協調元下關于無網格比的收斂性,有許多問題需要進行更深層次的研究.

首先,為了提高有限元解的逼近精度,超收斂的思想已成為了一個重要的研究途徑.事實上,在理論分析和實際計算中,若有好的網格,有限元解與有限元插值的誤差在某種范數的意義下比有限元解與真解的誤差要小得多,即超逼近現象.從上個世紀80年代開始,以林群院士為代表的眾多學者專家都在此方面取得了許多有出色的成果,所以如何將無網格比收斂的結果推廣到無網格比超收斂上去是我們感興趣的話題.但是,為了達到超收斂的結果,如果我們把文獻[1–9]中所考慮的區域換成更具一般性的矩形區域(不再滿足C2的條件),則由橢圓的正則性可以看到,引入的時間離散方程解的有界性就很難達到H3-模.因此,如何在時間離散方程解的有界性較弱的前提下,探討無網格比的超逼近結果就顯得尤為重要.

其次,由于非協調元方法在大多數情況下對方程解的正則性要求比較低,因此人們對非協調元的研究一直保持著較高的熱度(見以石鐘慈院士為代表的眾多學者專家所得到的具有特色的工作).這樣就很有必要研究怎樣利用自由度少、精度高的低階非協調單元研究非線性發展方程的無網格比的超收斂性.

再次,傳統的有限元方法對解的光滑度要求都比較高,這會給實際計算造成很多困難,因此混合有限元方法受到了高度的關注.事實上,混合有限元方法的關鍵性問題在于如何構造出合適的空間對,使其滿足LBB條件,這其實是不容易做到的.因此構造特別的格式來降低LBB條件的難度成為了一個熱點,比如:文獻[10–13]對二階橢圓問題提出了一種混合元格式,它具有當兩個逼近空間滿足一個簡單的包含關系時,LBB條件容易滿足且能避勉因涉及散度算子帶來的麻煩等優點.另一方面,直接繞開LBB穩定性條件(如最小二乘法、穩定化有限元方法等)也成為了大家另一個關注的方面.事實上,1998年,Pani在文獻[14]中提出了一種稱之為H1-Galerkin混合有限元方法．這種方法不需要所選取的混合元空間滿足LBB相容性條件,并被廣泛應用在各種方程上.例如,長波方程[15],雙曲方程[16],帶有記憶項的方程[17],積分微分方程[18–20],拋物方程[21].因此,怎樣利用H1-Galerkin方法得到非線性發展方程無網格比的超收斂結果是值得深思的.

最后,對于線性化的全離散格式來說,由于當時時刻的時間層分析需要用到上一時刻時間層的結論,我們往往會選擇數學歸納法進行證明.但是對于每一個時間層的結果到最后都應該由一個統一的系數來控制這個問題顯然就不是一件容易的事了.更進一步地,根據不同非線性問題的具體特點,針對不同方程的逼近格式,設計新的高效的有限元數值算法來驗證理論分析的正確性也是必須且困難的.

最近,我們在文獻[22–34]中,在分裂思想的基礎上,博采眾家之長,創新性的把無網格比、高精度分析與非協調單元、新的線性化離散格式等特色和優勢有機結合起來,形成非線性發展方程在全離散格式上無網格比約束的有限元超收斂分析的一套新理論體系.與此同時,更是嘗試著探索、研究一些特殊的方程,考慮繞過分裂的方法也達到無網格比的超收斂結果.

近期所做的工作,主要的創新點集中表現在以下幾個方面:

(1)超收斂結果對方程解的光滑性要求比較高,但構造時間離散輔助問題(即時間離散方程系統)時,在多邊形區域(例如矩形)下,就無法保證其解較強模的有界性.因此我們利用了一些特殊的、不同以往的技巧,在其解空間較弱的條件下得到無網格比超收斂的結論;巧用Taylor展開式對非線性項進行處理,以保證對時間步長τ的階不丟失.

(2)由于選擇的全離散格式是線性化的形式,在利用數學歸納法分析第n層的結果時需要用到第n?1層的結論,我們用一個統一的系數來控制每一個時間層的結果,這也是其數學歸納法成立的關鍵所在.

(3)構造了非線性雙曲方程新的二階格式,以此得到無網格比超收斂結果.而以往對非線性雙曲方程的無網格比研究甚至連收斂性也沒有見到報道.

(4)對一些特殊的非線性發展方程,拋棄分裂誤差思想,采用一些新的技巧也證明了其無網格比超收斂性.

本文的目的是在前期我們所做的工作的基礎上,挑揀出有特色的創新點給予說明,以期窺探出對非線性發展方程無網格比超收斂分析的重要方法和思路,起到拋磚引玉的作用.

2 非線性拋物方程的無網格比超收斂分析

非線性拋物方程有著深刻的物理背景,它的有限元方法也越來越受人們關注.例如:文獻[35]針對一般的非線性拋物方程建立了兩種線性化的格式,當τ≤h時,利用線性三角形元得到了L2-模意義下的收斂結果.文獻[36]在限制下利用兩層時間離散的方法討論了非線性拋物方程的最優誤差估計.文獻[37]采用了一個非線性H1投影,當τ4=O(hq),q≤3時,得到了其解在L2-模和H1-模意義下的最優誤差估計.文獻[1]利用分裂技巧擺脫了此類限制,給出了一類稱之為Joule Heating的非線性拋物型方程的協調元無網格比收斂性分析.我們看到,一方面,一般的非線性拋物方程中的非線性項?·(a(u)?u)的處理對于超收斂的分析是很有挑戰的,特別是在分析空間誤差的時候,怎樣處理a(u)才可以在不降階的情況下使其結果能提出空間網格參數h,才能在使用逆不等式的時候不產生網格比,但同時還得能維持數學歸納法所需要要的系數統一性?另一方面,當非線性拋物方程右端的非線性項是局部Lipschitz連續時,對有限元解的正則性要求可能會更苛刻,如何把這些限制考慮進去且得到無網格比超收斂結果是我們想要研究的方向之一.

考慮如下非線性拋物方程

其中??R2是一個矩形,其邊界為??,0

2.1 非協調有限元方法

首先,設f(u)是?上整體Lipschitz連續的函數,?是一個四條邊都平行于坐標軸的矩形,Γh是一個擬一致正則矩形剖分.對于給定的K∈Γh,令其四個頂點和四條邊分別為ai,i=1～4和.記.定義非協調有限元空間:

其中[vh]表示vh跨過單元邊界F的跳度,而當F???時,[vh]=vh.令Ih:H1(?)→Vh為相對應的插值算子,且Ih=Ih|K滿足

則文獻[21,38–40]證明了下面重要引理.

引理1若,則對于任意的vh∈Vh,有

這里?h表示分片梯度,且是Vh上的一個能量模.文獻[40]證明了對于任意正整數m,vh∈Vh,

設{tn:tn=nτ;0≤n≤N}是[0,T]上的一個等距剖分,時間步長為τ=T/N,設,且u(X,tn)=un,若為一列函數.記

利用這些記號,考慮(2.1)式的線性化Galerkin有限元逼近:尋找,使得對于任意的vh∈Vh0,

第一步 建立一個時間離散系統,當n>1時求Un滿足

當n=1時,利用以下式子計算U1:

和

其中U1,0(X)|??=0,U1(X)|??=0.

令e1,0,u1?U1,0,en,un?Un(n=0,1,2,···,N). 通過分析時間誤差,給出U1,0,Un(n=0,1,2,···,N) 的正則性. 設u和Um(m=0,1,2,···,N) 分別為 (2.1)和(2.9)–(2.11)式的解,u∈L∞(0,T;H3(?)),ut,utt∈L∞(0,T;H2(?)),uttt∈L∞(0,T;L2(?)),則對于m=1,···,N,存在τ0>0,使得當τ≤τ0時,有

和

其中C0是一個與m,h和τ無關的正數.

此時,注意到由于?是矩形,其邊界不屬于C1,那么就不容易得到Un的H3-模有界性.因此隨后的無網格比超收斂分析需要利用新的方法得到.

第二步討論空間誤差,也為最終無網格比超逼近結果kIhun?做好準備.給出記號

第三步令

在這個過程中注意到,將τ從內積的一端轉向另一端的恒等變化

化簡過之后需要估計誤差

如果按照傳統的方法,則有

這樣最終的結果會降一階.但是若利用Taylor展開則有

其中

這樣就可以保持到想要的結果.

2.2 協調有限元方法

限制f(u)為局部Lipschitz連續的,利用雙線性協調單元可研究(2.1)式的無網格比超逼近性質(區域及剖分如前面一樣).定義其有限元空間Vh0為

其中Ih:H2(?)→Vh0是相對應的插值算子,且對于以上雙線性元有以下高精度結果[41].

引理2若,則對于任意的vh∈Vh0,有

引入時間離散方程:當n≥1時求Un滿足

接下來,我們分別就時間誤差和空間誤差中的新技巧予以說明.

時間誤差記,則有

注意到在第n?1層中歸納假設(2.21)式成立后,進一步地得到是必需的.主要表現在以下兩個方面.

1.在第n層估計中,誤差方程左端有,右端部分在的前提下有估計項

和

可以看到,此時當τ充分小時,在第n層誤差方程的右端才可以去掉.

2.由于f的局部Lipschitz連續,要想估計誤差方程右端項,則必須得有做前提,這樣就保證了

空間誤差注意到,由于f的局部Lipschitz連續,在估計下面誤差時,

利用前面的結論以及協調元的性質,有

這樣可以估計得到

注1更進一步地,我們在文獻[34]中將文獻[10–12]中的混合元和無網格比的思想有機的結合起來,利用分裂內積等思想,得到了(2.1)式關于原始變量u的H1-模和~p=?u的L2-模的無網格比超收斂結果.

3 非線性Schrdinger方程的無網格比超收斂分析

其中?同(2.1)式,i是虛數單位,u0(X)是已知復值函數.另外,f(s)是一個實值函數,且關于s是二階可導連續的.

選用上一部分的區域,剖分和雙線性元單元,且仍定義其有限元空間為Vh0.令Rh:是定義在Vh0上的相對應的Ritz投影算子

有

且

更進一步地,當u∈H3(?),又由文獻[13]可知

其中Ih是定義在Vh0上相對應的插值算子.下面我們仍采用了分裂技巧,分別給出時間誤差和空間誤差上分析時所遇到的困難和解決方法.

時間誤差對于CN格式將誤差方程相鄰兩層相減,則

就變成

至此,得到的結果kenk2=O(τ2)相關估計,可以比文獻[8]的結果高二分之一階,也就是這樣的結果導出了,為后面的無網格比超逼近結果奠定了基礎.

對于BE格式得到結果kenk2=O(τ)比文獻[7]中的結果階高二分之一階,這也導出了,也為下面的空間誤差做出了鋪墊.

空間誤差1.對誤差方程相鄰兩層相減,對比文獻[45]中的結果,可以看到不需要的有界性也得到了無網格比高精度的結果,這就改進了已有結論.

2.使用插值算子和投影算子相結合的思想:若僅使用插值算子,為了得到高精度結果,避免不了利用高精度結果(?(un?Ihun),?vh)=O(h2)kunk3kvhk1或者(?(un?Ihun),?vh)=O(h2)kunk4kvhk0,則對Un和un的正則性要求過于苛刻.然而,在?為一個矩形的前提下,目前只能得到kUnk2的有界性,此時選用投影算子Rh是合適的.另一方面,若僅僅使用投影算子Rh,則不能構造相應于Rh的插值后處理算子,也就不能得到整體超收斂結果了.

4 非線性Sobolev方程H1-Galerkin混合有限元方法的無網格比超收斂分析

Sobolev方程起源于流體通過裂隙巖石的流動、二階流體的熱力學剪切和粘土的固結等物理現象.到目前為止,已有很多文獻研究了它的數值方法.例如:文獻[46]得到了當τ=O(hd/3)(d≤3)時,在三種情況下關于H1-模最優誤差估計結果.而文獻[47]利用混合有限元方法,在條件τ=O(h)下得到了最優誤差估計.文獻[38]控制條件為τ=O(h1+ε)(ε>0)時分別利用協調有限元和非協調有限元討論了其特征有限元方法,也得到了H1(?)-模和L2(?)-模的最優誤差估計.

大家都知道,H1-Galerkin方法是一個不需要滿足LBB條件的混合有限元方法,加上一些技巧的應用,還可得到關于流量~p=?u散度模的誤差估計.但是由于H1-Galerkin混合有限元方法需要的時間離散方程解的正則性較高,在矩形區域下不容易得到,所以對于一般的諸如非線性拋物方程利用上述分裂技巧直接處理暫時還不能去掉h和τ的比值.非線性Sobolev方程有著其自身的特點,它比非線性拋物方程多了一個非線性的導數項,正是多了這一項,使得我們考慮在分析的時候可以不用以上的分裂法就得到無網格比超收斂結果.因此我們通過與前面不同的分析,不引進時間離散方程,即在不必考慮所謂的時間誤差的前提下,避免由于時間離散方程解的正則性達不到相應的要求而帶來的麻煩,給出了無網格比的超逼近結論.

考慮如下非線性Sobolev方程:

這里對于正數b1,有|b(u)|≤b1,其余同(2.1)式中的假設.

當n=1時,

和

下面先給出一個新的引理.

引理3對于任意的,則有

圖1

其中l1,l3分別為K的下邊和上邊,l2,l4分別為K右邊和左邊,則有

由引理3,并利用數學歸納法,可分以下幾步分析說明非線性Sobolev方程的無網格比超收斂結果

第一步.

第二步利用,得到

利用Gronwall引理,當τ充分小時有

所以

第三步進一步地,

將(4.9)式代入(4.10)式,當τ充分小,利用Gronwalls引理有,再利用(4.9)式得到kξnkh≤Ch2+Cτ2.

這里強調以下三點

(1)如果直接估計kξnkh,對τ和h的比例限制將不可避免;

5 非線性雙曲方程

在物理上,雙曲方程是一類一直很受關注的偏微分方程,它可以用來描述聲波和電磁波的傳播等現象,其中也有很多文獻關注其非線性問題的有限元方法.例如,文獻[48]和[49]分析了非線性雙曲方程的全離散格式,其中文獻[48]討論了混合有限元方法,達到了最優誤差估計.文獻[49]利用Galerkin交替方向法討論了一類三維非線性雙曲方程,利用先驗估計的結果得到了誤差的H1(?)-模和L2(?)-模.但上述結果也都沒有擺脫h和τ的比值限制,在文獻[48]和[49]中分別需要假設條件τ=O(h),hr=O(τ)(1≤r≤k+1,k≥0)和τ=O(h2).因此,如何有效的對非線性雙曲方程展開無網格比的研究具有相當重要的科學價值.另一方面,在現有的參考文獻中對非線性雙曲方程的二階線性化格式討論的非常少,怎樣構造新的非線性雙曲方程的線性化全離散格式,使得其有更好的穩定性以及超收斂結果也值得進行深入的探討.

考慮如下非線性雙曲方程

關于方程的一些基本假設如同第二節.

我們將對(5.1)式創造性地構造一個新的線性化二階格式,技巧性地證明其截斷誤差的二階性質,給出其相對應的時間離散方程解的正則性,并由此得到在非協調單元下無網格比的超逼近結果.

則有

和

時間誤差引入時間離散方程,利用其解Un分裂誤差,通過估計時間誤差得到Un的正則性.令,則有

和

注意到,在這一節里針對非線性雙曲方程構造了一個新的線性化的二階格式,可以看到要證明其截斷誤差為O(τ2)是非常不容易的.另一方面,在得到誤差結果時,由于C0在第n層和第n+1層必須統一,則在估計誤差時,我們需要利用ēm+1≤τ(m≤n?1),而不是ēm+1≤C0τ2(m≤n?1)來得到結果.

空間誤差利用以上結果得到超逼近結果,令,則有

在此過程中,有以下幾點需要特別關注

2.(5.8)式的右端不能直接被估計成C1(h2+τ2),否則h和τ的比值將無法避免.事實上,也不能將其估計為C1h,因為利用數學歸納法,我們需要統一n層和n+1層的,所以利用了.

進一步地,在誤差估計過程中會出現以下項

如果利用類似前面的估計方法,將τ從內積的一端轉向另一端,則有結果,這時將不可避免的出現網格比.

為了克服關鍵問題,重新分裂內積為

則有

最后可導出以下結果

6 總結與展望

首先,把原來文獻中已有的對非線性發展方程無網格比收斂的結果,進一步延伸到對非線性發展方程無網格比的超逼近和超收斂的研究中.由于要想得到高精度的結果,對原始方程解正則性的要求往往都會比較高.那么怎樣繞過在矩形區域下,引入的時間離散方程解的有界性達不到H3(?)時,有技巧的得到無網格比超逼近結果是之前文獻所不曾見過的.

其次,對于非線性發展方程,其非線性項中a(u)的處理也是之前很少有的.我們創新性的利用Taylor展開式,保持了最后時間方向上的誤差不丟失階.

再次,在分析過程中發現,對于一些方程,若僅僅利用插值算子將會提高對方程解的正則性要求;僅僅利用投影算子,沒有辦法進行插值后處理,進而沒有辦法得到整體超收斂的性質.所以在對非線性Schrdinger方程的分析當中,利用了插值算子和投影算子相結合的處理方式得到了其無網格比超逼近和超收斂的結果.

而后,由于分裂誤差的做法使得引入時間離散方程解在矩形區域下正則性結果不是很好,而此時想要利用H1-Galerkin有限元方法得到無網格比超逼近結果還是非常困難的.對于非線Sobolev方程來說,我們發現了此類方程自身的特點,創新性地拋棄了分裂誤差的方法,利用一些特殊的技巧也得到了到了其無網格比超收斂結果.

最后,由于還沒有文獻對非線性雙曲方程的無網格比討論過,特別是對其二階線性化的全離散格式研究還很少有文獻進行報道.我們嘗試著構造一個新的線性化二階格式,證明了其截斷誤差的二階性質,并給出了無網格比超逼近性質.

對未來的進一步工作,我們提出以下幾點供參考

(1)各向異性有限元的研究是目前該領域獨具特色和挑戰性的又一前沿熱點和難點(參見文獻[50–53]).怎樣把原來已有的擬一致的正則剖分下的無網格比超收斂結果拓展到各向異性剖分上去.另一方面,如何將目前只對矩形區域進行的無網格比超收斂研究拓展到一般區域.

(2)二重網格算法是處理非線性問題的有效方法(參見文獻[54,55]),但在誤差估計中粗網格的網格參數或者細網格的網格參數都會不可避免的和時間步長產生一個比值.事實上,我們在文獻[22]中已經討論了半線性拋物方程二重網格算法的無網格比收斂性.那么怎樣才可以得到非線性發展方程二重網格算法的無網格比高精度分析呢?由于二重網格算法在粗網格上往往是給出一個非線性的逼近方程,要去掉網格比約束條件,便要每一層都引入一個非線性的時間離散方程.這樣,非線性系數的有界性在當時時間層是未知的,得到它的有界性便是非常困難的事了.另一方面,二重網格算法需要利用粗網格的結果得到在細網格上的逼近方程,使得它每一層都是一個線性問題,那么此時使用光滑邊界來提高此處引入的另一組時間離散方程系統解的正則性也許是可行的辦法,但這樣再處理超收斂分析就是一個挑戰了.

(3)很多數學物理問題都需要利用積分微分方程進行求解.數學、工程技術和自然科學的許多領域,例如在系統識別、流體力學、信號重構與圖像恢復、電磁場理論與靜電學、地球物理勘探等方面都將歸結為求解積分微分方程的問題.積分微分方程已經成為眾多學者研究的一份重要方向.更值得一提的是,對于我們比較難以處理的雙曲方程,也可以利用ut=q使其變成一個拋物型的積分微分方程.對于積分微分方程的無網格比收斂性研究,現在還沒有涉及到,尤其是該如何解決其積分項的估計都還是待解決的問題,因此這方面的工作還有待進一步地研究.

(4)怎樣將無網格比超收斂高精度結果進一步推廣應用到更為復雜的二階問題上去,比如說非定常不可壓縮Navier-Stokes方程?事實上,人們對于該方程的有限元方法做了大量研究,如協調有限元方法、全離散加罰有限元方法、多尺度方法、穩定化有限元方法等[56–58].關于無網格比,文獻[59]利用特征有限元方法繞過了Navier-Stokes方程的非線性項的繁瑣,再利用分裂技巧得到了無網格比收斂的結果.那么如果換成一般的有限元方法,處理它的非線性項使得無網格比結果不掉階是一個麻煩的事情.關于非線性四階問題[59–61](如EFK方程、Cahn-Hilliard方程、四階非線性雙曲方程等)的無網格比超收斂研究也是我們所關心的話題.事實上,對于非線性的Navier-Stokes或者是四階問題,混合有限元方法是常用的一個方法,那么怎樣利用它得到這些問題的無網格比超收斂結果也是我們未來研究的方向之一.

以上這些問題都是全新的、懸而未決的關鍵性問題,作為延伸和后續的工作,亟待建立一套新的框架性理論體系和分析模式去深入展開對這些問題的研究和探索.可以說,任何實質性的創新和突破都將進一步豐富有限元方法的內涵,提升其品味.