Abstract
In this paper, a discrete orthogonal spline collocation method combining with a secondorder CrankNicolson weighted and shifted Grünwald integral (WSGI) operator is proposed for solving timefractional wave equations based on its equivalent partial integrodifferential equations. The stability and convergence of the schemes have been strictly proved. Several numerical examples in one variable and in two space variables are given to demonstrate the theoretical analysis.
1 Introduction
Recently, fractional partial differential equations (FPDEs) have attracted more and more attention, which can be used to describe some physical and chemical phenomenon more accurately than the classical integerorder differential equations. For example, when studying universal electromagnetic responses involving the unification of diffusion and wave propagation phenomena, there are processes that are modeled by equations with time fractional derivatives of order γ ∈ (1, 2) [1]. Generally, the analytical solutions of fractional partial differential equations are difficult to obtain, so many authors have resorted to numerical solution techniques based on convergence and stability. Various kinds of numerical methods for solving FPDEs have been proposed by researchers, such as finite element method [2, 3], finite difference method [4, 5, 6], meshless method [7, 8], wavelets method [9], spline collocation method [10, 11, 12] and so forth.
In this study, we consider the following twodimensional timefractional diffusionwave equation
subject to the initial condition
and the boundary condition
where Δ is Laplace operator, Ω = [0, 1] × [0, 1] with boundary ∂ Ω, ϕ(x, y), φ(x, y) and f(x, y, t) are given sufficiently smooth functions in their respective domains and
in which Γ(⋅) is the Gamma function. Without loss of generality, we assume that ϕ(x, y) ≡ 0 in(1.2), since we can solve the equation for v(x, y, t) = u(x, y, t) − ϕ(x, y) in general.
Most of the numerical algorithms in [1, 2, 3, 4, 5, 6, 7, 8] employed the L_{1} scheme to approximate fractional derivatives. Recently, Tian et al. [13] proposed secondand thirdorder approximations for RiemannLiouville fractional derivative via the weighted and shifted Grünwald difference (WSGD) operators. Thereafter, some related research work covering the WSGD idea were done by many scholars. In [14], Liu et al developed a highorder local discontinuous Galerkin method combined with WSGD approximation for a Caputo timefractional subdiffusion equation. In [15], Chen considered the numerical solutions of the multiterm time fractional diffusion and diffusionwave equations with variable coefficients, which the time fractional derivative was approximated by WSGD operator. In [16], Yang proposed a new numerical approximation, using WSGD operator with second order in time direction and orthogonal spline collocation method in spatial direction, for the twodimensional distributedorder time fractional reactiondiffusion equation. Following the idea of WSGD operator, Wang and Vong [17] used compact finite difference WSGI scheme for the temporal Caputo fractional diffusionwave equation. However, the numerical methods with WSGI approximation have been rarely studied. Cao et al.[18] applied the idea of WSGI approximation combining with finite element method to solve the time fractional wave equation.
Orthogonal spline collocation (OSC) method has evolved as a valuable technique for solving different types of partial differential equations [19, 20, 21, 22, 23]. The popularity of OSC is due to its conceptual simplicity, wide applicability and easy implementation. Comparing with finite difference method and the Galerkin finite element method, OSC method has the following advantages: the calculation of the coefficients in the equation determining the approximate solution is fast since there is no need to calculate the integrals; and it provides approximations to the solution and spatial derivatives. Moreover, OSC scheme always leads to the almost block diagonal linear system, which can be solved by the software packages efficiently [24]. Another feature of OSC method lies in its superconvergence [25].
Motivated and inspired by the work mentioned above, the main goal of this paper is to propose a highorder OSC approximation method combined with second order WSGI operator for solving twodimensional timefractional wave equation, which is abbreviated as WSGIOSC in forthcoming sections. The remainder of the paper is organized as follows. In Section 2, some notations and preliminaries are presented. In Section 3, the fully discrete scheme combining WSGI operator with second order and orthogonal spline collocation scheme is formulated. Stability and convergence analysis of WSGIOSC scheme are presented in Section 4. Section 5 provides detailed description of the WSGIOSC scheme. In Section 6, several numerical experiments are carried out to confirm the convergence analysis. Finally, the conclusion is drawn in Section 7.
2 Discretetime OSC scheme
2.1 Preliminaries
In this section, we will introduce some notations and basic lemmas. For some positive integers N_{x} and N_{y}, δ_{x} and δ_{y} are two uniform partitions of I = [0, 1] which are defined as follows:
and
where P_{r} denotes the set of polynomial of degree at most r. It is easy to know that the dimension of the spaces M_{x}(δ_{x}) and M_{y}(δ_{y}) are (r − 1)N_{x} := M_{x} and (r − 1)N_{y} := M_{y}, respectively.
Let δ = δ_{x} ⊗ δ_{y} be a quasiuniform partition of Ω, and M_{r}(δ) = M_{r}(δ_{x}) ⊗ M_{r}(δ_{y}) with the dimension of M>_{x} × M_{y}. Let
as the sets of Gauss points in x and y direction, respectively, where
Let 𝓖 = {ξ = (ξ^{x}, ξ^{y}) : ξ^{x} ∈ 𝓖_{x}, ξ^{y} ∈ 𝓖_{y}}. For the functions u and v defined on 𝓖, the inner product 〈u, v〉 and norm ∥v∥_{Mr} are respectively defined by
For m a nonnegative integer, let H^{m}(Ω) denotes the usual Sobolev space with norm
where the norm ∥⋅∥ denotes the usual L_{2} norm, sometimes it is written as ∥⋅∥_{H0} for convenience. The following important lemmas are required in our forthcoming analysis. First, we introduce the differentiable (resp. twice differentiable) map W : [0, T] → M_{r}(δ) by
where u is the solution of the Eqs.(1.1)(1.3) . Then we have the following estimates for u − W and its time derivatives.
Lemma 2.1
[26] If ∂^{l} u/∂ t^{l} ∈ H^{r+3−j}, for all t ∈ [0, T], l = 0, 1, 2, j = 0, 1, 2, and W is defined by (2.1), then there exists a constant C such that
Lemma 2.2
[26] If ∂^{i} u/∂ t^{i} ∈ H^{r+3}, for t ∈ [0, T], i = 0, 1, then
where 0 ≤ l = l_{1} + l_{2} ≤ 4.
and there exists a positive constant C such that
Lemma 2.4
[28] The norms ∥⋅∥_{Mr} and ∥⋅∥ are equivalent on M_{r}(δ).
Throughout the paper, we denote C > 0 a constant which is independent of mesh sizes h and τ. The following Young's inequality will also be used repeatedly,
2.2 Construction of the fully discrete orthogonal spline collocation scheme
In this subsection, we consider discretetime OSC schemes for solving the Eqs.(1.1)(1.3). Our main idea of the proposed method is to transform the time fractional diffusionwave equation into its equivalent partial integrodifferential equation. To construct the continuoustime OSC scheme to the solution u of (1.1), we introduce the RiemannLiouville fractional integral which is defined by
where 0 < α = γ − 1 < 1.
We integrate the equation(1.1) using RiemannLiouville fractional integral operator
Let t_{k} = kτ, k = 0, 1, ⋯, N, where τ = T/N is the time step size. For the convenience of description, we define
where Ẽ = O(τ^{2}) and
here
By using the CrankNicolson difference scheme and WSGI approximation formula to discretize the equation (2.8), we obtain the semidiscrete scheme in time direction
where
For the needs of analysis, we give the following equivalent Galerkin weak formulation of the equation(2.12) by multiplying the equation with v ∈
We take the space M_{r}(δ) ⊂
3 Stability and convergence analysis of WSGIOSC scheme
In this section, we will give the stability and convergence analysis for fullydiscrete WSGIOSC scheme (2.13). To this end, we further need the following lemmas.
Lemma 3.1
[17] Let
Lemma 3.2
(Gronwall’s ineqality) [29] Assume that k_{n} and p_{n} are nonnegative sequence, and the sequence ϕ_{n} satisfies
where, g_{0} ≥ 0. Then the sequence ϕ_{n} satisfies
Theorem 3.1
The fullydiscrete WSGIOSC scheme (2.15) is unconditionally stable for sufficiently small τ > 0, it holds
Proof
Taking
Summing (3.2) for n from 0 to L(0 ≤ n ≤ N − 1), we obtain
Multiplying the above equation by 2τ, also using Lemma 1, then dropping the nonnegative terms
we have
Then, it gives that,
Provided the time step τ is sufficiently small, there exists a positive constant C such that
Using Gronwall’s Lemma 3.2, we get
The proof is complete.
Theorem 3.2
Suppose u is the exact solution of (1.1)(1.3), and
Proof
With W defined in (2.1), we set
thus we have
Because the estimate of η^{n} are provided by Lemma 2.2, it is sufficient to bound ζ^{n}, then use the triangle inequality to bound u^{n} −
where
Multiplying (3.12) by 2τ, and summing from n = 0 to n = L − 1 (1 ≤ n ≤ N + 1), it follows that
Next, we will give the estimate of I_{1}, I_{2} and I_{3}, respectively.
Taking advantages of mean value theorem and CauchySchwarz inequality as well as Young inequality, we have t_{n} ≤ t_{n+θ} ≤ t_{n+1}
Using Lemma 1, we obtain
Substituting (3.14), (3.15), (3.16) in (3.13) and removing the nonnegative terms, we attain
that is
Using the Gronwall’s inequality, Lemma 2.2 and triangle inequality, in the case that the time step τ is sufficiently small, there exists a positive constant C such that
and
which completes the proof.
4 Description of the WSGIOSC scheme
It can be observed from the fully discrete scheme (2.13) that we need to handle a twodimensional partial differential equation for each time level, that is
We denote
For applying the numerical schemes, firstly, we usually represent
then
where
then the equation (4.2) can be written in the following form by Kronecker product
where
and
The matrices A^{x}, B^{x}, A^{y} and B^{y} are M_{x} × M_{y} having the following structure,
We carry out the WSGIOSC scheme in piecewise Hermite cubic spline space M_{3}(δ), which satisfies zero boundary condition. Detailedly, we choose the basis of cubic Hermite polynomials [30], namly, for 1 ≤ i ≤ K − 1, it follows that
and
Note that functions ϕ_{i}(x), ψ_{i}(x) satisfy zero boundary conditions ϕ_{i}(0) = ϕ_{i}(1) = ψ_{i}(0) = ψ_{i}(1) = 0. Renumber the basis functions and let
then
In order to recover the coefficient matrix of the equations (4.3), we need to calculate the values of the basis functions at the Gauss point and their secondorder derivatives. They are defined as follows:
where
It can be seen from the tensor product calculation that the WSGIOSC scheme requires the solution of an almost block diagonal linear system at each time level, which can be solved efficiently by the software package COLROW [24].
5 Numerical experiments
In this section, four examples are given to demonstrate our theoretical analysis. In our implementations, we adopt the space of piecewise Hermite bicubics(r = 3) on uniform partitions of I in both x and y directions with N_{x} = N_{y} = K. The forcing term f(x, y, t) is approximated by the piecewise Hermite interpolant projection in the Guass points. To check the accuracy of WSGIOSC scheme, we present L_{∞} and L_{2} errors at T = 1 and the corresponding convergence order defined by
where h_{m} = 1/K is the time step size and e_{m} is the norm of the corresponding error.
Example 1
We consider the following onedimensional timefractional diffusionwave equation
where
From the theoretical analysis, the numerical convergence order of WSGIOSC (4.2) is expected to be O(τ^{2} + h^{4}) when r = 3. In order to check the second order accuracy in time direction, we select τ = h so that the error caused by the spatial approximation can be negligible. Table 1 lists L_{∞} and L_{2} errors and the corresponding convergence orders of WSGIOSC scheme for γ ∈ (1, 2). We observe that our scheme generates the temporal accuracy with the order 2. To test the spatial approximation accuracy, Table 2 shows that our scheme has the accuracy of 4 in spatial direction, where the temporal step size τ = h^{2} is fixed. Numerical solution and global error for γ = 1.3, h = 1/32, τ = 1/32 are shown in Figure 1.
Figure 1
γ  τ  L_{∞} error  Convergence order  L_{2} error  Convergence order 

