栗雪娟，王丹

(西安建筑科技大學(xué) 理學(xué)院，陜西 西安 710055)

本文考慮的四階Cahn-Hilliard方程為

(1)

1958年,Cahn和Hilliard提出Cahn-Hilliard方程,該方程最早被用來(lái)描述在溫度降低時(shí)兩種均勻的混合物所發(fā)生的相分離現(xiàn)象。隨著學(xué)者對(duì)該方程的研究越來(lái)越深入,該方程的應(yīng)用也越來(lái)越廣泛,特別是在材料科學(xué)和物理學(xué)等領(lǐng)域中有廣泛的應(yīng)用[1-3]。

Cahn-Hilliard方程的數(shù)值解法目前已有很多研究,文獻(xiàn)[4]使用了全離散有限元方法,文獻(xiàn)[5]使用了一類二階穩(wěn)定的Crank-Nicolson/Adams-Bashforth離散化的一致性有限元逼近方法,文獻(xiàn)[6-7]使用了有限元方法,文獻(xiàn)[8]使用了不連續(xù)伽遼金有限元方法,文獻(xiàn)[9]使用了Cahn-Hilliard方程的完全離散譜格式,文獻(xiàn)[10]使用了高階超緊致有限差分方法,文獻(xiàn)[11]使用了高階優(yōu)化組合型緊致有限差分方法。

綜上所述,本文擬對(duì)Cahn-Hilliard方程構(gòu)造一種新的超緊致差分格式,將空間組合型超緊致差分方法和修正的時(shí)間四階Runge-Kutta方法相結(jié)合,求解Cahn-Hilliard方程的數(shù)值解,得到相對(duì)于現(xiàn)有廣義格式精度更高的數(shù)值求解格式,并對(duì)組合型超緊致差分格式進(jìn)行誤差估計(jì),最后通過(guò)數(shù)值算例驗(yàn)證該方法的可行性。

1 高階精度數(shù)值求解方法

1.1 空間組合型超緊致差分格式

早期的緊致差分格式是在Hermite多項(xiàng)式的基礎(chǔ)上構(gòu)造而來(lái)的,Hermite多項(xiàng)式中連續(xù)三個(gè)節(jié)點(diǎn)的一階導(dǎo)數(shù)、二階導(dǎo)數(shù)和函數(shù)值的數(shù)值關(guān)系可以表示為

(2)

1998年,Krishnan提出如下緊致差分格式:

(3)

式(2)對(duì)應(yīng)f(x)展開(kāi)以xi為鄰域的泰勒級(jí)數(shù)為

(4)

差分格式的各項(xiàng)系數(shù)由式(3)決定,可得到如下的三點(diǎn)六階超緊致差分格式:

(5)

為優(yōu)化三點(diǎn)六階緊致差分格式,并保持較好的數(shù)值頻散,將迎風(fēng)機(jī)制[12]引入式(5),構(gòu)造出如下三點(diǎn)五階迎風(fēng)型超緊致差分格式:

(6)

左右邊界可達(dá)到三階精度緊致格式:

(7)

(8)

上述組合型超緊致差分格式只需要相鄰的三個(gè)節(jié)點(diǎn)便可以同時(shí)求得一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的五階精度近似值,比普通差分格式的節(jié)點(diǎn)更少,降低了計(jì)算量。

為便于編程計(jì)算,將上述構(gòu)造的組合型超緊致差分格式重寫為矩陣表達(dá)形式。假設(shè)U為位移矩陣,其大小為m×n,則求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的離散過(guò)程可以用矩陣運(yùn)算表示為

AF=BU,

(9)

結(jié)合內(nèi)點(diǎn)的三點(diǎn)五階迎風(fēng)型超緊致差分格式和邊界點(diǎn)的三點(diǎn)三階差分格式,組成式(9)中等式左邊的矩陣A和等式右邊的矩陣B,大小分別為2m×2n和2m×n;F為奇數(shù)行為空間一階導(dǎo)數(shù)和偶數(shù)行為空間二階導(dǎo)數(shù)組成的矩陣,大小為2m×n。以上矩陣分別為:

