A fast second-order implicit difference method for time-space fractional advection-diffusion equation

TL;DR

This paper introduces a fast, second-order implicit difference method for solving time-space fractional advection-diffusion equations, combining stability, convergence, and efficient Krylov subspace solvers to reduce computational costs.

Contribution

It develops a novel second-order implicit scheme with proven stability and convergence, and designs efficient Krylov solvers with circulant preconditioners for large linear systems.

Findings

Method achieves optimal convergence order of O(τ^2 + h^2)

Krylov solvers reduce memory from O(N^2) to O(N)

Computational complexity reduced to O(N log N)

Abstract

In this paper, we consider a fast and second-order implicit difference method for approximation of a class of time-space fractional variable coefficients advection-diffusion equation. To begin with, we construct an implicit difference scheme, based on $L2-1_{\sigma}$ formula [A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, \emph{J. Comput. Phys.}, 280 (2015)] for the temporal discretization and weighted and shifted Gr\"{u}nwald method for the spatial discretization. Then, unconditional stability of the implicit difference scheme is proved, and we theoretically and numerically show that it converges in the $L_2$-norm with the optimal order $\mathcal{O}(\tau^2 + h^2)$ with time step $\tau$ and mesh size $h$. Secondly, three fast Krylov subspace solvers with suitable circulant preconditioners are designed to solve the discretized linear systems with the Toeplitz matrix. In each iterative step, these methods reduce the memory requirement of the resulting linear equations from $\mathcal{O}(N^2)$ to $\mathcal{O}(N)$ and the computational complexity from $\mathcal{O}(N^3)$ to $\mathcal{O}(N \log N)$, where $N$ is the number of grid nodes. Finally, numerical experiments are carried out to demonstrate that these methods are more practical than the traditional direct solvers of the implicit difference methods, in terms of memory requirement and computational cost.

Tables

Table 1. Table 1: L 2 subscript 𝐿 2 L_{2} -norm and maximum norm error behavior versus τ 𝜏 \tau -grid size reduction when α = 0.6 𝛼 0.6 \alpha=0.6 , β = 1.8 𝛽 1.8 \beta=1.8 , and h = 1 / 3000 ℎ 1 3000 h=1/3000 in Example 1.

$\theta$	$\tau$	${\Vert \xi \Vert }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	CO in $\parallel \cdot {\parallel }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	max${}_{0\le n\le M}{}^{}\parallel {\xi }^n{\parallel }_0$	CO in $\parallel \cdot {\parallel }_0$
0.10	1/8	5.2321e-05	–	3.3293e-05	–
	1/16	1.3096e-05	1.9982	8.3340e-06	1.9981
	1/32	3.2603e-06	2.0061	2.0748e-06	2.0061
	1/64	7.9868e-07	2.0293	5.0828e-07	2.0292
	1/128	1.8626e-07	2.1003	1.1850e-07	2.1007
0.50	1/8	3.3405e-04	–	2.1261e-04	–
	1/16	8.3929e-05	1.9928	5.3418e-05	1.9928
	1/32	2.1014e-05	1.9978	1.3374e-05	1.9979
	1/64	5.2426e-06	2.0030	3.3366e-06	2.0030
	1/128	1.2947e-06	2.0177	8.2404e-07	2.0176
0.90	1/8	6.5352e-04	–	4.1601e-04	–
	1/16	1.6365e-04	1.9976	1.0417e-04	1.9976
	1/32	4.0921e-05	1.9997	2.6048e-05	1.9997
	1/64	1.0216e-05	2.0020	6.5027e-06	2.0021
	1/128	2.5376e-06	2.0093	1.6153e-06	2.0093
0.99	1/8	7.1584e-04	–	4.5567e-04	–
	1/16	1.7899e-04	1.9998	1.1395e-04	1.9996
	1/32	4.4733e-05	2.0004	2.8479e-05	2.0005
	1/64	1.1167e-05	2.0021	7.1095e-06	2.0021
	1/128	2.7756e-06	2.0084	1.7671e-06	2.0084

Table 2. Table 2: L 2 subscript 𝐿 2 L_{2} -norm and maximum norm error behavior versus grid size reduction when α = 0.6 𝛼 0.6 \alpha=0.6 , β = 1.8 𝛽 1.8 \beta=1.8 , and τ = h 𝜏 ℎ \tau=h in Example 1.

$\theta$	$\tau$	${\Vert \xi \Vert }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	CO in $\parallel \cdot {\parallel }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	max${}_{0\le n\le M}{}^{}\parallel {\xi }^n{\parallel }_0$	CO in $\parallel \cdot {\parallel }_0$
0.10	1/20	6.0326e-04	–	4.1160e-04	–
	1/40	1.4719e-04	2.0351	9.9593e-05	2.0471
	1/80	3.5781e-05	2.0404	2.4068e-05	2.0489
	1/160	8.6941e-06	2.0411	5.8187e-06	2.0483
	1/320	2.1128e-06	2.0409	1.4087e-06	2.0464
0.50	1/20	5.5002e-04	–	3.7743e-04	–
	1/40	1.3407e-04	2.0365	9.1172e-05	2.0496
	1/80	3.2522e-05	2.0435	2.1985e-05	2.0520
	1/160	7.8829e-06	2.0446	5.3038e-06	2.0514
	1/320	1.9110e-06	2.0444	1.2815e-06	2.0492
0.90	1/20	4.8566e-04	–	3.3574e-04	–
	1/40	1.1823e-04	2.0384	8.0954e-05	2.0522
	1/80	2.8599e-05	2.0475	1.9472e-05	2.0557
	1/160	6.9105e-06	2.0491	4.6857e-06	2.0551
	1/320	1.6699e-06	2.0490	1.1296e-06	2.0524
0.99	1/20	4.7183e-04	–	3.2657e-04	–
	1/40	1.1479e-04	2.0392	7.8721e-05	2.0526
	1/80	2.7752e-05	2.0484	1.8925e-05	2.0564
	1/160	6.7015e-06	2.0500	4.5520e-06	2.0558
	1/320	1.6181e-06	2.0502	1.0969e-06	2.0530

Table 3. Table 3: L 2 subscript 𝐿 2 L_{2} -norm and maximum norm error behavior versus τ 𝜏 \tau -grid size reduction when α = 0.99 𝛼 0.99 \alpha=0.99 , β = 1.99 𝛽 1.99 \beta=1.99 , and h = 1 / 3000 ℎ 1 3000 h=1/3000 in Example 1.

$\theta$	$\tau$	${\Vert \xi \Vert }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	CO in $\parallel \cdot {\parallel }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	max${}_{0\le n\le M}{}^{}\parallel {\xi }^n{\parallel }_0$	CO in $\parallel \cdot {\parallel }_0$
0.10	1/8	5.2548e-05	–	3.3453e-05	–
	1/16	1.3149e-05	1.9987	8.3697e-06	1.9989
	1/32	3.2701e-06	2.0076	2.0799e-06	2.0086
	1/64	7.9770e-07	2.0354	5.0584e-07	2.0398
	1/128	1.8525e-07	2.1064	1.1708e-07	2.1112
0.50	1/8	3.3536e-04	–	2.1354e-04	–
	1/16	8.4257e-05	1.9929	5.3650e-05	1.9929
	1/32	2.1094e-05	1.9980	1.3430e-05	1.9982
	1/64	5.2595e-06	2.0038	3.3469e-06	2.0045
	1/128	1.2956e-06	2.0213	8.2291e-07	2.0240
0.90	1/8	6.5583e-04	–	4.1765e-04	–
	1/16	1.6423e-04	1.9976	1.0458e-04	1.9976
	1/32	4.1065e-05	1.9997	2.6150e-05	1.9998
	1/64	1.0249e-05	2.0024	6.5248e-06	2.0028
	1/128	2.5427e-06	2.0111	1.6172e-06	2.0124
0.99	1/8	7.1823e-04	–	4.5737e-04	–
	1/16	1.7958e-04	1.9998	1.1437e-04	1.9997
	1/32	4.4877e-05	2.0005	2.8580e-05	2.0006
	1/64	1.1200e-05	2.0025	7.1310e-06	2.0028
	1/128	2.7804e-06	2.0101	1.7688e-06	2.0113

Table 4. Table 4: L 2 subscript 𝐿 2 L_{2} -norm and maximum norm error behavior versus grid size reduction when α = 0.99 𝛼 0.99 \alpha=0.99 , β = 1.99 𝛽 1.99 \beta=1.99 , and τ = h 𝜏 ℎ \tau=h in Example 1.

$\theta$	$\tau$	${\Vert \xi \Vert }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	CO in $\parallel \cdot {\parallel }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	max${}_{0\le n\le M}{}^{}\parallel {\xi }^n{\parallel }_0$	CO in $\parallel \cdot {\parallel }_0$
0.10	1/20	6.0735e-04	–	4.4270e-04	–
	1/40	1.5155e-04	2.0027	1.1041e-04	2.0035
	1/80	3.7819e-05	2.0026	2.7528e-05	2.0039
	1/160	9.4336e-06	2.0032	6.8625e-06	2.0041
	1/320	2.3531e-06	2.0033	1.7107e-06	2.0041
