齐次分数阶扩散方程的加权C-N格式及其修正
Weighted C-N Scheme of HomogeneousFractional Di?usion Equations and Its Correction

作者: 陈 著 , 黄凤辉 :华南理工大学数学学院,广东 广州;

关键词: 时间分数阶扩散方程Crank-Nicolson格式收敛性Time Fractional Di?usion Equation Crank-Nicolson Scheme Convergence

摘要:
对于齐次时间分数阶扩散方程,在精确解不光滑时,数值方法精度会下降。针对这种情况,本文提出加权Crank-Nicolson格式(简记为加权C-N格式)及其修正格式,在精确解不光滑时,修正原格式的第1步后,可恢复方法的时间2阶精度。本文接着给出详细的收敛性分析,并且数值算例验证了方法的有效性。

Abstract: For the homogeneous time fractional diffusion equation, the accuracy of the numerical method will decrease when the exact solution is not smooth enough. In this case, we consider a weighted Crank-Nicolson scheme (masked as weighted C-N scheme) and its correction. After correcting the first step of the weighted C-N scheme, the second-order time accuracy can be restored. Then we give the detailed convergence analysis, and numerical examples verify the effectiveness of the scheme.

1. 引言

分数阶微分方程的数值方法的精度常依赖于精确解的光滑性。针对精确解不光滑的情况,许多学者开始研究数值方法的修正格式,以保持离散格式的高精度。例如,Lubich [1] [2] [3] 给出基于1阶和2阶向后差分格式的两种修正方法,并给出收敛性分析。Yan [4] [5] 考虑L1格式的修正格式,并提出基于分段2次插值 的新型离散格式,然后给出其修正格式。Tadjeran [6] 提出分数阶扩散方程的2阶C-N格式,Jin [7] 接着考虑分数阶C-N格式的修正。在本文中,我们考虑齐次分数阶扩散方程的加权C-N格式,并针对精确解不光滑的情况提出比较简单的修正方法,只需修正原格式的第1步,即可保持格式的2阶时间精度。本文接着给出修正格式的收敛性分析,最后通过数值算例验证方法的有效性。

本文考虑如下齐次分数阶扩散方程:

{ D 0 C u t α ( x , t ) A u ( x , t ) = 0 , L < x < R , 0 < t T u ( x , 0 ) = u 0 , L < x < R u ( L , t ) = b L ( t ) , u ( R , t ) = b R ( t ) , 0 < t T (1)

其中 0 < α 1 D 0 C u t α ( x , t ) 是Caputo时间分数阶导数,其定义为:

D 0 C u t α ( x , t ) = { 1 Γ ( 1 α ) 0 t ( t s ) α u ( s ) d s , 0 < α < 1 d u ( x , t ) d t , α = 1 (2)

算子A表示有界正则区域上的自伴正定二阶椭圆偏微分算子 [4],满足

( z A ) 1 C | z | 1 , z Σ ω = { z 0 : | arg z | < ω } , ω ( π 2 , π ) (3)

其中 · 表示 L 2 范数,记 D 0 R u t α ( x , t ) 为Riemman-Liouville导数,则有 D 0 R ( u ( x , t ) u 0 ) t α = D 0 C u t α ( x , t )

2. 加权C-N格式及其修正格式

V ( t ) = u ( t ) u 0 ,则方程(1)可表示为:

D 0 C V t α ( x , t ) A V ( x , t ) = A u 0 (4)

τ 为时间步长, t n = n τ , n = 0 , 1 , 2 , , N , 0 t n T ,h为空间步长, x i = L + i h , i = 0 , 1 , , M , L x i R ¯ t α V ( t n ) : = τ α j = 0 n b n j V j 表示 D 0 R u t α ( x , t ) 的向后Euler卷积逼近,其生成函数为:

b ˜ ( ξ ) = j = 0 b j ξ j = ( 1 ξ ) α (5)

分数阶导数的逼近格式在 t n 时刻的时间精度为1阶,在 t n α 2 τ 时刻的时间精度为2阶,对 A V ( t ) t n 时刻和 t n 1 时刻作线性拉格朗日插值,可得 t n α 2 τ 时刻的逼近格式 ( 1 α 2 ) A V n + α 2 A V n 1 。带入方程(4),可得加权C-N格式:

τ α j = 0 n b n j V j ( 1 α 2 ) A V n α 2 A V n 1 = A u 0 , n = 1 , 2 , , N (6)

在精确解不光滑的情况下,加权C-N格式达不到2阶时间精度。我们对加权C-N格式(6)的第1步得初值条件添加一个权系数 c 0 ,适当选取可使离散格式保持2阶精度,加权C-N修正格式如下:

{ τ α j = 0 n b n j V j ( 1 α 2 ) A V n α 2 A V n 1 = ( 1 + c 0 ) A u 0 , n = 1 τ α j = 0 n b n j V j ( 1 α 2 ) A V n α 2 A V n 1 = A u 0 , n = 2 , , N (7)

3. 加权C-N修正格式的收敛性分析

为了证明加权C-N修正格式的收敛性分析,我们先给出3个引理。

引理1:定义式子:

μ ( z ) = 1 e τ z ( 1 α 2 + α 2 e τ z ) 1 α c 0 e τ z + e τ z 1 e τ z 1 α 2 + α 2 e τ z (8)

本文统一令 z Γ τ Γ τ = { z Γ : | z | π τ } Σ ω [1],当权系数 c 0 = 1 α 2 时,则有:

| μ ( z ) 1 | C τ 2 | z | 2 (9)

其中C为正常数。

证明:令 w = τ z f ( w ) = μ ( z ) 1 ,则 lim x 0 f ( w ) = 0 lim x 0 f ( w ) = c 0 + α 1 2 = 0 lim x 0 f ( w ) 0 ,所以 μ ( w ) 1 = O ( w 2 ) ,即可证得 | μ ( z ) 1 | C τ 2 | z | 2

引理2:定义式子:

K ( z ) = z 1 ( z α A ) 1 A (10)

z τ = 1 e τ z ( 1 α 2 + α 2 e τ z ) 1 α (11)

则下式成立:

K ( z ) C | z | 1 (12)

K ( z τ ) C | z | 1 (13)

其中 · 表示 L 2 范数,C为正常数。

证明:令 w = τ z ,则 w 0 时, τ z τ = O ( w ) ,所以 c | τ z | | τ z τ | C | τ z | ,即 c | z | | z τ | C | z |

由引理1可得 z Σ ω , z α Σ ω ,将 z , z α 带入式(3)可得:

( z α A ) 1 C | z | α (14)

( z τ α A ) 1 C | z | α (15)

且下式成立:

( z α A ) 1 A = z α ( z α A ) 1 I (16)

带入式(10)即可证得式(12)和(13):

K ( z ) = | z | 1 z α ( z α A ) 1 I C | z | 1 (17)

K ( z τ ) = | z τ | 1 z τ α ( z τ α A ) 1 I C | z | 1 (18)

引理3: z , z α 分别由引理1和引理2定义,则下式成立:

K ( z τ ) K ( z ) C τ 2 | z | (19)

证明:令 w = τ z ,则 lim x 0 g ( w ) = lim x 0 g ( w ) = lim x 0 g ( w ) = 0 lim x 0 g ( w ) 0 ( 0 < α < 1 ) ,所以可得 g ( w ) = O ( w 3 ) ,即:

| z τ z | C τ 2 z 3 (20)

| z τ α z α | C | z | α 1 | z τ z | C τ 2 | z | 2 + α ( [8] 引理B.1),且由式(16)及可得:

K ( z τ ) K ( z ) = z τ 1 ( z τ α A ) 1 A z 1 ( z α A ) 1 A | z τ | 1 ( z τ α A ) 1 A ( z α A ) 1 A + | z τ 1 z 1 | ( z α A ) 1 A | z τ | 1 z τ α ( z τ α A ) 1 z α ( z α A ) 1 + | z τ 1 z 1 | z α ( z α A ) 1 I C | z | 1 | z τ α z α | ( z τ α A ) 1 + C | z | 1 | z | α ( z τ α A ) 1 ( z α A ) 1 + C | z z τ | | z z τ | C τ 2 | z | 1 | z | 2 + α | z | α + C | z | α 1 ( z α z τ α ) ( z τ α A ) 1 ( z α A ) 1 + C τ 2 | z | C τ 2 | z | (21)

其中 ( z τ α A ) 1 ( z α A ) 1 = ( z α z τ α ) ( z τ α A ) 1 ( z α A ) 1 ,即证得引理3。

为了给出收敛性分析,接下来我们分别借助Laplace变换和Cauchy积分公式,给出齐次分数阶扩散方程(4)的精确解和加权C-N修正格式(7)的数值解。对方程(4)作Laplace变换可得:

V ^ ( z ) = z 1 ( z α A ) 1 A u 0 (22)

对式(21)作Laplace逆变换可得方程(4)的精确解:

V ( t ) = 1 2 π i Γ e t z z 1 ( z α A ) 1 A u 0 d z = 1 2 π i Γ e t z K ( z ) u 0 d z (23)

其中 Γ = { z : | arg z | = θ } Σ ω θ ( α 2 , α 1 + α ) [4]。

考虑加权C-N修正格式(7),式子两边同乘 ξ n ,并关于n求和, n = 1 , 2 , ,令 V ˜ ( ξ ) = n = 1 V n ξ n ,于是 n = 1 ( j = 0 n b n j V j ) ξ n = ( j = 0 b j ξ j ) ( n = 1 V n ξ n ) ,由式(5)可得:

τ α b ˜ ( ξ ) V ˜ ( ξ ) ( 1 α 2 + α 2 ξ ) A V ˜ ( ξ ) = ( c 0 ξ + ξ 1 ξ ) A u 0 (24)

V ˜ ( ξ ) = c 0 ξ + ξ 1 ξ 1 α 2 + α 2 ξ [ ( 1 ξ τ ( 1 α 2 + α 2 ξ ) 1 α ) α A ] 1 A u 0 (25)

ξ = e τ z ,由Cauchy积分公式及式(8)、(10)和(11)可得加权C-N修正格式的数值解:

V n = 1 2 π i | ξ | = ρ ξ n 1 V ˜ ( ξ ) d ξ = τ 2 π i Γ τ e t n z V ˜ ( e τ z ) d z = 1 2 π i Γ τ e t n z μ ( z ) K ( z τ ) u 0 d z (26)

定理1 V ( t n ) , V n 分别为 t = t n 时刻方程(4)的精确解和加权C-N修正格式(7)的数值解,则下式成立:

V ( t n ) V n C τ 2 t n 2 u 0 (27)

证明: V ( t n ) V n = 1 2 π i Γ τ e t n z [ K ( z ) μ ( z ) K ( z τ ) ] u 0 d z + 1 2 π i Γ / Γ τ e t n z K ( z ) u 0 d z I + I I ,而 μ ( z ) K ( z τ ) K ( z ) | μ ( z ) 1 | K ( z τ ) + K ( z τ ) K ( z ) ,由引理1和引理2可得:

μ ( z ) K ( z τ ) K ( z ) C τ 2 | z | (28)

接着考虑 I , I I ,令 r = t n | z | ,c为正常数,则:

I 1 2 π Γ τ | e t n z | μ ( z ) K ( z τ ) K ( z ) u 0 d z C Γ τ | e t n z | τ 2 | z | u 0 d z C τ 2 0 e c r r t n 1 t n 1 u 0 d r C τ 2 t n 2 u 0 (29)

I 1 2 π Γ / Γ τ | e t n z | K ( z ) u 0 d z C 1 τ e c t n | z | | z | 2 | z | d z C τ 2 0 e c r r t n 1 t n 1 u 0 d r C τ 2 t n 2 u 0 (30)

则可证得 ,即权系数取 c 0 = 1 α 2 时,加权C-N修正格式为时间2阶精度。

4. 数值算例

数值算例1:考虑如下齐次分数阶扩散方程:

{ D 0 C u t α ( x , t ) 1 4 π 2 2 u t 2 = 0 , 0 < x < 1 , 0 < t T u ( x , 0 ) = sin ( 2 π x ) , 0 < x < 1 u ( 0 , t ) = u ( 1 , t ) = 0 , 0 < t T (31)

该方程的精确解为 u ( x , t ) = k = 0 ( t α ) k Γ ( α k + 1 ) sin ( 2 π x ) ,u在 t = 0 时不光滑。不同时间步长的误差如表1所示,不同 α 下的误差与对应步长的对数关系如图1所示,其中空间步长取 h = 10 4 , T = 1 ,统一取离散误差为 L 2 误差 V ( t N ) V N

Table 1. L 2 error of weighted C-N scheme (6) and weighted C-N modified scheme (7)

表1. 加权C-N格式(6)和加权C-N修正格式(7)的 L 2 误差

表1给出不同 α 及不同时间步长取值下,两种方法所得的误差。方法(a)为加权C-N格式离散,方法(b)为加权C-N修正格式离散,由表1可以看出,两种格式误差均收敛,且修正后的误差更小,方法更精确。

Figure 1. The relationship between error and time step

图1. 误差与时间步长的关系

图1给出不同 α 下误差与时间步长的关系,两坐标均为对数坐标,其中图1(a)为加权C-N格式,图1(b)为加权C-N修正格式。正如理论证明的结果一样,加权C-N格式只有1阶时间精度,而修正后的格式可达到2阶时间精度。

数值算例2:考虑如下齐次分数阶扩散方程:

{ D 0 C u t α ( x , t ) 1 4 π 2 2 u t 2 = 0 , 0 < x < 1 , 0 < t T u ( x , 0 ) = x ( 1 x ) , 0 < x < 1 u ( 0 , t ) = u ( 1 , t ) = 0 , 0 < t T (32)

该方程的精确解为 u ( x , t ) = 8 π 3 n = 1 1 ( 2 n + 1 ) 3 E α ( ( 2 n + 1 ) 2 π 2 t α ) sin ( ( 2 n + 1 ) π x ) ,u在 t = 0 时不光滑。取 α = 0. 1 ,不同时间步长下的两种格式的误差和精度如表2所示。

Table 2. Error and time accuracy of weighted C-N scheme (6) and weighted C-N modified scheme (7)

表2. 加权C-N格式(6)和加权C-N修正格式(7)的误差和时间精度

表2可以看出,随着时间步长减小,两方法的误差均减小,但加权C-N格式只有1阶时间精度,而修正格式的精度可达到2阶。

基金项目

国家自然科学基金青年科学基金项目(61802129),广东省国家青年基金纵向协同项目(2018A030310381)。

文章引用: 陈 著 , 黄凤辉 (2019) 齐次分数阶扩散方程的加权C-N格式及其修正。 应用数学进展, 8, 1611-1618. doi: 10.12677/AAM.2019.810189

参考文献

[1] Lubich, C., Sloan, I.H. and Thomee, V. (1996) Nonsmooth Data Error Estimates for Approximations of an Evolution Equation with a Positive-Type Memory Team. Mathematics of Computation, 65, 1-17.
https://doi.org/10.1090/S0025-5718-96-00677-1

[2] Lubich, C. (2004) Convolution Quadrature Revisited. Bit Numerical Mathematics, 44, 503-514.
https://doi.org/10.1023/B:BITN.0000046813.23911.2d

[3] Cuesta, E., Lubich, C. and Palencia, C. (2006) Convolution Quadrature Time Discretization of Fractional Diffusion-Wave Equations. Mathematics of Computation, 75, 673-696.
https://doi.org/10.1090/S0025-5718-06-01788-1

[4] Yan, Y.B., Khan, M. and Ford, N.J. (2018) An Analysis of the Modified L1 Scheme for Time-Fractional Partial Differential Equations with Nonsmooth Data. SIAM Journal on Numerical Analysis, 56, 210-227.
https://doi.org/10.1137/16M1094257

[5] Xing, Y.Y. and Yan, Y.B. (2018) A Higher Order Numerical Method for Time Fractional Partial Differential Equations with Nonsmooth Data. Journal of Computational Physics, 357, 305-323.
https://doi.org/10.1016/j.jcp.2017.12.035

[6] Tadjeran, C., Meerschaert, M.M. and Scheffler, H.P. (2006) A Second-Order Accurate Numerical Approximation for the Fractional Diffusion Equation. Journal of Computational Physics, 213, 205-213.
https://doi.org/10.1016/j.jcp.2005.08.008

[7] Jin, B.T., Lazarov, R. and Zhou, Z. (2018) An Analysis of the Crank-Nicolson Method for Subdiffusion. IMA Journal of Numerical Analysis, 38, 518-541.
https://doi.org/10.1093/imanum/drx019

[8] Jin, B.T., Lazarov, R. and Zhou, Z. (2017) Correction of High-Order BDF Convolution Quadrature for Fractional Evolution Equations. SIAM Journal on Scientific Computing, 39, A3129-A3152.
https://doi.org/10.1137/17M1118816

分享
Top