李 娟
(南京审计大学金审学院 基础部, 江苏 南京 210023)
材料在介观和宏观尺度上的性质是由其复杂微观结构空缺、晶界和位错等拓扑缺陷决定的。晶体相场模型是描述材料微观结构中拓扑缺陷的一个有效方法[1-2]。对于固态晶体材料,其原子的位置呈规则周期性排列, 故考虑晶体相场模型的周期初边值问题:
φt=(φ3+(1-ε)φ+2φxx+φxxxx)xx,
x∈R,t>0
(1)
φ(x,0)=φ0(x),x∈R
(2)
其中,φ为场密度;ε<1为正常数。
上述方程是非线性高阶反应扩散方程,其精确解难以求出,因此发展其数值方法求数值解具有实际应用价值。文献[3]利用凸分解方法讨论了晶体相场模型的无条件稳定数值格式, 并给出了L2范数下的时间方向一阶、空间方向二阶的无条件收敛性;文献[4]利用凸分解方法给出了二阶非线性三层差分格式, 并利用非线性多重网格方法求解数值格式;文献[5]给出了一些二阶截断误差的能量稳定数值格式, 但未进行收敛性分析;文献[6]建立了无条件稳定的二阶差分格式, 但未讨论全局收敛性;文献[7]研究了新的一阶无条件梯度稳定、二阶无条件弱能量稳定凸分解方法, 并给出了数值格式的质量守恒性、唯一可解性、能量稳定性和截断误差的阶。这些研究均对晶体相场方程建立了稳定的具有二阶收敛精度的数值格式, 但非线性数值方法在数值计算过程中会因迭代运算而耗费大量的时间;文献[8]对该方程建立了三层线性化差分格式, 并证明了差分格式的唯一可解性及在L2范数下的二阶收敛性;文献[9]讨论了晶体相场模型的一阶、二阶无条件稳定的线性化格式, 但未对收敛性进行证明。目前对该问题数值方法的研究在收敛精度上仅在时间、空间方向上二阶收敛。
为提高数值收敛精度和数值计算速度,本文提出空间方向四阶、时间方向二阶的线性化差分方法,并证明其收敛性。利用泰勒公式得到求解第1时间层的显式差分格式,对其余时间层建立三层线性化隐式紧差分格式;并利用能量分析方法对差分格式进行理论分析,证明差分格式是唯一可解的,在L2范数下是收敛的,通过数值算例 验证了差分格式的有效性。
本节讨论问题(1)、(2)数值格式的推导。考虑周期边界条件, 不失一般性, 只讨论一个周期[0,L]内的数值方法。
取正整数M、N,时间区间[0,T]。记h=L/M,τ=T/N,Ωh={xi|xi=ih,0≤i≤M,i∈Z},Ωτ={tn|tn=nτ,0≤n≤N,n∈Z},Uh={u|u={ui},ui+M=ui,0≤i≤M,i∈Z}。对任意网格函数u∈Uh,记
15ui+1-20ui+15ui-1-6ui-2+ui-3)。
对网格函数u∈Uh, 定义如下平均值算子:
A2ui=A(Aui),A3ui=A(A2ui),
对定义在Ωτ上的网格函数w=(w0,w1,w2,…,wN),引入以下记号:
为便于差分格式的推导,引入下面的引理。
引理1[10]若g(x)∈C6[xi-1,xi+1], 则有:
其中,θi∈(-1,1);α(s)=(1-s)3[5-3(1-s)2]。
下面应用降阶法对问题(1)、(2)建立数值格式。令
v=2φ+φxx,μ=φ3+(1-ε)φ+vxx,
则问题(1)、(2)等价于以下耦合方程组:
φt=μxx,x∈[0,L],t∈[0,T]
(3)
μ=φ3+(1-ε)φ+vxx,
x∈[0,L],t∈[0,T]
(4)
v=2φ+φxx,x∈[0,L],t∈[0,T]
(5)
φ(x,0)=φ0(x),x∈[0,L]
(6)
1≤i≤M, 1≤n≤N-1
(7)
1≤i≤M, 1≤n≤N-1
(8)
1≤i≤M, 1≤n≤N-1
(9)
(10)
再用算子A2作用(7)式, 并将(10)式代入, 可得:
1≤i≤M, 1≤n≤N-1
(11)
其中
1≤i≤M, 1≤n≤N-1
(12)
本文对第1时间层建立显式差分格式。在方程(1)中令t=0, 结合(6)式, 可得:
2(φ0)xxxx+(φ0)xxxxxx
(13)
应用带积分型余项的泰勒公式,可得:
1≤i≤M
(14)
其中
(15)
由初值条件(6)可得:
(16)
若问题(1)、(2)的精确解适当光滑, 则存在正常数c0,使得:
(17)
1≤i≤M, 1≤n≤N-1
(18)
(19)
(20)
综上可知, 第0层数值解由初值条件(20)式给出, 第1时间层的数值解由显式格式(19)式求得。当第n、n-1层的数值解un、un-1(n≥1)已知时, (18)式为关于第n+1层数值解un+1的线性方程组,可通过解线性方程组求得其余时间层的数值解。
本节将对差分格式进行理论分析, 讨论唯一可解性和收敛性。对任意的u,v∈Uh, 引入如下内积和范数:
为便于差分格式唯一可解性及收敛性的证明,引入如下引理。
引理2[11]对于u,v∈Uh, 有
(u,δxv)=-(δxu,v),(u,Av)=(Au,v)。
由平均值算子A的定义及引理2知:
引理3 对于u∈Uh,有
下面利用能量分析法证明差分格式的唯一可解性。
定理1差分格式(18)~(20)式是唯一可解的。
证明用数学归纳法证明。由(19)式、(20)式易知,u0、u1已唯一确定。
假设第n-1、n层的数值解已唯一确定, 此时, (18)式为关于第n+1层数值解un+1的线性方程组。欲证其唯一可解性, 需证其对应的下列线性方程组仅有零解。
1≤i≤M, 1≤n≤N-1
(21)
用un+1与(21)式作内积, 可得:
(22)
再应用引理2,可得:
(23)
根据引理3, 由(23)式可得:
(24)
(25)
(26)
从而‖un+1‖2=0,即un+1是唯一确定的。
由数学归纳法知, 差分格式是唯一可解的。
下面讨论差分格式(18)~(20)式数值解的收敛性。定义误差函数为:
用(11)式、(14)式、(16)式分别减去(18)~(20)式,可得如下误差方程:
1≤i≤M, 1≤n≤N-1
(27)
(28)
定理2 假设问题(1)、(2)的解适当光滑。差分格式(18)~(20)式的解依L2范数收敛于问题的精确解, 收敛阶在时间方向上二阶收敛, 空间方向上四阶收敛,即存在正常数c,使得:
‖en‖≤c(τ2+h4), 0≤n≤N
(29)
证明用数学归纳法证明定理2。
由(28)式和(17)式知, 当n=0,1时,‖en‖≤c(τ2+h4)显然成立。
假设当n=0,1,2,…,m时, (29)式均成立。 设h-1(τ2+h4)≤c-1,则有:
1≤i≤M, 0≤n≤m
(30)
从而有:
1≤i≤M, 0≤n≤m
(31)
‖(Φn)3-(un)3‖=‖en[(Φn)2+Φnun+(un)2]‖≤
c1‖en‖, 0≤n≤m
(32)
(33)
根据引理2, 由(33)式可得:
(34)
由引理3, 应用柯西不等式, 结合(32)式, 由(34)式可得:
(35)
(36)
将(36)式代入(35)式, 有
(37)
在(37)式中, 将n改写成l, 并对l从1~n求和, 可得:
‖A2e1‖2-‖A2e0‖2)≤
(38)
(39)
由离散的Gronwall不等式, 可得:
1≤n≤m
(40)
从而有:
(41)
结合(41)式和引理3, 可得:
(42)
即当n=m+1时, (29)式成立。由数学归纳法知, 定理成立。
从定理2可以看出,差分格式(18)~(20)式的收敛阶为时间方向二阶收敛、空间方向四阶收敛, 提高了此前工作的收敛精度。
对于充分小的空间步长h, 定义时间收敛阶:
对于充分小的时间步长τ, 定义空间收敛阶:
对于时间和空间步长τ、h, 定义如下收敛阶:
在问题(1)、(2)中,取初值φ0(x)=0.5×sin(2πx/32),在空间区间[0,L]、时间区间[0,T]内,利用Matlab编程进行数值计算。
首先, 取参数ε=0.5, 空间步长h=0.002L,时间步长τ取T/4、T/8、T/16、T/32, 计算时刻T=1时的最大模误差和时间收敛order1, 具体结果见表1所列。
其次, 取参数ε=0.5, 时间网格点数N=5 000, 空间网格点数M取20、40、80、160,计算时刻T=1时的最大模误差和空间收敛阶order2, 具体结果见表2所列。
最后, 取参数ε=0.5, 空间步长取1/4、1/8、1/16,时间步长取1/16、1/64、1/256,计算时刻T=1时的最大模误差和收敛阶order3, 具体结果见表3所列。
表1 差分格式(18)~(20)式的最大模误差和时间收敛阶
表2 差分格式(18)~(20)式的最大模误差和空间收敛阶
表3 差分格式(18)~(20)式的最大模误差和收敛阶
从表1~表3可以看出, 差分格式在时间方向上二阶收敛、空间方向上四阶收敛,符合理论分析的结论。
本文研究了晶体相场模型的周期初边值问题的有限差分方法。首先引入中间变量,将方程降阶为耦合的二阶微分方程组,然后利用中心差商进行离散,建立空间方向四阶、时间方向二阶收敛的线性化紧差分格式,提高了数值算法的收敛阶;对差分格式进行理论分析,利用能量分析方法证明差分格式的唯一可解性及L2范数下的收敛性。数值算例表明,差分格式的收敛精度在最大模意义下达到空间四阶、时间二阶收敛。后续,将充分利用相场模型中非线性项的特点[13],进一步分析数值格式在最大模意义下的收敛性。