0.50	1/20	5.5527e-04	–	4.0960e-04	–
	1/40	1.3860e-04	2.0022	1.0219e-04	2.0030
	1/80	3.4589e-05	2.0026	2.5475e-05	2.0041
	1/160	8.6263e-06	2.0035	6.3494e-06	2.0044
	1/320	2.1512e-06	2.0036	1.5825e-06	2.0045
0.90	1/20	4.9287e-04	–	3.6976e-04	–
	1/40	1.2317e-04	2.0006	9.2308e-05	2.0021
	1/80	3.0734e-05	2.0027	2.3012e-05	2.0041
	1/160	7.6637e-06	2.0037	5.7346e-06	2.0046
	1/320	1.9108e-06	2.0039	1.4289e-06	2.0048
0.99	1/20	4.7956e-04	–	3.6117e-04	–
	1/40	1.1989e-04	2.0000	9.0184e-05	2.0017
	1/80	2.9914e-05	2.0028	2.2483e-05	2.0041
	1/160	7.4589e-06	2.0038	5.6023e-06	2.0047
	1/320	1.8596e-06	2.0040	1.3958e-06	2.0049

Table 5. Table 5: L 2 subscript 𝐿 2 L_{2} -norm and maximum norm error behavior versus τ 𝜏 \tau -grid size reduction when α = 0.6 𝛼 0.6 \alpha=0.6 , β = 1.8 𝛽 1.8 \beta=1.8 and h = 1 / 800 ℎ 1 800 h=1/800 in Example 2.

$\theta$	$\tau$	${\Vert \xi \Vert }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	CO in $\parallel \cdot {\parallel }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	max${}_{0\le n\le M}{}^{}\parallel {\xi }^n{\parallel }_0$	CO in $\parallel \cdot {\parallel }_0$
0.10	1/8	6.7718e-04	–	4.3108e-04	–
	1/16	1.7300e-04	1.9688	1.1009e-04	1.9693
	1/32	4.1798e-05	2.0493	2.6564e-05	2.0511
	1/64	8.3378e-06	2.3257	5.2789e-06	2.3312
	1/128	1.6425e-06	2.3438	4.9654e-07	3.4103
0.50	1/8	2.8747e-03	–	1.8297e-03	–
	1/16	7.3181e-04	1.9739	4.6572e-04	1.9741
	1/32	1.8254e-04	2.0033	1.1613e-04	2.0037
	1/64	4.3648e-05	2.0642	2.7733e-05	2.0661
	1/128	8.7404e-06	2.3201	5.5326e-06	2.3256
0.90	1/8	4.0951e-03	–	2.6098e-03	–
	1/16	1.0327e-03	1.9875	6.5724e-04	1.9894
	1/32	2.5667e-04	2.0084	1.6330e-04	2.0089
	1/64	6.2038e-05	2.0487	3.9436e-05	2.0499
	1/128	1.3330e-05	2.2185	8.4480e-06	2.2228
0.99	1/8	4.1885e-03	–	2.6801e-03	–
	1/16	1.0594e-03	1.9832	6.7470e-04	1.9900
	1/32	2.6328e-04	2.0086	1.6752e-04	2.0099
	1/64	6.3658e-05	2.0482	4.0472e-05	2.0493
	1/128	1.3737e-05	2.2123	8.7084e-06	2.2164

Table 6. Table 6: L 2 subscript 𝐿 2 L_{2} -norm and maximum norm error behavior versus grid size reduction when α = 0.6 𝛼 0.6 \alpha=0.6 , β = 1.8 𝛽 1.8 \beta=1.8 and τ = h 𝜏 ℎ \tau=h in Example 2.

$\theta$	$\tau$	${\Vert \xi \Vert }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	CO in $\parallel \cdot {\parallel }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	max${}_{0\le n\le M}{}^{}\parallel {\xi }^n{\parallel }_0$	CO in $\parallel \cdot {\parallel }_0$
0.10	1/20	5.2772e-03	–	3.6233e-03	–
	1/40	1.2913e-03	2.0309	8.7967e-04	2.0423
	1/80	3.1470e-04	2.0368	2.1331e-04	2.0440
	1/160	7.6674e-05	2.0372	5.1742e-05	2.0435
	1/320	1.8683e-05	2.0370	1.2566e-05	2.0418
0.50	1/20	4.8930e-03	–	3.3744e-03	–
	1/40	1.1943e-03	2.0346	8.1793e-04	2.0446
	1/80	2.9056e-04	2.0393	1.9797e-04	2.0467
	1/160	7.0647e-05	2.0401	4.7934e-05	2.0462
	1/320	1.7181e-05	2.0398	1.1621e-05	2.0443
0.90	1/20	4.6652e-03	–	3.2255e-03	–
	1/40	1.1378e-03	2.0357	7.8171e-04	2.0448
	1/80	2.7670e-04	2.0399	1.8908e-04	2.0476
	1/160	6.7216e-05	2.0414	4.5748e-05	2.0472
	1/320	1.6330e-05	2.0413	1.1083e-05	2.0454
0.99	1/20	4.6390e-03	–	3.2055e-03	–
	1/40	1.1310e-03	2.0362	7.7718e-04	2.0442
	1/80	2.7508e-04	2.0397	1.8801e-04	2.0474
	1/160	6.6816e-05	2.0416	4.5488e-05	2.0473
	1/320	1.6231e-05	2.0414	1.1020e-05	2.0454

Table 7. Table 7: L 2 subscript 𝐿 2 L_{2} -norm and maximum norm error behavior versus τ 𝜏 \tau -grid size reduction when α = 0.99 𝛼 0.99 \alpha=0.99 , β = 1.99 𝛽 1.99 \beta=1.99 and h = 1 / 800 ℎ 1 800 h=1/800 in Example 2.

$\theta$	$\tau$	${\Vert \xi \Vert }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	CO in $\parallel \cdot {\parallel }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	max${}_{0\le n\le M}{}^{}\parallel {\xi }^n{\parallel }_0$	CO in $\parallel \cdot {\parallel }_0$
0.10	1/8	6.7919e-04	–	4.3233e-04	–
	1/16	1.7333e-04	1.9703	1.1015e-04	1.9727
	1/32	4.1678e-05	2.0562	2.6300e-05	2.0663
	1/64	8.0974e-06	2.3638	4.9278e-06	2.4160
	1/128	9.5365e-07	3.0859	6.4747e-07	2.9281
0.50	1/8	2.8857e-03	–	1.8374e-03	–
	1/16	7.3455e-04	1.9740	4.6750e-04	1.9746
	1/32	1.8305e-04	2.0046	1.1631e-04	2.0070
	1/64	4.3571e-05	2.0708	2.7500e-05	2.0805
	1/128	8.5104e-06	2.3561	5.1893e-06	2.4058
0.90	1/8	4.1053e-03	–	2.6172e-03	–
	1/16	1.0367e-03	1.9855	6.5998e-04	1.9875
	1/32	2.5752e-04	2.0092	1.6373e-04	2.0111
	1/64	6.2046e-05	2.0533	3.9263e-05	2.0601
	1/128	1.3122e-05	2.2414	8.1204e-06	2.2735
0.99	1/8	4.1843e-03	–	2.6789e-03	–
	1/16	1.0619e-03	1.9783	6.7724e-04	1.9839
	1/32	2.6411e-04	2.0074	1.6791e-04	2.0120
	1/64	6.3651e-05	2.0529	4.0288e-05	2.0593
	1/128	1.3524e-05	2.2347	8.3771e-06	2.2658

Table 8. Table 8: L 2 subscript 𝐿 2 L_{2} -norm and maximum norm error behavior versus grid size reduction when α = 0.99 𝛼 0.99 \alpha=0.99 , β = 1.99 𝛽 1.99 \beta=1.99 and τ = h 𝜏 ℎ \tau=h in Example 2.

$\theta$	$\tau$	${\Vert \xi \Vert }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	CO in $\parallel \cdot {\parallel }_{𝒞\left({\overline{\omega }}_{h\tau }\right)}$	max${}_{0\le n\le M}{}^{}\parallel {\xi }^n{\parallel }_0$	CO in $\parallel \cdot {\parallel }_0$
0.10	1/20	5.1023e-03	–	3.7233e-03	–
	1/40	1.2730e-03	2.0029	9.2852e-04	2.0036
	1/80	3.1774e-04	2.0023	2.3151e-04	2.0039
	1/160	7.9261e-05	2.0032	5.7718e-05	2.0040
	1/320	1.9772e-05	2.0032	1.4389e-05	2.0041
0.50	1/20	4.7154e-03	–	3.4794e-03	–
	1/40	1.1767e-03	2.0026	8.6769e-04	2.0036
	1/80	2.9361e-04	2.0028	2.1629e-04	2.0042
	1/160	7.3226e-05	2.0035	5.3912e-05	2.0043
	1/320	1.8263e-05	2.0034	1.3438e-05	2.0043
0.90	1/20	4.4905e-03	–	3.3363e-03	–
	1/40	1.1217e-03	2.0012	8.3257e-04	2.0026
	1/80	2.7992e-04	2.0026	2.0758e-04	2.0039
	1/160	6.9818e-05	2.0033	5.1743e-05	2.0042
	1/320	1.7413e-05	2.0034	1.2897e-05	2.0043