1.1 

7.0727×10^{−5}  4.4681×10^{−5}  

1.7932×10^{−5}  1.9798  1.1012×10^{−5}  2.0206  

4.5623×10^{−6}  1.9747  2.7487×10^{−6}  2.0022  

1.1483×10^{−6}  1.9903  6.8758×10^{−7}  1.9992  
1.3 

2.6081×10^{−4}  1.7238×10^{−4}  

6.6648×10^{−5}  1.9684  4.2518×10^{−5}  2.0194  

1.6825×10^{−6}  1.9860  1.0577×10^{−5}  2.0072  

4.2263×10^{−7}  1.9931  2.6387×10^{−6}  2.0030  
1.5 

4.1657×10^{−4}  2.7911×10^{−4}  

1.0633×10^{−4}  1.9701  6.8593×10^{−5}  2.0247  

2.6736×10^{−5}  1.9916  1.7020×10^{−5}  2.0108  

6.7115×10^{−6}  1.9941  4.2405×10^{−6}  2.0050  
1.7 

5.3422×10^{−4}  3.6265×10^{−4}  

1.3701×10^{−4}  1.9632  8.9343×10^{−5}  2.0212  

3.4419×10^{−5}  1.9930  2.2160×10^{−5}  2.0114  

8.6292×10^{−6}  1.9959  5.5175×10^{−6}  2.0059  
1.9 