(10)

(11)

(12)

(13)

由式(9)可得

F=A-1BU。

(14)

解線性代數(shù)方程組(9)可得Cahn-Hilliard方程的空間一階導(dǎo)數(shù)和二階導(dǎo)數(shù)。對(duì)于四階導(dǎo)數(shù),可將已求得的二階導(dǎo)數(shù)替代式(14)中的U,再次使用式(14)進(jìn)行求取。

1.2 時(shí)間離散格式

在對(duì)很多偏微分方程的數(shù)值求解中不僅需要高精度的空間離散格式,同時(shí)還需要高精度的時(shí)間離散格式。普通的一階精度時(shí)間離散格式顯然滿足不了高精度計(jì)算要求,因此本文選用時(shí)間四階Runge-Kutta格式進(jìn)行時(shí)間離散。Runge-Kutta方法是基于歐拉方法改進(jìn)后的求解偏微分方程的常用方法,這種方法不僅計(jì)算效率高,而且穩(wěn)定性好。格式的推算過(guò)程如下:

假設(shè)求解方程為

(15)

式中F是對(duì)空間變量的微分算子,則修正的四階Runge-Kutta格式為

(16)

1.3 誤差估計(jì)

(17)

將式(17)代入式(6),所求得組合型超緊致差分格式的一階導(dǎo)數(shù)及二階導(dǎo)數(shù)對(duì)應(yīng)的截?cái)嗾`差為:

(18)

(19)

使用組合型超緊致差分格式的好處是在每一個(gè)網(wǎng)格點(diǎn)上存在一個(gè)一階和二階連續(xù)導(dǎo)數(shù)的多項(xiàng)式。本文比較了組合型超緊致差分格式和現(xiàn)有廣義格式的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的截?cái)嗾`差:

(20)

式中參數(shù)α,β,a,b,c在各種格式中取不同的值(表1,表2)。本文發(fā)現(xiàn)在各種方案中,組合型超緊致差分格式的截?cái)嗾`差最小。

表1 不同格式一階導(dǎo)數(shù)的截?cái)嗾`差

表2 不同格式二階導(dǎo)數(shù)的截?cái)嗾`差

2 數(shù)值算例

誤差范數(shù)L1和L2的定義為:

對(duì)給出的數(shù)值算例,計(jì)算誤差范數(shù)L1和L2,并采用四種方法進(jìn)行數(shù)值模擬,對(duì)其數(shù)值結(jié)果進(jìn)行誤差分析和對(duì)比,結(jié)果見(jiàn)表3,本文所使用方法效果最佳,由此證明所提方法的有效性和可行性。

表3 0.5 s時(shí)刻精確度測(cè)試結(jié)果(N=10)

用本文提出的式(6)—式(8)和式(16)計(jì)算算例, 圖1—圖3給出了不同時(shí)刻數(shù)值解與精確解的對(duì)比圖,可以看出,數(shù)值解與精確解吻合很好,表明本文給出的數(shù)值格式是可行的,并且精度較高。

(a)精確解(b)數(shù)值解圖2 0.5 s的精確解與數(shù)值解

(a)精確解(b)數(shù)值解圖3 1 s的精確解與數(shù)值解

3 結(jié)論

本文研究了組合型超緊致差分方法和四階Runge-Kutta方法,并將其運(yùn)用于四階Cahn-Hilliard方程的數(shù)值求解,通過(guò)研究與分析,得到如下結(jié)論:

1)使用泰勒級(jí)數(shù)展開(kāi)鎖定差分格式系數(shù),得到本文的組合型超緊致差分格式精度更高,誤差更小。

2)在邊界點(diǎn)處有效地達(dá)到了降階,并提高了精度。

3)通過(guò)數(shù)值算例驗(yàn)證了數(shù)值格式的有效性。

4)預(yù)估該方法可應(yīng)用于高階偏微分方程的數(shù)值求解。