0.99	1/20	4.4668e-03	–	3.3193e-03	–
	1/40	1.1155e-03	2.0016	8.2848e-04	2.0023
	1/80	2.7838e-04	2.0026	2.0657e-04	2.0038
	1/160	6.9435e-05	2.0033	5.1492e-05	2.0042
	1/320	1.7317e-05	2.0035	1.2834e-05	2.0044

Table 9. Table 9: CPU time in seconds for solving Example 1 with θ = 0.9 𝜃 0.9 \theta=0.9 , α = 0.8 𝛼 0.8 \alpha=0.8 , β = 1.9 𝛽 1.9 \beta=1.9 .

		$\mathrm{PBiCGSTAB}$		$\mathrm{PGPBiCOR}(3,1)$
$h=\tau$	Time1	Time2(Iter)	Speed-up	Time3(Iter)	Speed-up
$2^{-6}$	0.031	0.161(6.9)	0.19	0.213(5.6)	0.14
$2^{-7}$	0.145	0.488(7.0)	0.30	0.612(6.1)	0.24
$2^{-8}$	1.624	1.217(7.0)	1.33	1.500(6.3)	1.08
$2^{-9}$	17.375	4.828(7.0)	3.60	5.592(6.5)	3.11
$2^{-10}$	143.968	17.544(7.0)	8.21	19.393(6.7)	7.42

Table 10. Table 10: CPU time in seconds for solving Example 1 with θ = 0.7 𝜃 0.7 \theta=0.7 , α = 0.6 𝛼 0.6 \alpha=0.6 , β = 1.5 𝛽 1.5 \beta=1.5 .

		$\mathrm{PBiCGSTAB}$		$\mathrm{PGPBiCOR}(3,1)$
$h=\tau$	Time1	Time2(Iter)	Speed-up	Time3(Iter)	Speed-up
$2^{-6}$	0.031	0.167(6.8)	0.18	0.154(6.0)	0.20
$2^{-7}$	0.149	0.509(7.0)	0.29	0.485(6.3)	0.31
$2^{-8}$	1.701	1.335(7.8)	1.27	1.281(7.2)	1.33
$2^{-9}$	17.251	5.325(8.0)	3.24	5.275(7.6)	3.27
$2^{-10}$	143.739	19.166(8.1)	7.50	18.261(7.8)	7.87

Table 11. Table 11: CPU time in seconds for solving Example 2 with θ = 0.9 𝜃 0.9 \theta=0.9 , α = 0.8 𝛼 0.8 \alpha=0.8 , β = 1.9 𝛽 1.9 \beta=1.9 .

		$\mathrm{PBiCGSTAB}$		$\mathrm{PGPBiCOR}(3,1)$
$h=\tau$	Time1(Iter1)	Time2(Iter1, Iter2)	Speed-up	Time3(Iter1,Iter2)	Speed-up
$2^{-6}$	0.038(5.0)	0.470(5.0, 7.0)	0.08	0.482(5.0, 6.0)	0.08
$2^{-7}$	0.164(4.1)	1.414(4.1, 7.0)	0.12	1.420(4.1, 6.0)	0.12
$2^{-8}$	1.922(4.0)	3.325(4.0, 7.0)	0.58	3.350(4.0, 6.0)	0.57
$2^{-9}$	19.122(4.0)	13.923(4.0, 7.0)	1.37	13.879(4.0, 6.0)	1.38
$2^{-10}$	159.561(3.5)	42.497(3.5, 7.4)	3.75	41.005(3.5, 6.0)	3.89

Table 12. Table 12: CPU time in seconds for solving Example 2 with θ = 0.7 𝜃 0.7 \theta=0.7 , α = 0.6 𝛼 0.6 \alpha=0.6 , β = 1.7 𝛽 1.7 \beta=1.7 .

		$\mathrm{PBiCGSTAB}$		$\mathrm{PGPBiCOR}(3,1)$
$h=\tau$	Time1(Iter1)	Time2(Iter1, Iter2)	Speed-up	Time3(Iter1, Iter2)	Speed-up
$2^{-6}$	0.037(5.0)	0.480(5.0, 7.0)	0.08	0.476(5.0, 6.0)	0.08
$2^{-7}$	0.171(5.0)	1.666(5.0, 7.0)	0.10	1.866(5.0, 7.0)	0.09
$2^{-8}$	1.802(4.3)	3.975(4.3, 8.1)	0.45	4.022(4.3, 7.0)	0.45
$2^{-9}$	18.950(4.0)	15.604(4.0, 8.0)	1.21	15.873(4.0, 7.2)	1.19
$2^{-10}$	164.387(4.0)	50.313(4.0, 8.0)	3.27	54.279(4.0, 7.8)	3.03

Equations