5.7600×10^{−4}  3.9339×10^{−4}  

1.4884×10^{−4}  1.9523  9.7244×10^{−5}  2.0163  

3.7391×10^{−5}  1.9930  2.4112×10^{−5}  2.0119  

9.3633×10^{−6}  1.9976  5.9996×10^{−6}  2.0068  
1.95 

5.6941×10^{−4}  3.8862×10^{−4}  

1.4696×10^{−4}  1.9540  9.6061×10^{−5}  2.0163  

3.6917×10^{−5}  1.9931  2.3812×10^{−5}  2.0123  

9.2425×10^{−6}  1.9979  5.9232×10^{−6}  2.0072 
γ  h  L_{∞} error  Convergence order  L_{2} error  Convergence order 

1.1 

2.4371×10^{−6}  1.7740×10^{−6}  

1.5377×10^{−7}  3.9863  1.0837×10^{−7}  4.0329  

9.6290×10^{−9}  3.9972  6.6928×10^{−9}  4.0172  

6.0225×10^{−10}  3.9989  4.1576×10^{−10}  4.0088  
1.3 

3.8377×10^{−6}  2.6750×10^{−6}  

2.4364×10^{−7}  3.9774  1.6332×10^{−7}  4.0338  

1.5241×10^{−8}  3.9987  1.0085×10^{−8}  4.0174  

9.5308×10^{−10}  3.9992  6.2644×10^{−10}  4.0089  
1.5 

4.7527×10^{−6}  3.2535×10^{−6}  

3.0159×10^{−7}  3.9781  1.9851×10^{−7}  4.0347  

1.8863×10^{−8}  3.9990  1.2256×10^{−8}  4.0177  

1.1798×10^{−9}  3.9990  7.6129×10^{−10}  4.0089  
1.7 

5.1530×10^{−6}  3.4857×10^{−6}  

3.2579×10^{−7}  3.9834  2.1258×10^{−7}  4.0354  

2.0382×10^{−8}  3.9986  1.3123×10^{−8}  4.0178  

1.2754×10^{−9}  3.9982  8.1509×10^{−10}  4.0090  
1.9 

4.6730×10^{−6}  3.0735×10^{−6}  

2.9311×10^{−7}  3.9948  1.8735×10^{−7}  4.0361  

1.8412×10^{−8}  3.9927  1.1563×10^{−8}  4.0181  

1.1509×10^{−9}  3.9999  7.1819×10^{−10}  4.0090  
1.95 

4.3316×10^{−6}  2.8280×10^{−6}  

2.7151×10^{−7}  3.9958  1.7235×10^{−7}  4.0364  

1.7062×10^{−8}  3.9922  1.0637×10^{−8}  4.0182  

1.0665×10^{−9}  3.9999  6.6066×10^{−10}  4.0091 
Example 2
Consider the following onedimensional fractional diffusionwave equation
where
In order to test the temporal accuracy of the proposed method, we choose τ = h to avoid contamination of the spatial error. The maximum L_{∞}, L_{2} errors and temporal convergence orders are shown in Table 3. To check the convergence order in space, the time step τ and space step h are chosen such that τ = h^{2}, and γ = 1.1, 1.3, 1.5, 1.7, 1.9, 1.95. Table 4 presents the maximum L_{∞}, L_{2} errors and spatial convergence orders. The results in Tables 3 and 4 demonstrate the expected convergence rates of 2 order in time and 4 order in space simultaneously. Numerical solution and global error at T = 1 with γ = 1.5, h = 1/32, τ = 1/32 are shown in Figure 2.
Figure 2
γ  τ  L_{∞} error  Convergence order  L_{2} error  Convergence order 

1.1 

2.7779×10^{−5}  1.8686×10^{−5}  

6.9405×10^{−6}  2.0009  4.5452×10^{−6}  2.0395  

1.7225×10^{−6}  2.0105  1.1135×10^{−6}  2.0292  

4.2704×10^{−7}  2.0121  2.7427×10^{−7}  2.0215  
1.3 

6.8399×10^{−5}  4.6042×10^{−5}  

1.6912×10^{−5}  2.0159  1.1079×10^{−5}  2.0551  

4.1818×10^{−6}  2.0158  2.7032×10^{−6}  2.0352  

1.0358×10^{−6}  2.0134  6.6503×10^{−7}  2.0232  
1.5 

1.0251×10^{−4}  6.9519×10^{−5}  

2.5114×10^{−5}  2.0292  1.6555×10^{−5}  2.0701  

6.2025×10^{−6}  2.0176  4.0327×10^{−6}  2.0375  

1.5384×10^{−6}  2.0114  9.9335×10^{−7}  2.0214  
1.7 

1.4424×10^{−4}  9.9816×10^{−5}  

3.5642×10^{−5}  2.0168  2.4001×10^{−5}  2.0562  

8.8717×10^{−6}  2.0063  5.8868×10^{−6}  2.0275  

2.2116×10^{−6}  2.0041  1.4572×10^{−6}  2.0143  
1.9 

1.9061×10^{−4}  1.2932×10^{−4}  

4.7290×10^{−5}  2.0110  3.1561×10^{−5}  2.0347  

1.1810×10^{−5}  2.0015  7.7937×10^{−6}  2.0178  

2.9519×10^{−6}  2.0003  1.9367×10^{−6}  2.0087  
1.95 

2.0110×10^{−4}  1.3455×10^{−4}  

4.9930×10^{−5}  2.0099  3.2930×10^{−5}  2.0306  

1.2482×10^{−5}  2.0001  8.1369×10^{−6}  2.0169  