$$\{\begin{array}{l}\begin{array}{rl}{D}_{0,t}^{\theta }u(x,t)=& -d_{+}(t)D_{0,x}^{\alpha }u(x,t)-d_{-}(t)D_{x,L}^{\alpha }u(x,t)\\ & +e_{+}(t)D_{0,x}^{\beta }u(x,t)+e_{-}(t)D_{x,L}^{\beta }u(x,t)+f(x,t),\end{array}\\ u(x,0)=u_0(x),0\le x\le L,\\ u(0,t)=u(L,t)=0,0\le t\le T,\end{array}$$

$$D_{0,t}^{\theta }u(x,t)=\frac{1}{\Gamma (1-\theta )}{\int }_0^t\frac{\partial u(x,\xi )}{\partial \xi }\frac{d\xi }{(t-\xi {)}^{\theta }},$$

$$D_{0,x}^{\gamma }u(x,t)=\frac{1}{\Gamma (n-\gamma )}\frac{d^n}{d{x}^n}{\int }_0^x\frac{u(\eta ,t)}{(x-\eta {)}^{\gamma -n+1}}d\eta ,$$

$$D_{x,L}^{\gamma }u(x,t)=\frac{(-1{)}^n}{\Gamma (n-\gamma )}\frac{d^n}{d{x}^n}{\int }_x^L\frac{u(\eta ,t)}{(\eta -x{)}^{\gamma -n+1}}d\eta ,$$

$$\left|D_{0,t}^{\theta }y(t)-{\Delta }_{0,t_{j+\sigma }}^{\theta }y(t)\right|=\mathcal{O}({\tau }^{3-\theta }),$$

$$\Delta^{\theta}_{0,t_{j+\sigma}}y(t)=\frac{\tau^{-\theta}}{\Gamma(2-\theta)}\sum\limits^{j}_{s=0}c^{(\theta,\sigma)}_{j-s}\big{[}y(t_{s+1})-y(t_{s})\big{]},$$

$$c_s^{(\theta ,\sigma )}=\{\begin{array}{ll}{a}_0^{(\theta ,\sigma )}+b_1^{(\theta ,\sigma )},& s=0,\\ a_s^{(\theta ,\sigma )}+b_{s+1}^{(\theta ,\sigma )}-b_s^{(\theta ,\sigma )},& 1\le s\le j-1,\\ a_j^{(\theta ,\sigma )}-b_j^{(\theta ,\sigma )},& s=j,\end{array}$$

$$D_{0,x}^{\gamma }u(x,t)=\frac{1}{h^{\gamma }}\sum _{k=0}^{\left[\frac{x}{h}\right]+1}{\omega }_k^{(\gamma )}u(x-(k-1)h,t)+\mathcal{O}(h^2),$$

$$D_{x,L}^{\gamma }u(x,t)=\frac{1}{h^{\gamma }}\sum _{k=0}^{\left[\frac{L-x}{h}\right]+1}{\omega }_k^{(\gamma )}u(x+(k-1)h,t)+\mathcal{O}(h^2).$$

$${\delta }_{x,+}^{\gamma }{u}_i^j=\frac{1}{h^{\gamma }}\sum _{k=0}^{i+1}{\omega }_k^{(\gamma )}{u}_{i-k+1}^j,{\delta }_{x,-}^{\gamma }{u}_i^j=\frac{1}{h^{\gamma }}\sum _{k=0}^{N-i+1}{\omega }_k^{(\gamma )}{u}_{i+k-1}^j,$$

$${\omega }_0^{(\gamma )}=\frac{\gamma }{2}{g}_0^{(\gamma )},{\omega }_k^{(\gamma )}=\frac{\gamma }{2}{g}_k^{(\gamma )}+\frac{2-\gamma }{2}{g}_{k-1}^{(\gamma )},k\ge 1,$$

$$g_{0}^{(\gamma)}=1,\qquad g_{k}^{(\gamma)}=\Big{(}1-\frac{\gamma+1}{k}\Big{)}g_{k-1}^{(\gamma)},~{}k=1,2,\ldots$$

$$u_i^{(j+\sigma )}=\sigma u_i^{j+1}+(1-\sigma )u_i^j,d_{\pm }^{j+\sigma }=d_{\pm }(t_{j+\sigma }),e_{\pm }^{j+\sigma }=e_{\pm }(t_{j+\sigma }),f_i^{j+\sigma }=f(x_i,t_{j+\sigma }),$$

$${\delta }_h^{\alpha ,\beta }{u}_i^{(j+\sigma )}=-d_{+}^{j+\sigma }{\delta }_{x,+}^{\alpha }{u}_i^{(j+\sigma )}-d_{-}^{j+\sigma }{\delta }_{x,-}^{\alpha }{u}_i^{(j+\sigma )}+e_{+}^{j+\sigma }{\delta }_{x,+}^{\beta }{u}_i^{(j+\sigma )}+e_{-}^{j+\sigma }{\delta }_{x,-}^{\beta }{u}_i^{(j+\sigma )}.$$

$${\Delta }_{0,t_{j+\sigma }}^{\theta }{u}_i={\delta }_h^{\alpha ,\beta }{u}_i^{(j+\sigma )}+f_i^{j+\sigma }.$$

$$\{\begin{array}{ll}{\Delta }_{0,t_{j+\sigma }}^{\theta }{u}_i={\delta }_h^{\alpha ,\beta }{u}_i^{(j+\sigma )}+f_i^{j+\sigma },& 1\le i\le N-1,0\le j\le M-1,\\ u_i^0=u_0(x_i),& 1\le i\le N-1,\\ u_0^j=u_N^j=0,& 0\le j\le M-1.\end{array}$$

$$\{\begin{array}{l}{\omega }_0^{(\alpha )}=\frac{\alpha }{2}>0,{\omega }_1^{(\alpha )}=\frac{2-\alpha -{\alpha }^2}{2}>0,{\omega }_2^{(\alpha )}=\frac{\alpha ({\alpha }^2+\alpha -4)}{4}<0,\\ {\omega }_2^{(\alpha )}<{\omega }_3^{(\alpha )}<{\omega }_4^{(\alpha )}<\cdots <0,{\omega }_0^{(\alpha )}+{\omega }_2^{(\alpha )}<0,\\ \sum _{k=0}^{\infty }{\omega }_k^{(\alpha )}=0,\sum _{k=0}^N{\omega }_k^{(\alpha )}>0,N\ge 1.\end{array}$$

$$\{\begin{array}{l}{\omega }_0^{(\beta )}=\frac{\beta }{2}>0,{\omega }_1^{(\beta )}=\frac{2-\beta -{\beta }^2}{2}<0,{\omega }_2^{(\beta )}=\frac{\beta ({\beta }^2+\beta -4)}{4},\\ 1\ge {\omega }_0^{(\beta )}\ge {\omega }_3^{(\beta )}\ge {\omega }_4^{(\beta )}\ge \cdots \ge 0,{\omega }_0^{(\beta )}+{\omega }_2^{(\beta )}>0,\\ \sum _{k=0}^{\infty }{\omega }_k^{(\beta )}=0,\sum _{k=0}^N{\omega }_k^{(\beta )}<0,N\ge 2.\end{array}$$

$$(\mathbf{u},\mathbf{v})=h\sum _{i=1}^{N-1}{u}_i{v}_i,\Vert \mathbf{u}\Vert =\sqrt{(\mathbf{u},\mathbf{u})},for\forall \mathbf{u},\mathbf{v}\in V_h,$$

$$V_h=\{\mathbf{v}\mid \mathbf{v}=\{v_i\}isagridfunctionon{\bar{\omega }}_{h}and{v}_i=0ifi=0,N\}.$$

$$(\delta^{\alpha}_{x,+}\bm{v},\bm{v})=(\delta^{\alpha}_{x,-}\bm{v},\bm{v})\geq\Big{(}\frac{1}{h^{\alpha}}\sum\limits_{k=0}^{N-1}\omega_{k}^{(\alpha)}\Big{)}\left\|\bm{v}\right\|^{2}.$$

$$({\delta }_{x,+}^{\alpha }\mathbf{v},\mathbf{v})=({\delta }_{x,-}^{\alpha }\mathbf{v},\mathbf{v})$$

$$\begin{split}(\delta^{\alpha}_{x,+}\bm{v},\bm{v})&=h^{1-\alpha}\sum\limits_{i=1}^{N-1}\Big{(}\sum\limits_{k=0}^{i+1}\omega_{k}^{(\alpha)}v_{i-k+1}\Big{)}v_{i}\\ &=\omega_{1}^{(\alpha)}h^{1-\alpha}\sum\limits_{i=1}^{N-1}(v_{i})^{2}+(\omega_{0}^{(\alpha)}+\omega_{2}^{(\alpha)})h^{1-\alpha}\sum\limits_{i=1}^{N-2}v_{i}v_{i+1}+\sum\limits_{k=3}^{N-1}\omega_{k}^{(\alpha)}h^{1-\alpha}\sum\limits_{i=1}^{N-k}v_{i+k-1}v_{i}\\ &\geq\omega_{1}^{(\alpha)}h^{-\alpha}\left\|\bm{v}\right\|^{2}+(\omega_{0}^{(\alpha)}+\omega_{2}^{(\alpha)})h^{1-\alpha}\sum\limits_{i=1}^{N-2}\frac{(v_{i})^{2}+(v_{i+1})^{2}}{2}\\ &~{}~{}~{}+\sum\limits_{k=3}^{N-1}\omega_{k}^{(\alpha)}h^{1-\alpha}\sum\limits_{i=1}^{N-k}\frac{(v_{i+k-1})^{2}+(v_{i})^{2}}{2}\geq\Big{(}\frac{1}{h^{\alpha}}\sum\limits_{k=0}^{N-1}\omega_{k}^{(\alpha)}\Big{)}\left\|\bm{v}\right\|^{2}.\end{split}$$

$$({\delta }_{x,+}^{\alpha }\mathbf{v},\mathbf{v})=({\delta }_{x,-}^{\alpha }\mathbf{v},\mathbf{v})>c_1\ln 2{\Vert \mathbf{v}\Vert }^2.$$

$$\sum _{k=N}^{2N+2}{\omega }_k^{(\alpha )}=\sum _{k=N}^{2N}{g}_k^{(\alpha )}+\frac{2-\alpha }{2}{g}_{N-1}^{(\alpha )}+\frac{\alpha }{2}{g}_{2N+2}^{(\alpha )}+g_{2N+1}^{(\alpha )}.$$

$$\sum _{k=N}^{2N+2}{\omega }_k^{(\alpha )}<\sum _{k=N}^{2N}{g}_k^{(\alpha )},N\ge 2.$$

$$\begin{split}\frac{1}{h^{\alpha}}\sum\limits_{k=N}^{\infty}\Big{(}-\omega_{k}^{(\alpha)}\Big{)}&>\frac{1}{h^{\alpha}}\sum\limits_{k=N}^{2N+2}\Big{(}-\omega_{k}^{(\alpha)}\Big{)}>\frac{1}{h^{\alpha}}\sum\limits_{k=N}^{2N}\Big{(}-g_{k}^{(\alpha)}\Big{)}\\ &\geq\tilde{c}_{1}\sum\limits_{k=N}^{2N}k^{-(\alpha+1)}N^{\alpha}>\tilde{c}_{1}\sum\limits_{k=N}^{2N}k^{-(\alpha+1)}\Big{(}\frac{k}{2}\Big{)}^{\alpha}\\ &=\frac{\tilde{c}_{1}}{2^{\alpha}}\sum\limits_{k=N}^{2N}\frac{1}{k}>c_{1}\int\limits_{N}^{2N+1}\frac{1}{x}dx\geq c_{1}\int\limits_{N}^{2N}\frac{1}{x}dx=c_{1}\ln 2,\quad N\geq 2.\end{split}$$

$$(\delta^{\alpha}_{x,+}\bm{v},\bm{v})=(\delta^{\alpha}_{x,-}\bm{v},\bm{v})\geq\frac{1}{h^{\alpha}}\sum\limits_{k=N}^{\infty}\Big{(}-\omega_{k}^{(\alpha)}\Big{)}\left\|\bm{v}\right\|^{2}>c_{1}\ln 2\left\|\bm{v}\right\|^{2}.$$

$$(\delta^{\beta}_{x,+}\bm{v},\bm{v})=(\delta^{\beta}_{x,-}\bm{v},\bm{v})\leq\Big{(}\frac{1}{h^{\beta}}\sum\limits_{k=0}^{N-1}\omega_{k}^{(\beta)}\Big{)}\left\|\bm{v}\right\|^{2}.$$

$$(-{\delta }_{x,+}^{\beta }\mathbf{v},\mathbf{v})=(-{\delta }_{x,-}^{\beta }\mathbf{v},\mathbf{v})>c_2\ln 2{\Vert \mathbf{v}\Vert }^2.$$

Full text

Key words: Linear/Nonlinear fractional advection-diffusion equation; Weighted and shifted Grünwald scheme; $L2$-$1_{\sigma}$ formula; Krylov subspace method; Toeplitz matrix; Fast Fourier transform; Circulant preconditioner.

1 Introduction

This manuscropt is concerned with a fast second-order implicit difference method (IDM) for solving the initial-boundary value problem of the time-space fractional advection-diffusion equation (TSFADE) [1, 2]:

[TABLE]

where $\theta,\alpha\in(0,1]$, $\beta\in(1,2]$, $0<x<L$, and $0<t\leq T$. Here, the parameters $\theta,\alpha$ and $\beta$ are the order of the TSFADE, $f(x,t)$ is the source term, and diffusion coefficient functions $d_{\pm}(t)$ and $e_{\pm}(t)$ are non-negative under the assumption that the flow is from left to right. The TSFADE (1.1) can be regarded as generalizations of classical advection-diffusion equations with the first-order time derivative replaced by the Caputo fractional derivative of order $\theta\in(0,1]$, and the first-order and second-order space derivatives replaced by the two-sided Riemann-Liouville fractional derivatives of order $\alpha\in(0,1]$ and $\beta\in(1,2]$. Namely, the time fractional derivative in (1.1) is the Caputo fractional derivative[3, 4] of order $\theta$ denoted by

[TABLE]

and the left ($D_{0,x}^{\alpha}$, $D_{0,x}^{\beta}$) and right ($D_{x,L}^{\alpha}$, $D_{x,L}^{\beta}$) space fractional derivatives in (1.1) are the Riemann-Liouville fractional derivatives[3, 4] defined as

(1) left Riemann-Liouville fractional derivative:

[TABLE]

(2) right Riemann-Liouville fractional derivative:

[TABLE]

where $\Gamma(\cdot)$ denotes the Gamma function and $n-1\leq\gamma<n$ ($n$ is a positive integer). In reality, when $\theta=\alpha=1$ and $\beta=2$, the above equation (1.1) reduces to the classical advection-diffusion equation (ADE).

There have been many studies and applications of the fractional ADE. It should be point out that the fractional ADE provides more adequate and accurate description of the movement of solute in an aquifer than traditional second-order ADEs do [5, 6], and also used to approach the description of transport dynamics in complex systems governed by anomalous diffusion and non-exponential relaxation patterns[7]. Moreover, the fractional ADE is a more suitable model for many problems, such as engineering, chemistry, entropy[8] and hydrology[9]. To obtain the analytical solutions of fractional partial differential equations [3], numerous analytical methods, such as the Fourier transform method, Adomian s decomposition method, the Laplace transform method, shifted Legendre polynomials[10] and the Mellin transform method, have been developed. Since there are very few cases in which the closed-form analytical solutions are available, or the obtained analytical solutions are less practical (expressed by the transcendental functions or infinite series). Researches on numerical approximation and techniques for the solution of fractional difference equations (FDEs) have attracted intensive interest; see [11, 12, 13, 14, 15, 16, 17, 18, 19] and references therein.

Regarding numerical methods for fractional advection-diffusion problems. In the last years, most early established numerical methods are developed for handling the space or time fractional ADE. For the time fractional ADE, many early implicit numerical methods are derived by the combination of the $L1$ approximate formula[20, 21]. These methods are unconditionally stable, but the time accuracy can not meet second order. In addition, some other related numerical methods have already been proposed for handling the time fractional ADE; see [25, 26, 27] for details. On the other hand, for the space fractional ADE, many researchers exploited the conventional shifted Grünwald discretization[22] and the implicit Euler (or Crank-Nicolson) time-step discretization for two-sided Riemann-Liouville fractional derivatives and the first order time derivative, respectively. Later, Sousa and Li[23] derived a weighted finite difference scheme for producing the novel numerical methods, which achieved the second order accuracy in both time and space for space fractional ADE. Qu and Lei et al. employed circulant and skew-circulant splitting iteration, which could reduce the required algorithm storage, for implementing the above mentioned second order numerical method; see [24] for details.

Contrarily, although the numerical methods for space or time fractional ADE were extensively investigated in the past researches, the work about numerically handling the TSFADE is not too much. Firstly, Liu et al. [22] and Zhang [28, 29] worked out a series of studies about constructing the implicit difference scheme (IDS) for TSFADE, however all these numerical schemes achieve the convergence with first order accuracy in both space and time from both the theoretical and numerical perspectives. Liu et al. [1] also considered a space-time fractional advection dispersion equation and the implicit proposed difference method has convergence $\mathcal{O}(\tau+h)$. Then, Shen et al.[30] presented an explicit difference approximation with convergence $\mathcal{O}(\tau+h)$ and an implicit difference approximation with convergence rate $\mathcal{O}(\tau^{2-\alpha}+h)$. Inspired by Liu’s paper [1], Zhao et al.[2] introduced preconditioned iterative methods to reduce the amount of computation while solving TSFADE, but its convergence can only reach $\mathcal{O}(\tau^{2-\alpha}+h)$. Later, Gu et al. [50] established fast solution techniques to solve the discretization linear system of the time-space fractional convection-diffusion equation, which has convergence rate $\mathcal{O}(\tau^{2}+h^{2})$. Furthermore, most of these methods have no complete theoretical analysis for both stability and convergence, refer to [31, 32, 33, 34] for details.

Since the fractional differential operator is nonlocal, it shows that a naive discretization of the FDE, even though implicit, leads to unconditionally unstable [14, 15]. Moreover, traditional methods for solving FDEs tend to generate full coefficient matrices, which require computational cost of $\mathcal{O}(N^{3})$ and storage of $\mathcal{O}(N^{2})$ [35]. There is no doubt that if $N$ is not small, both calculation time and storage requirements are tremendous, and this situation is that we try to avoid. To optimize the computational complexity, Meerschaet and Tadjeran [14, 15] proposed a shifted Grünwald discretization to approximate FDEs, which has been proved to be unconditionally stable. Later, Wang, Wang and Sircar [36] discovered that the full coefficient matrix derived by the shifted Grünwald scheme holds a Toeplitz-like structure. More precisely, this matrix can be written as a sum of diagonal-multiply-Toeplitz matrices. This means that the matrix-vector multiplication for Toeplitz matrix can be computed by the fast Fourier transform (FFT) with $\mathcal{O}(N\log N)$ operations [37], and the storage requirement reduces to $\mathcal{O}(N)$. Thanks to the Toeplitz-like structure, Wang and Wang [38] employed the CGNR method to solve discretized linear system of FDE by the Meerschaet-Tadjeran’s method, and the cost of each iteration is $\mathcal{O}(N\log N)$. Numerical experiments indicate that the CGNR method is fast when the diffusion coefficients are small enough, namely, the discretized system is well-conditioned. Nevertheless, the discretized system will become ill-conditioned when the diffusion coefficients are not small, hence the CGNR method converges very slowly. To overcome this shortcoming, Zhao et al. extended the precondition technique, for handling the Toeplitz-like matrix with structure as the sum of diagonal-multiply-Toeplitz matrices [2]. Their results are related to the promising acceleration of the convergence of the iterative methods, while solving (1.1), see [39, 40] for other preconditioning methods. For the same reason, we propose two preconditioned iterative methods, i.e., the preconditioned biconjugate gradient stabilized (BiCGSTAB) method [41] and the preconditioned generalized product-type solvers based on BiCOR (GPBiCOR($m$, $\ell$)) method [42], and observe results related to the acceleration of the convergence of the iterative methods to solve (1.1).

The rest of this paper is arranged as follows. In Section 2, the implicit difference scheme for (1.1) is presented, and we prove that this scheme is unconditionally stable and convergent with the accuracy of $\mathcal{O}(\tau^{2}+h^{2})$. Then, we extend the problem to nonlinear situation through utilizing the same technique. In Section 3, the BiCGSTAB method and the GPBiCOR($m$, $\ell$) method with suitable circulant preconditioners are proposed to solve the discretized system. In Section 4, numerical results are reported to demonstrate the efficiency of our numerical approaches, and concluding remarks are given in Section 5.

2 An implicit difference scheme for TSFADE

In this section, we present an implicit difference method for discretizing the TSFADE defined by (1.1). Unlike the former numerical approaches [1, 2, 30], we exploit henceforth two-sided fractional derivatives to approximate the Riemann-Liouville derivatives in (1.3) and (1.4). we can show that, by two-sided fractional derivatives, this proposed method is also unconditionally stable and convergent second order accuracy in both time and space.

2.1 Discretization of the TSFADE

In order to derive the proposed scheme, we first introduce the mesh $\bar{\omega}_{h\tau}=\bar{\omega}_{h}\times\bar{\omega}_{\tau}$, where $\bar{\omega}_{h}=\{x_{i}=ih,i=0,1,\cdots,N;x_{0}=0,x_{N}=L\}$ and $\bar{\omega}_{\tau}=\{t_{j}=j\tau,j=0,1,\cdots,M;t_{M}=T\}$. Besides, $\bm{v}=\{v_{i}\mid 0\leq i\leq N\}$ be any grid function. Then the following lemma introduced in [43] gives a description on the Caputo fractional derivative discretization.

Lemma 2.1**.**

suppose $0<\theta<1$, $\sigma=1-\frac{\theta}{2}$, $y(t)\in\mathcal{C}^{3}[0,T]$, and $t_{j+\sigma}=(j+\sigma)\tau$. Then

[TABLE]

where

[TABLE]

and for $j=0$, $c^{(\theta,\sigma)}_{0}=a^{(\theta,\sigma)}_{0}$, for $j\geq 1$,

[TABLE]

in which $a^{(\theta,\sigma)}_{0}=\sigma^{1-\theta},~{}a^{(\theta,\sigma)}_{l}=(l+\sigma)^{1-\theta}-(l-1+\sigma)^{1-\theta},~{}for~{}l\geq 1$; and b^{(\theta,\sigma)}_{l}=\frac{1}{2-\theta}\Big{[}(l+\sigma)^{2-\theta}-(l-1+\sigma)^{2-\theta}\Big{]}-\frac{1}{2}\Big{[}(l+\sigma)^{1-\theta}-(l-1+\sigma)^{1-\theta}\Big{]}.

On the other hand, the Riemann-Liouville derivatives (1.3) and (1.4) can be approximated by the weighted and shifted Grünwald difference (WSGD) operator in Deng’s paper[44] (in this paper $(p,q)=(1,0)$), for parameter $\gamma$, i.e.,

[TABLE]

[TABLE]

Let $u^{j}_{i}$ represent the numerical approximation of $u(x_{i},t_{j})$, then we denote that

[TABLE]

where

[TABLE]

and

[TABLE]

Suppose $u(x,t)\in\mathcal{C}^{4,3}_{x,t}[0,L]\times[0,T]$ is a solution of the TSADE (1.1). For simplicity, we introduce some symbols

[TABLE]

[TABLE]

Using (2.1)-(2.3), we shall see that the solution of (1.1) can be approximated by the following implicit difference scheme for $(x,t)=(x_{i},t_{j+\sigma})\in\bar{\omega}_{h\tau}$, $i=1,2,\ldots,N-1$, $j=0,1,\ldots,M-1$:

[TABLE]

Then we derive the implicit difference scheme with the approximation order $\mathcal{O}(\tau^{2}+h^{2})$:

[TABLE]

It is interesting to note that for $\theta\rightarrow 1$, we obtain the Crank-Nicolson difference scheme.

2.2 Analysis of the implicit difference scheme

In this subsection, we analyze the stability and convergence for the implicit difference scheme (2.4). Before proving the stability and convergence of the implicit difference scheme (2.4), we provide some properties of the coefficients $\omega_{k}^{(\alpha)}$ and $\omega_{k}^{(\beta)}$, which are derived from $g_{k}^{(\alpha)}$ and $g_{k}^{(\beta)}$ [1, 45, 2], respectively.

Lemma 2.2**.**

[45]** Suppose that $0<\alpha<1$, then the coefficients $\omega_{k}^{(\alpha)}$ satisfy

[TABLE]

Lemma 2.3**.**

[45]** Suppose that $1<\beta<2$, then the coefficients $\omega_{k}^{(\beta)}$ satisfy

[TABLE]

In the rest of this paper, we define the discrete inner product and the corresponding discrete $L_{2}$-norm as follows,

[TABLE]

where

[TABLE]

Based on the above lemmas, we will provide some other lemmas and corollaries, which are important properties to proof the stability and the convergence of the implicit difference scheme (2.4).

Lemma 2.4**.**

For $0<\alpha<1$, and any $\bm{v}\in V_{h}$, it holds that

[TABLE]

Proof.

According to the discrete definitions of $\delta^{\alpha}_{x,+}$ and $\delta^{\alpha}_{x,-}$, and $v_{0}=v_{N}=0$, we have

[TABLE]

and

[TABLE]

∎

Further, the following conclusion is true.

Corollary 1**.**

For $0<\alpha<1$, and any $\bm{v}\in V_{h},~{}N\geq 2$, there exists a positive constant $c_{1}$, such that

[TABLE]

Proof.

Since

[TABLE]

According to the nature of $g_{k}^{(\alpha)}$ [1, 45, 2], the following inequality is established:

[TABLE]

Then, there exist two positive constants $\tilde{c}_{1}$ and $c_{1}$, such that

[TABLE]

Here, the penultimate and antepenult inequalities are true due to the fact that $J_{1}(x)=1/x$ is a lower convex function and $J_{2}(x)=\ln(x)>0,~{}x\in[N,2N+1]$ is an increasing function.

Combining with Lemmas 2.2 and 2.4, we obtain

[TABLE]

∎

Using the same argument as in the proof of Corollary 1, we can easily carry out the following corollary.

Corollary 2**.**

(1) For $1<\beta<2$, and any $\bm{v}\in V_{h}$, it holds that

[TABLE]

(2) For $1<\beta<2$, and any $\bm{v}\in V_{h},~{}N\geq 3$, there exists a positive constant $c_{2}$, such that

[TABLE]

With the help of the preceding lemmas and corollaries, we can now establish the following theorem, which is essential for analyzing the stability of the proposed implicit difference scheme.

Theorem 2.1**.**

For any $\bm{v}\in V_{h}$, it holds that

[TABLE]

where $c$ is a positive constant independent of the spatial step size $h$.

Proof.

The concrete expression of $(\delta^{\alpha,\beta}_{h}\bm{v},\bm{v})$ can be written as

[TABLE]

According to Corollaries 1-2, there exist two positive constants $c_{1}$ and $c_{2}$ independent of the spatial step size $h$, such that for any non-vanishing vector $\bm{v}\in V_{h}$, we have

[TABLE]

For simplicity, we may take c=c_{1}\big{(}d_{+}^{j+\sigma}+d_{-}^{j+\sigma}\big{)}+c_{2}\big{(}e_{+}^{j+\sigma}+e_{-}^{j+\sigma}\big{)}. Hence, the targeted result is immediately completed. ∎

Now we can proof the stability and convergence of the implicit difference scheme (2.4). For convenience, in our proof, we denote $q_{s}^{j+1}=\frac{\tau^{-\theta}c_{j-s}^{(\theta,\sigma)}}{\Gamma(2-\theta)}$. Then \Delta_{0,t_{j+\sigma}}^{\theta}=\sum_{s=0}^{j}\big{(}u^{s+1}-u^{s}\big{)}q_{s}^{j+1}.

Theorem 2.2**.**

Denote $u^{j+1}=(u_{1}^{j+1},u_{2}^{j+1},\ldots,u_{N-1}^{j+1})^{T}$ and $\left\|f^{j+\sigma}\right\|^{2}=h\sum\limits_{i=1}^{N-1}f^{2}(x_{i},t_{j+\sigma})$. Then the implicit difference scheme (2.4) is unconditionally stable for $N\geq 3$, and the following a priori estimate holds:

[TABLE]

Proof.

Taking the inner product of (2.4) with $u^{(j+\sigma)}$, we have

[TABLE]

Using Corollary 1 in [43] and Theorem 2.1, obtains

[TABLE]

From (2.6), at $\varepsilon=c\ln 2$ we get

[TABLE]

Let us rewrite inequality (2.7) in the form

[TABLE]

noticing that $q_{0}^{j+1}>\frac{1}{2T^{\theta}\Gamma(1-\theta)}$ (cf.[47]).

With the inequality (2.8) in hand, the inequality (2.5) can be reached through utilizing the mathematical induction method.

The proof of Theorem 2.5 is completed. ∎

With the above proof, the convergence of the implicit difference scheme (2.4) is easy to obtain.

Theorem 2.3**.**

Suppose that $u(x,t)$ is the solution of (1.1) and $\{u_{i}^{j}\mid x_{i}\in\bar{\omega}_{h},~{}~{}0\leq j\leq M\}$ is the solution of the implicit difference scheme (2.4). Denote

[TABLE]

Then there exists a positive constant $\tilde{c}$ such that

[TABLE]

Proof.

It can easily obtain that $\xi^{j}$ satisfies the following error equation

[TABLE]

where $R_{i}^{j+\sigma}=\mathcal{O}(\tau^{2}+h^{2})$. By exploiting Theorem 2.5 we get

[TABLE]

which implies the convergence in the mesh $L_{2}$-norm with rate $\mathcal{O}(\tau^{2}+h^{2})$. ∎

The unconditionally stability and convergence of IDS (2.4) have been proved, and a numerical test is displayed to verify our results in Section 4.

2.3 A nonlinear extension of TSFADE

In this subsection, we extend the TSFADE to nonlinear case:

[TABLE]

where $g(u(x,t),x,t)=f(x,t)+y(u(x,t))$.

Considering (2.9) at point $(x_{i},t_{j+\sigma})$ and doing some simple manipulations, we have

[TABLE]

in which $g_{i}^{j+\sigma}=g(\sigma u_{i}^{j+1}+(1-\sigma)u_{i}^{j},x_{i},t_{j+\sigma})$, and the truncation error is $\mathcal{O}(\tau^{2}+h^{2})$.

We assume $U_{i}^{j}$ is the approximation solution of $u_{i}^{j}$ given in (2.10), and let $\xi_{i}^{j}=u_{i}^{j}-U_{i}^{j},~{}i=1,\cdots,N-1;~{}j=1,\cdots,M-1$, be the error satisfying

[TABLE]

Correspondingly, we denote $E^{j}=\left[\zeta_{1}^{j},\zeta_{2}^{j},\cdots,\zeta_{N-1}^{j}\right]^{T}~{}(j=1,\cdots,M-1)$, and suppose the nonlinear term $y(u(x,t))$ satisfies:

[TABLE]

with $0<L_{1}<\frac{a_{0}^{(\theta,\sigma)}}{2\left[\sigma^{2}+(1-\sigma)^{2}\right]\tau^{2}\Gamma(2-\theta)}$. Now, we proof the stability of (2.10).

Theorem 2.4**.**

The implicit difference scheme (2.10) is stable, i.e.:

[TABLE]

where $\tilde{C}$ is a positive constant, which may independent of $h$.

Proof.

Taking the inner product of (2.11) with $E^{(j+\sigma)}$, we have

[TABLE]

Utilizing Corollary 1 in [43] and noticing $(\delta^{\alpha,\beta}_{h}E^{(j+\sigma)},E^{(j+\sigma)})<0$, gets

[TABLE]

The right hand side can be estimated as below

[TABLE]

Bringing the above estimate to inequality (2.12), obtains

[TABLE]

Employing the mathematical induction method, the targeted result is immediately completed. ∎

Next, we study the convergence of (2.10).

Theorem 2.5**.**

Suppose that $u(x,t)$ is the solution of (2.9) and $\{u_{i}^{j}\mid x_{i}\in\bar{\omega}_{h},~{}~{}0\leq j\leq M\}$ is the solution of the implicit difference scheme (2.10). Denote

[TABLE]

Then there exists a positive constant $\hat{C}$ such that

[TABLE]

Proof.

Subtract (2.10) from (2.9), obtains

[TABLE]

where $\tilde{R}_{i}^{j+\sigma}=\mathcal{O}(\tau^{2}+h^{2})$. For convenience, we assume that there is a positive constant $\bar{C}$ such that: $\tilde{R}_{i}^{j+\sigma}\leq\bar{C}(\tau^{2}+h^{2}),~{}1\leq i\leq N-1,~{}0\leq j\leq M-1$.

Taking inner product of this error equation with $\psi^{j+\sigma}$, we have

[TABLE]

Similar as the proof of Theorem 2.2, we arrive at

[TABLE]

Utilizing the mathematical induction method, the priori estimate follows that

[TABLE]

The proof of this theorem is completed. ∎

The two assumptions about $y(u(x,t))$ are extremely strict, which is not very satisfactory. So our further work is to seek a more suitable method to prove the stability of IDS (2.10).

3 Fast implementation of IDS with the circulant precondtioner

Before moving into the investigation of fast solution techniques, the matrix form of implicit difference scheme (2.4) must be established first. In order to facilitate our discussion, we use notations in Section 2. Then, the matrix form of (2.4) corresponding to each time layer $j$ can be written as follows:

[TABLE]

we assume that the sums to be equal to zero if the upper summation index is less than the lower one, where

[TABLE]

[TABLE]

and $I$ is an identity matrix of order $N-1$, the two matrices $A_{\alpha}$ and $A_{\beta}$ are defined by

$A_{\alpha}=\begin{bmatrix}\omega_{1}^{(\alpha)}&\omega_{0}^{(\alpha)}&&&\\ \omega_{2}^{(\alpha)}&\omega_{1}^{(\alpha)}&\omega_{0}^{(\alpha)}&&\\ \vdots&\ddots&\ddots&\ddots&\\ \omega_{N-2}^{(\alpha)}&\cdots&\ddots&\ddots&\omega_{0}^{(\alpha)}\\ \omega_{N-1}^{(\alpha)}&\omega_{N-2}^{(\alpha)}&\cdots&\omega_{2}^{(\alpha)}&\omega_{1}^{(\alpha)}\end{bmatrix}$,     $A_{\beta}=\begin{bmatrix}\omega_{1}^{(\beta)}&\omega_{0}^{(\beta)}&&&\\ \omega_{2}^{(\beta)}&\omega_{1}^{(\beta)}&\omega_{0}^{(\beta)}&&\\ \vdots&\ddots&\ddots&\ddots&\\ \omega_{N-2}^{(\beta)}&\cdots&\ddots&\ddots&\omega_{0}^{(\beta)}\\ \omega_{N-1}^{(\beta)}&\omega_{N-2}^{(\beta)}&\cdots&\omega_{2}^{(\beta)}&\omega_{1}^{(\beta)}\end{bmatrix}$.

It is apparent that $A_{\alpha}$ and $A_{\beta}$ are Toeplitz matrices. Therefore, they can be stored with $2N+2$ entries. Krylov subspace methods with suitable circulant preconditioners [37, 48] can be used to efficiently solve Toeplitz or Toeplitz-like linear systems with a fast convergence rate. In this case, it also remarked that the algorithmic complexity of preconditioned Krylov subspace methods is only in $\mathcal{O}(N\log N)$ arithmetic operations per iteration step.

Inspired by the above consideration, we propose a circulant preconditioner to solve (3.1), which is generated from the Strang’s circulant preconditioner [49], through two preconditioned Krylov subspace methods. The Strang’s circulant matrix $s(Q)=[s_{j-k}]_{0\leq j,k<N}$ for a real Toeplitz matrix $Q=[q_{j-k}]_{0\leq j,k<N}$ is obtained by copying the central diagonals of $Q$ and bringing them around to complete the circulant requirement. More precisely, the diagonals of $s(Q)$ are given by

[TABLE]

Then our circulant preconditioner is defined as

[TABLE]

where the first columns of $s(A_{\gamma})$ and $s(A_{\gamma}^{T})$ ($\gamma=\alpha,\beta$) are given by

$\begin{bmatrix}\omega_{1}^{(\gamma)}\\ \vdots\\ \omega_{\left\lfloor\frac{N}{2}\right\rfloor}^{(\gamma)}\\ 0\\ \vdots\\ 0\\ \omega_{0}^{(\gamma)}\end{bmatrix}$ and $\begin{bmatrix}\omega_{1}^{(\gamma)}\\ \omega_{0}^{(\gamma)}\\ 0\\ \vdots\\ 0\\ \omega_{\left\lfloor\frac{N}{2}\right\rfloor}^{(\gamma)}\\ \vdots\\ \omega_{2}^{(\gamma)}\end{bmatrix}$, respectively.

To make sure the preconditioner defined in (3.2) is well-defined, let us illustrate that $P^{(j+\sigma)}$ are nonsingular. Before that, we need to give the following theorem, which is essential to check the nonsingularity of $P^{(j+\sigma)}$.

Theorem 3.1**.**

1. The real parts of all eigenvalues of $s(A_{\alpha})$ and $s(A_{\alpha}^{T})$ are strictly positive for all $N$;

2. The real parts of all eigenvalues of $s(A_{\beta})$ and $s(A_{\beta}^{T})$ are strictly negative for all $N$.

Proof.

Firstly, we proof the real parts of all eigenvalues of $s(A_{\alpha})$ are strictly positive for all $N$.

Recall that the real parts of all eigenvalues of $s(A_{\alpha})$ are equivalent to the eigenvalues of $\frac{s(A_{\alpha})+[s(A_{\alpha})]^{*}}{2}$ (where $S^{*}$ represents conjugate transpose of matrix $S$). Therefore, it is sufficient to show that all eigenvalues of $\frac{s(A_{\alpha})+s^{*}(A_{\alpha})}{2}$ are strictly positive for all $N$.

According to the Gershgorin circle theorem, all the Gershgorin disc of the circulant matrix $\frac{s(A_{\alpha})+s^{*}(A_{\alpha})}{2}$ are centered at $\omega_{1}^{(\alpha)}>0$ with radius

[TABLE]

by the Lemma 2.2, the first part of this theorem is proved.

Then, we prove the second part of the theorem. Similarly, all the Gershgorin disc of the circulant matrix $\frac{s(A_{\beta})+s^{*}(A_{\beta})}{2}$ are centered at $\omega_{1}^{(\beta)}<0$ with radius

[TABLE]

by the Lemma 2.3, the second part of this theorem is proved. Finally, we complete the proof of the theorem. ∎

With the help of Theorem 3.1, we will check the circulant preconditioners $P^{(j+\sigma)}$ are nonsingular.

Theorem 3.2**.**

The circulant preconditioners $P^{(j+\sigma)}$ defined in (3.2) are nonsingular.

Proof.

As we know, a circulant matrix can be diagonalized by the Fourier matrix F [49]. Then it follows that $s(A_{\gamma})=F^{*}\Lambda_{\gamma}F$, $s(A_{\gamma}^{T})=F^{*}\bar{\Lambda}_{\gamma}F$, where $\bar{\Lambda}_{\gamma}$ is the complex conjugate of $\Lambda_{\gamma}$. Decompose the circulant matrix $P^{(j+\sigma)}=F^{*}\Lambda_{P}F$ with the diagonal matrix

[TABLE]

Then the real part of $\Lambda_{P}$ can be written as

[TABLE]

Combining Theorem 3.1, we obtain

[TABLE]

Consequently, $P^{(j+\sigma)}$ are invertible. ∎

Unfortunately, due to the properties of coefficients $\omega_{k}^{(\alpha)}$ and $\omega_{k}^{(\beta)}$, it is difficult to theoretically investigate the eigenvalue distributions of preconditioned matrix $(P^{(j+\sigma)})^{-1}\mathcal{M}^{j+\sigma}$, but we still can give some figures to illustrate the clustering eigenvalue distributions of several specified preconditioned matrices in section 4.

At the end of this section, we study the fast implementation of IDS (2.10). For convenience, (2.10) is recast as

[TABLE]

in which

[TABLE]

Then the following algorithm is utilized to solve (3.3),

[TABLE]

with

[TABLE]

and

[TABLE]

At each iteration, the above system is linearized, and compared with (3.1), the system has not changed significantly. Thus, the previously mentioned fast method can be applied directly to this system. In the next section, we provide two experiments to test the effectiveness of our fast method, even for nonlinear situation.

4 Numerical results

In this section, we carry out several numerical experiments, which are given in Examples 1-2, to illustrate that our proposed IDM can indeed convergent with second order accuracy in both space and time. Some other numerical experiments, given in Examples 3-4, are reported to illustrate the effectiveness of the fast solution techniques. For direct solver, we choose LU factorization of MATLAB in Examples 3-4, and its CPU time is represented. For two Krylov subspace methods (BiCGSTAB and GPBiCOR($m$, $\ell$)) with circulant preconditioners, number of iterations required for convergence and CPU time of those methods are reported. We denote PBiCGSTAB and PGPBiCOR($m$, $\ell$) as the preconditioned version of BiCGSTAB and GPBiCOR($m$, $\ell$) respectively. The stopping criterion of those methods is

[TABLE]

where $r^{(k)}$ is the residual vector of the linear system after $k$ iterations, and the initial guess at each time step is chosen as the zero vector. All experiments were performed on a Windows 7 (32 bit) PC-Intel(R) Core(TM) i3-2130 CPU 3.40GHz, 4GB of RAM using MATLAB R2015b.

4.1 Verification of convergence rate

Example 1. We consider the linear problem (1.1), take coefficients

[TABLE]

and the source term

[TABLE]

Then the causal solution is $u(x,t)=t^{\theta+2}x^{2}(1-x)^{2}$.

As can be seen in Tables 1 and 3, when $h=1/3000$, the maximum error decreases steadily with the shortening of time step, and the convergence order of time is the expected $\mathcal{O}(\tau^{2})$, where the convergence order (CO) is given by the following formula: $\rm{CO}=\log_{\tau_{1}/\tau_{2}}\frac{\left\|\xi_{1}\right\|}{\left\|\xi_{2}\right\|}$ ($\xi_{i}$ is the error corresponding to $h_{i}$). On the other hand, Tables 2 and 4 illustrate that if $h=\tau$, a reduction in the maximum error occurs along with the decrease of space step and time step, and the spatial convergence order is $\mathcal{O}(h^{2})$, where the convergence order (CO) is given by the formula: $\rm{CO}=\log_{h_{1}/h_{2}}\frac{\left\|\xi_{1}\right\|}{\left\|\xi_{2}\right\|}$. Figs. 1-2 are plotted to further explain the reliability of our proposed scheme, and indicate that the ‘quadratic-type’ order of accuracy can achieve the desired order $\mathcal{O}(\tau^{2}+h^{2})$.

Example 2. In this example, we consider the nonlinear equation on space interval $[0,L]=[0,1]$ and time interval $[0,T]=[0,1]$ with advection coefficients $d_{+}(t)=sech(t)$, $d_{-}(t)=4sech(-t)$, diffusion coefficients $e_{+}(t)=(2+\cos(t))^{2}$, $e_{-}(t)=2+\cos^{2}(t)$, and nonlinear term

[TABLE]

in which $E_{\mu,\nu}(z)$ is the Mittag-Leffler function with two parameters defined by

[TABLE]

The exact solution is $u(x,t)=\left(t^{2+\theta}+\exp(2t)\right)x^{2}(1-x)^{2}$. For the finite difference discretization, the space step and time step are taken to be $h=1/N$ and $\tau=h$, respectively. The errors ($\xi=U-u$) and convergence order (CO) in the norms $\left\|\cdot\right\|_{0}$ and $\left\|\cdot\right\|_{\mathcal{C}(\bar{\omega}_{h\tau})}$, where $\left\|U\right\|_{\mathcal{C}(\bar{\omega}_{h\tau})}=$ max${}_{(x_{i},t_{j})\in\bar{\omega}_{h\tau}}|U|$, are given in Example 1.

As can be seen in Tables 5-8, the numerical solution provided by the difference approximation is in good agreement with our theoretical analysis. In Tables 5 and 7, we take $h=1/800$, the errors in maximum norm and $L_{2}$-norm decrease steadily with the shortening of time step, and the convergence order of time is the expected $\mathcal{O}(\tau^{2})$. While in Tables 6 and 8, the mesh size $\tau=h$ is chosen and the spatial convergence rates of the scheme (2.10) are also nearly two, for $\theta=0.1,0.5,0.9,0.99$. Furthermore, Figs. 3 and 4 are plotted to further illustrate the reliability of our proposed scheme, i.e., the slopes of the error curves in this log-log plot are 2, for $\theta=0.1,0.5,0.9,0.99$.

4.2 Fast implementations

Example 3. In this example, we consider the case of linear situation (1.1), and the exact solution is given in Example 1. In MATLAB code, we take the commands of generating matrices as the start point of Time1, and in the rest tables, “Speed-up” denotes ${\frac{\texttt{Time1}}{\texttt{Time}k}~{}(k=2,3)}$ and “Iter” is the average number of iterations required for solving the TSADE problem, namely,

[TABLE]

where $\texttt{Iter}(m)$ represents the number of iterations required for solving (3.1). As for GPBiCOR($m$, $\ell$) method [42] with circulant preconditioners in this example, a large number of experiments reveal that the results are better than others, when $m=3$ and $\ell=1$. Some eigenvalue plots about both original and preconditioned matrices are drawn in Figs. 5-6. These two figures confirm that for circulant preconditioning, the eigenvalues of preconditioned matrices are clustered at 1 except for a few of them. That is to say, the vast majority of the eigenvalues are well separated away from 0. We validate the effectiveness and robustness of the designed circulant preconditioner from the respective of clustering spectrum distribution.

The numerical results in Example 3 are presented in Tables 9-10, reflecting that the average iterations of the two Krylov subspace methods with circulant preconditioners have little difference. In other words, the number of average iterations of these methods is about 6.6. On the other hand, the CPU time by the two Krylov subspace methods with circulant preconditioners are much less than that by LU factorization method as $M$ and $N$ become large. When $M=N=2^{10}$ in Table 9, the CPU time is about 18.5 seconds, and the Speed-up is almost 7.4, especially for PBiCGSTAB method, the Speed-up can reach about 8.21. At the same case, in Table 10 the CPU time is about 18.7 seconds, and the Speed-up is almost 7.7. Meanwhile, although Time1 required by LU factorization method for small test problems ($M=N=2^{6},2^{7}$) is the cheapest among other two methods (i.e. PBiCGSTAB and PGPBiCOR($m$, $\ell$)), our proposed methods are still more attractive in aspects of the lower memory requirement.

Example 4. In this example, the exact solution and nonlinear source term are given in Example 2. For performing algorithm (3.4)-(3.6), the outer iteration tolerance is chosen as $10^{-12}$. As for GPBiCOR($m$, $\ell$) method [42] with circulant preconditioners in this example, a large number of experiments reveal that the results are better than others, when $m=3$ and $\ell=1$. Some eigenvalue plots about both original and preconditioned matrices are drawn in Figs. 7-8. These two figures confirm that for circulant preconditioning, the eigenvalues of preconditioned matrices are clustered at 1 except for a few of them. That is to say, the vast majority of the eigenvalues are well separated away from 0. We validate the effectiveness and robustness of the designed circulant preconditioner from the respective of clustering spectrum distribution.

In Tables 11-12, “Iter1” represents the average number of iterations required for solving the nonlinear problem with algorithm (3.4)-(3.6), i.e.,

[TABLE]

in which $\texttt{Iter1}(m)$ is the number of outer iteration required to solve (3.4), and “Iter2” denotes the average iterative number of a preconditioned Krylov subspace method required to solve (3.4), namely,

[TABLE]

in which $\texttt{Iter2}(k)$ is the iterative number of the preconditioned Krylov subspace method required at $k$th outer iteration. In view of the overall, the CPU time by the two Krylov subspace methods with circulant preconditioners are much less than that by LU factorization method as $M$ and $N$ become large. When $M=N=2^{10}$ in Table 11, the CPU time is about 41.8 seconds, and the Speed-up is almost 3.8. At the same case, Table 12 indicates that the CPU time is about 52.3 seconds, and the Speed-up is almost 3.2. Meanwhile, although Time1 required by LU factorization method for small test problems ($M=N=2^{6},2^{7},2^{8}$) is the cheapest among other two methods (i.e. PBiCGSTAB and PGPBiCOR($m$, $\ell$)), our proposed methods are still more attractive in aspects of the lower memory requirement. As a conclusion, our proposed IDS with fast solution techniques is still more practical than LU factorization method.

5 Conclusion

In this work, we constructed a difference scheme of the second order approximation in space and time for the linear TSADE with variable coefficients. Then, we utilized the same numerical approximation technique to simulate the nonlinear problem of TSADE and the convergence order still was $\mathcal{O}(\tau^{2}+h^{2})$. Numerical tests shown that our implicit difference schemes converge with second-order accuracy in both time and space. More significantly, we employed two reliable preconditioning iterative techniques, which have only computational cost and memory of $\mathcal{O}(N\log N)$ and $\mathcal{O}(N)$ respectively, to ameliorate the calculation skill. And numerical experiments strongly support the efficiency of the two preconditioning methods. In our future work, we will apply our IDM into more linear/nonlinear fractional partial differential equations and make a rigorous theoretical analysis of the stability for the developed difference schemes.

Acknowledgments

This work is supported by 973 Program (2013CB329404), NSFC (1370147, 11501085 and 61402082), the Fundamental Research Funds for the Central Universities (ZYGX2016J132 and ZYGX2016J138), the Scientific Research Fund of Sichuan Provincial Education Department (15ZA0288).

Appendix