3.1218×10^{−6}  1.9994  2.0222×10^{−6}  2.0085 
γ  h  L_{∞} error  Convergence order  L_{2} error  Convergence order 

1.1 

9.8873×10^{−8}  7.9806×10^{−8}  

6.3013×10^{−9}  3.9719  5.0124×10^{−9}  3.9929  

4.0052×10^{−10}  3.9757  3.1726×10^{−10}  3.9818  

2.5425×10^{−11}  3.9775  2.0139×10^{−11}  3.9776  
1.3 

2.0590×10^{−7}  1.1973×10^{−7}  

1.2348×10^{−8}  4.0596  7.0248×10^{−9}  4.0912  

7.5090×10^{−10}  4.0395  4.2273×10^{−10}  4.0547  

4.6084×10^{−11}  4.0263  2.5826×10^{−11}  4.0328  
1.5 

3.1378×10^{−7}  1.8224×10^{−7}  

1.9140×10^{−8}  4.0351  1.0838×10^{−8}  4.0716  

1.1827×10^{−9}  4.0165  6.6114×10^{−10}  4.0350  

7.3513×10^{−11}  4.0079  4.0836×10^{−11}  4.0170  
1.7 

4.2637×10^{−7}  2.5157×10^{−7}  

2.6414×10^{−8}  4.0127  1.5170×10^{−8}  4.0516  

1.6453×10^{−9}  4.0049  9.3246×10^{−10}  4.0241  

1.0270×10^{−10}  4.0019  5.7825×10^{−11}  4.0113  
1.9 

6.2873×10^{−7}  3.5996×10^{−7}  

3.9276×10^{−8}  4.0007  2.1898×10^{−8}  4.0389  

2.4548×10^{−9}  4.0000  1.3510×10^{−9}  4.0188  

1.5342×10^{−10}  4.0000  8.3898×10^{−11}  4.0092  
1.95 

6.7414×10^{−7}  3.9216×10^{−7}  

4.2262×10^{−8}  3.9956  2.3894×10^{−8}  4.0367  

2.6423×10^{−9}  3.9995  1.4745×10^{−9}  4.0184  

1.6515×10^{−10}  3.9999  9.1572×10^{−11}  4.0091 
Example 3
Consider the following twodimensional fractional diffusionwave equation
where
Similar to the selection of parameters in Examples 1 and 2, Tables 5 and 6 list the maximum L_{∞}, L_{2} errors and convergence orders, respectively. The similar convergence rates in time and space are also obtained. As we hope, the convergence order of all numerical results match that of the theoretical analysis. Figure 3 plots the numerical solution and global error at T = 1 with γ = 1.7, h = 1/32, τ = 1/32.
Figure 3
γ  N  L_{∞} error  Convergence order  L_{2} error  Convergence order 

1.1  10  1.6611×10^{−4}  8.3486×10^{−5}  
15  7.5461×10^{−5}  1.9461  3.7876×10^{−5}  1.9493  
20  4.2909×10^{−5}  1.9624  2.1509×10^{−5}  1.9669  
25  2.7589×10^{−5}  1.9792  1.3841×10^{−5}  1.9754  
1.3  10  5.1729×10^{−4}  2.6164×10^{−4}  
15  2.3249×10^{−4}  1.9724  1.1769×10^{−4}  1.9704  
20  1.3148×10^{−4}  1.9813  6.6585×10^{−5}  1.9799  
25  8.4565×10^{−5}  1.9779  4.2760×10^{−5}  1.9847  
1.5  10  7.7899×10^{−4}  3.9475×10^{−4}  
15  3.4829×10^{−4}  1.9853  1.7651×10^{−4}  1.9850  
20  1.9648×10^{−4}  1.9899  9.9627×10^{−5}  1.9882  
25  1.2607×10^{−4}  1.9886  6.3896×10^{−5}  1.9905  
1.7  10  9.8958×10^{−4}  5.0659×10^{−4}  
15  4.3990×10^{−4}  1.9995  2.2433×10^{−4}  2.0090  
20  2.4748×10^{−4}  1.9995  1.2609×10^{−4}  2.0028  
25  1.5841×10^{−4}  1.9994  8.0685×10^{−5}  2.0006  
1.9  10  1.1985×10^{−3}  6.2808×10^{−4}  
15  5.3173×10^{−4}  2.0044  2.7856×10^{−4}  2.0052  
20  2.9891×10^{−4}  2.0022  1.5657×10^{−4}  2.0026  
25  1.9123×10^{−4}  2.0015  1.0018×10^{−4}  2.0014 
γ  N  L_{∞} error  Convergence order  L_{2} error  Convergence order 

1.1  10  1.7277×10^{−6}  5.8164×10^{−7}  
15  3.7129×10^{−7}  3.7921  1.1583×10^{−7}  3.9799  
20  1.2237×10^{−7}  3.8582  3.6753×10^{−8}  3.9902  
25  5.1343×10^{−8}  3.8922  1.5074×10^{−8}  3.9942  
1.3  10  4.5383×10^{−6}  2.0806×10^{−6}  
15  8.9315×10^{−7}  4.0091  4.1188×10^{−7}  3.9946  
20  2.8185×10^{−7}  4.0092  1.3042×10^{−7}  3.9974  
25  1.1523×10^{−7}  4.0083  5.3439×10^{−8}  3.9984  
1.5  10  7.1532×10^{−6}  3.3974×10^{−6}  
15  1.4118×10^{−6}  4.0021  6.7176×10^{−7}  3.9976  
20  4.4624×10^{−7}  4.0036  2.1262×10^{−7}  3.9988  
25  1.8263×10^{−7}  4.0037  8.7104×10^{−8}  3.9993  
1.7  10  9.2188×10^{−6}  4.4527×10^{−6}  
15  1.8187×10^{−6}  4.0031  8.7956×10^{−7}  4.0000  
20  5.7483×10^{−7}  4.0038  2.7830×10^{−7}  3.9999  
25  2.3526×10^{−7}  4.0036  1.1399×10^{−7}  3.9999  
1.9  10  1.1444×10^{−5}  5.7230×10^{−6}  
15  2.2505×10^{−6}  4.0110  1.1299×10^{−6}  4.0011  
20  7.1020×10^{−7}  4.0091  3.5746×10^{−7}  4.0005  
25  2.9046×10^{−7}  4.0068  1.4641×10^{−7}  4.0003 
Example 4
Consider the following twodimensional fractional diffusionwave equation
where
Tables 7 and 8 display L_{∞} and L_{2} errors and the corresponding convergence orders in time and space for some γ ∈ (1, 2). Once again, the expected convergence rates with secondorder accuracy in time direction and fourthorder accuracy in spatial direction can be observed from two tables. Numerical solution and global error at T = 1 with γ = 1.9, h = 1/32, τ = 1/32 are displayed in Figure 4.
Figure 4
γ  N  L_{∞} error  Convergence order  L_{2} error  Convergence order 

1.1  10  8.8381×10^{−4}  4.4449×10^{−4}  
15  3.9978×10^{−4}  1.9566  2.0073×10^{−4}  1.9607  
20  2.2738×10^{−4}  1.9615  1.1385×10^{−4}  1.9711  
25  1.4625×10^{−4}  1.9775  7.3238×10^{−5}  1.9771  
1.3  10  3.1514×10^{−3}  1.5847×10^{−3}  
15  1.4225×10^{−3}  1.9617  7.1426×10^{−4}  1.9654  
20  8.0814×10^{−4}  1.9656  4.0464×10^{−4}  1.9752  
25  5.1939×10^{−4}  1.9811  2.6009×10^{−4}  1.9807  
1.5  10  5.3058×10^{−3}  2.6680×10^{−3}  
15  2.3861×10^{−3}  1.9709  1.1981×10^{−3}  1.9745  
20  1.3534×10^{−3}  1.9711  6.7766×10^{−4}  1.9808  
25  8.6906×10^{−4}  1.9851  4.3519×10^{−4}  1.9847  
1.7  10  7.2062×10^{−3}  3.6236×10^{−3}  
15  3.2347×10^{−3}  1.9755  1.6242×10^{−3}  1.9792  
20  1.8321×10^{−3}  1.9760  9.1737×10^{−4}  1.9857  
25  1.1754×10^{−3}  1.9893  5.8858×10^{−4}  1.9889  
1.9  10  8.0346×10^{−3}  4.0402×10^{−3}  
15  3.6198×10^{−3}  1.9665  1.8175×10^{−3}  1.9701  
20  2.0516×10^{−3}  1.9736  1.0273×10^{−3}  1.9833  
25  1.3162×10^{−3}  1.9893  6.5910×10^{−4}  1.9888 
γ  N  L_{∞} error  Convergence order  L_{2} error  Convergence order 

1.1  10  1.2725×10^{−5}  6.5169×10^{−5}  
15  2.5700×10^{−6}  3.9453  1.2858×10^{−5}  4.0029  
20  8.0847×10^{−7}  4.0202  4.0665×10^{−5}  4.0014  
25  3.3300×10^{−7}  3.9751  1.6653×10^{−5}  4.0009  
1.3  10  3.6079×10^{−5}  1.8230×10^{−5}  
15  7.1968×10^{−6}  3.9758  3.6064×10^{−6}  3.9964  
20  2.2773×10^{−6}  3.9997  1.1417×10^{−6}  3.9982  
25  9.3472×10^{−7}  3.9907  4.6773×10^{−7}  3.9989  
1.5  10  5.7773×10^{−5}  2.9133×10^{−5}  
15  1.1491×10^{−5}  3.9829  5.7620×10^{−6}  3.9968  
20  3.6402×10^{−6}  3.9959  1.8240×10^{−6}  3.9984  
25  1.4930×10^{−6}  3.9942  7.4725×10^{−7}  3.9991  
1.7  10  7.6488×10^{−5}  3.8541×10^{−5}  
15  1.5190×10^{−5}  3.9868  7.6189×10^{−6}  3.9981  
20  4.8133×10^{−6}  3.9948  2.4113×10^{−6}  3.9991  
25  1.9734×10^{−6}  3.9959  9.8780×10^{−7}  3.9994  
1.9  10  8.4694×10^{−5}  4.2666×10^{−5}  
15  1.6803×10^{−5}  3.9892  8.4290×10^{−6}  3.9997  
20  5.3240×10^{−6}  3.9951  2.6670×10^{−6}  3.9999  
25  2.1823×10^{−6}  3.9968  1.0924×10^{−6}  4.0000 
6 Conclusion
In this paper, we have constructed a CrankNicolson WSGIOSC method for the twodimensional timefractional diffusionwave equation. The original fractional diffusionwave equation is transformed into its equivalent partial integrodifferential equations, then CrankNicolson orthogonal spline collocation method with WSGI approximation is developed. The proposed method holds a higher convergence order than the convergence order O(τ^{3−γ}) of general L_{1} approximation. The stability and convergence analysis are derived. Some numerical examples are also given to confirm our theoretical analysis.
Acknowledgement
The authors would like to thank the referees for their very helpful and detailed comments, which have significantly improved the presentation of this paper. This work was supported by the National Natural Science Foundation of China (Grant No.11601076) and the Ph.D. Research Startup Fund Project of East China University of Technology (Grant No.DHBK2019213).
References
[1] A.A. Kilbas, H.M. Srivastava, and J.J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier, Amsterdam, 2006. Search in Google Scholar
[2] Y.D. Zhang, Y.M. Zhao, F.L. Wang, and Y.F. Tang, Highaccuracy finite element method for 2D time fractional diffusionwave equation on anisotropic meshes, Int. J. Comput. Math. 95 (2018), no. 1, 218–230. 10.1080/00207160.2017.1401708 Search in Google Scholar
[3] L.M. Li, D. Xu, and M. Luo, Alternating direction implicit Galerkin finite element method for the twodimensional fractional diffusionwave equation, J. Comput. Phys. 255 (2013), 471–485, 10.1016/j.jcp.2013.08.031 Search in Google Scholar
[4] Z.Z. Sun and X. Wu, A fully discrete difference scheme for a diffusionwave system, Appl. Numer. Math. 56 (2006), no. 2, 193–209, 10.1016/j.apnum.2005.03.003 Search in Google Scholar
[5] J.F. Huang, Y.F. Tang, L. Vázquez, and J.Y. Yang, Two finite difference schemes for time fractional diffusionwave equation, Numer. Algorithms 64 (2013), no. 4, 707–720, 10.1007/s1107501296890 Search in Google Scholar
[6] Y.N. Zhang, Z.Z. Sun, and X. Zhao, Compact alternating direction implicit scheme for the twodimensional fractional diffusionwave equation, SIAM J. Numer. Anal. 50 (2012), no. 3, 1535–1555, 10.1137/110840959 Search in Google Scholar
[7] M. Aslefallah and E. Shivanian, An efficient meshless method based on RBFs for the time fractional diffusionwave equation, Afr. Mat. 29 (2018), no. 78, 1203–1214, 10.1007/s133700180616y Search in Google Scholar
[8] M. Dehghan, M. Abbaszadeh, and A. Mohebbi, Analysis of a meshless method for the time fractional diffusionwave equation, Numer. Algorithms 73 (2016), no. 2, 445–476, 10.1007/s1107501601031 Search in Google Scholar
[9] M.H. Heydari, M.R. Hooshmandasl, F.M. Maalek Ghaini, and C. Cattanic, Wavelets method for the time fractional diffusionwave equation, Phys. Lett. A 379 (2015), no. 3, 71–76, 10.1016/j.physleta.2014.11.012 Search in Google Scholar
[10] G. Fairweather, X.H. Yang, D. Xu, and H.X. Zhang, An ADI CrankNicolson orthogonal spline collocation method for the twodimensional fractional diffusionwave equation, J. Sci. Comput. 65 (2015), no. 3, 1217–1239, 10.1007/s109150150003x Search in Google Scholar
[11] M. Yaseen, M. Abbas, T. Nazir, and D. Baleanu, A finite difference scheme based on cubic trigonometric Bsplines for a time fractional diffusionwave equation, Adv. Differ. Equ. 2017 (2017), 274, 10.1186/s136620171330z Search in Google Scholar
[12] A. Esen, O. Tasbozan, Y. Ucar, and N.M. Yagmurlu, A Bspline collocation method for solving fractional diffusion and fractional diffusionwave equations, Tbilisi Math. J. 8 (2015), no. 2, 181–193, 10.1515/tmj20150020 Search in Google Scholar
[13] W.Y. Tian, H. Zhou, and W.H. Deng, A class of second order difference approximations for solving space fractional diffusion equations, Math. Comput. 84 (2015), no. 294, 1298–1314, 10.1090/S002557182015029172 Search in Google Scholar
[14] Y. Liu, M. Zhang, H. Li, and J.C. Li, Highorder local discontinuous Galerkin method combined with WSGDapproximation for a fractional subdiffusion equation, Comput. Math. Appl. 73 (2017), no. 6, 1298–1314, 10.1016/j.camwa.2016.08.015 Search in Google Scholar
[15] H. Chen, S. Lü, and W. Chen, A unified numerical scheme for the multiterm time fractional diffusion and diffusionwave equations with variable coefficients, J. Comput. Appl. Math. 330 (2018), 380–397, 10.1016/j.cam.2017.09.011 Search in Google Scholar
[16] X. Yang, H. Zhang, and D. Xu, WSGDOSC Scheme for twodimensional distributed order fractional reactiondiffusion equation, J. Sci. Comput. 76 (2018), no. 3, 1502–1520, 10.1007/s1091501806723 Search in Google Scholar
[17] Z. Wang and S. Vong, Compact difference schemes for the modified anomalous fractional subdiffusion equation and the fractional diffusionwave equation, J. Comput. Phys. 277 (2014), 1–15, 10.1016/j.jcp.2014.08.012 Search in Google Scholar
[18] Y. Cao, B.L. Yin, Y. Liu, and H. Li, CrankNicolson WSGI difference scheme with finite element method for multidimensional timefractional wave problem, Comput. Appl. Math. 37 (2018), no. 4, 5126–5145, 10.1007/s4031401806262 Search in Google Scholar
[19] B. Bialecki and G. Fairweather, Orthogonal spline collocation methods for partial differential equations, J. Comput. Appl. Math. 128 (2001), no. 12, 55–82, 10.1016/S03770427(00)005094 Search in Google Scholar
[20] C.E. GreenwellYanik and G. Fairweather, Analyses of spline collocation methods for parabolic and hyperbolic problems in two space variables, SIAM J. Numer. Anal. 23 (1986), no. 2, 282–296, 10.1137/0723020 Search in Google Scholar
[21] C. Li, T.G. Zhao, W.H. Deng, and Y.J. Wu, Orthogonal spline collocation methods for the subdiffusion equation, J. Comput. Appl. Math. 255 (2014), 517–528, 10.1016/j.cam.2013.05.022 Search in Google Scholar
[22] L. Qiao and D. Xu, Orthogonal spline collocation scheme for the multiterm timefractional diffusion equation, Int. J. Comput. Math. 95 (2018), no. 8, 1478–1493, 10.1080/00207160.2017.1324150 Search in Google Scholar
[23] H.X. Zhang, X.H. Yang, and D. Xu, A highorder numerical method for solving the 2D fourthorder reactiondiffusion equation, Numer. Algorithms 80 (2019), no. 3, 849–877, 10.1007/s110750180509z Search in Google Scholar
[24] G. Fairweather and I. Gladwell, Algorithms for almost block diagonal linear systems, SIAM Rev. 46 (2004), no. 1, 49–58, 10.1137/S003614450240506X Search in Google Scholar
[25] A.K. Pani, G. Fairweather, and R.I. Fernandes, ADI orthogonal spline collocation methods for parabolic partial integrodifferential equations, IMA J. Numer. Anal. 30 (2010), no. 1, 248–276, 10.1093/imanum/drp024 Search in Google Scholar
[26] C.E. GreenwellYanik and G. Fairweather, Analyses of spline collocation methods for parabolic and hyperbolic problems in two space variables, SIAM J. Numer. Anal. 23 (1996), no. 2, 282–296, 10.1137/0723020 Search in Google Scholar
[27] M.P. Robinson and G. Fairweather, Orthogonal spline collocation methods for Schrödingertype equations in one space variable, Numer. Math. 68 (1994), no. 3, 355–376, 10.1007/s002110050067 Search in Google Scholar
[28] B. Li, G. Fairweather, and B. Bialecki, Discretetime orthogonal spline collocation methods for Schrödinger equations in two space variables, SIAM J. Numer. Anal. 35 (1998), no. 2, 453–477, 10.1137/S0036142996302396 Search in Google Scholar
[29] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, Springer, Berlin, 1997. Search in Google Scholar
[30] S. Arora, I. Kaur, H. Kumar, and V.K. Kukreja, A robust technique of cubic hermite collocation for solution of two phase nonlinear model, Journal of King Saud University – Engineering Sciences 29 (2017), no. 2, 159–165, 10.1016/j.jksues.2015.06.003 Search in Google Scholar
© 2020 Xiaoyong Xu and Fengying Zhou, published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.