Crank-Nicolson orthogonal spline collocation method combined with WSGI difference scheme for the two-dimensional time-fractional diffusion-wave equation

Abstract In this paper, a discrete orthogonal spline collocation method combining with a second-order Crank-Nicolson weighted and shifted Grünwald integral (WSGI) operator is proposed for solving time-fractional wave equations based on its equivalent partial integro-differential 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.


Introduction
Recently, fractional partial di erential 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 di erential equations. For example, when studying universal electromagnetic responses involving the uni cation of di usion and wave propagation phenomena, there are processes that are modeled by equations with time fractional derivatives of order γ ∈ ( , ) [1]. Generally, the analytical solutions of fractional partial di erential equations are di cult 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 nite element method [2,3], nite di erence method [4][5][6], meshless method [7,8], wavelets method [9],spline collocation method [10][11][12] and so forth.
Most of the numerical algorithms in [1][2][3][4][5][6][7][8] employed the L scheme to approximate fractional derivatives. Recently, Tian et al. [13] proposed second-and third-order approximations for Riemann-Liouville fractional derivative via the weighted and shifted Grünwald di erence (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 time-fractional sub-di usion equation. In [15], Chen considered the numerical solutions of the multi-term time fractional di usion and di usion-wave equations with variable coe cients, 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 two-dimensional distributed-order time fractional reaction-di usion equation. Following the idea of WSGD operator, Wang and Vong [17] used compact nite di erence WSGI scheme for the temporal Caputo fractional di usion-wave equation. However, the numerical methods with WSGI approximation have been rarely studied. Cao et al. [18] applied the idea of WSGI approximation combining with nite element method to solve the time fractional wave equation.
Orthogonal spline collocation (OSC) method has evolved as a valuable technique for solving di erent types of partial di erential equations [19][20][21][22][23]. The popularity of OSC is due to its conceptual simplicity, wide applicability and easy implementation. Comparing with nite di erence method and the Galerkin nite element method, OSC method has the following advantages: the calculation of the coe cients 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 e ciently [24]. Another feature of OSC method lies in its super-convergence [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 two-dimensional time-fractional wave equation, which is abbreviated as WSGI-OSC 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 WSGI-OSC scheme are presented in Section 4. Section 5 provides detailed description of the WSGI-OSC scheme. In Section 6, several numerical experiments are carried out to con rm the convergence analysis. Finally, the conclusion is drawn in Section 7.

Discrete-time OSC scheme . Preliminaries
In this section, we will introduce some notations and basic lemmas. For some positive integers Nx and Ny, δx and δy are two uniform partitions of I = [ , ] which are de ned as follows: δx : = x < x < · · · < x Nx = , δy : = y < y < · · · < y Ny = , . Let Mr(δx) and Mr(δy) be the space of piecewise polynomial of degree at most r ≥ , de ned by where Pr denotes the set of polynomial of degree at most r. It is easy to know that the dimension of the spaces For the functions u and v de ned on G, the inner product u, v and norm v Mr are respectively de ned by Throughout the paper, we denote C > a constant which is independent of mesh sizes h and τ. The following Young's inequality will also be used repeatedly, . Construction of the fully discrete orthogonal spline collocation scheme In this subsection, we consider discrete-time OSC schemes for solving the Eqs.(1.1)-(1.3). Our main idea of the proposed method is to transform the time fractional di usion-wave equation into its equivalent partial integro-di erential equation. To construct the continuous-time OSC scheme to the solution u of (1.1), we introduce the Riemann-Liouville fractional integral which is de ned by We integrate the equation(1.1) using Riemann-Liouville fractional integral operator I α t de ned in (2.7), then the problem is transformed into its equivalent partial integro-di erential equation as follows Let t k = kτ, k = , , · · · , N, where τ = T/N is the time step size. For the convenience of description, we de ne D t u n+ = u n+ −u n τ , and u n+ = u n+ +u n , where u n u(x, y, tn). Based on the idea of weighted and shifted Grünwald di erence operator, Wang and Vong ( [17]) established the second order accuracy approximation formula of the Riemann-Liouville fractional integral operator I α t u n+ , which is called as WSGI approximation, here By using the Crank-Nicolson di erence scheme and WSGI approximation formula to discretize the equation (2.8), we obtain the semi-discrete scheme in time direction For the needs of analysis, we give the following equivalent Galerkin weak formulation of the equation(2.12) by multiplying the equation with v ∈ H and integrating with respect to spatial domain Ω (2.14) We take the space Mr(δ) ⊂ H and obtain the fully discrete scheme as follows:

Stability and convergence analysis of WSGI-OSC scheme
In this section, we will give the stability and convergence analysis for fully-discrete WSGI-OSC scheme (2.13).
To this end, we further need the following lemmas. [29] Assume that kn and pn are nonnegative sequence, and the sequence ϕn satis es

Lemma 3.2. (Gronwall's ineqality)
where, g ≥ . Then the sequence ϕn satis es  (2.15) and applying the Cauchy-Schwarz inequality and Young inequality, it gives that Multiplying the above equation by τ, also using Lemma 1, then dropping the nonnegative terms Then, it gives that, Provided the time step τ is su ciently small, there exists a positive constant C such that The proof is complete.
where E n+ is de ned in (2.12). Taking v h = η n+ in (3.11), we have Next, we will give the estimate of I , I and I , respectively. (3.14) Taking advantages of mean value theorem and Cauchy-Schwarz inequality as well as Young inequality, we have tn ≤ t n+θ ≤ t n+ Using Lemma 1, we obtain and which completes the proof.

Description of the WSGI-OSC scheme
It can be observed from the fully discrete scheme (2.13) that we need to handle a two-dimensional partial di erential equation for each time level, that is We denote α = τ α+ λ (α) , β = τ α+ , then the above equation can be rewritten as are unknown coe cients to be determined. Settinĝ u = [û n , ,û n , , · · · ,û n ,My ,û n , ,û n , , · · · ,û n Mx ,My ] T , then the equation (4.2) can be written in the following form by Kronecker product and G n+ = [g n+ (ξ x , ξ y ), g n+ (ξ x , ξ y ), · · · , g n+ (ξ x , ξ y My ), g n+ (ξ x , ξ y ), · · · , g n+ (ξ x Mx , ξ y My )] T , (4.5) G n = [g n (ξ x , ξ y ), g n (ξ x , ξ y ), · · · , g n (ξ x , ξ y My ), g n (ξ x , ξ y ), g n (ξ x , ξ y ), · · · , g n (ξ x Mx , ξ y My )] T . (4.6) The matrices A x , B x , A y and B y are Mx × My having the following structure, We carry out the WSGI-OSC scheme in piecewise Hermite cubic spline space M (δ), which satis es zero boundary condition. Detailedly, we choose the basis of cubic Hermite polynomials [30], namly, for ≤ i ≤ K − , it follows that and In order to recover the coe cient matrix of the equations (4.3), we need to calculate the values of the basis functions at the Gauss point and their second-order derivatives. They are de ned as follows: It can be seen from the tensor product calculation that the WSGI-OSC scheme requires the solution of an almost block diagonal linear system at each time level, which can be solved e ciently by the software package COLROW [24].

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 = ) on uniform partitions of I in both x and y directions with Nx = Ny = 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 WSGI-OSC scheme, we present L∞ and L errors at T = and the corresponding convergence order de ned by where hm = /K is the time step size and em is the norm of the corresponding error.

Example 1
We consider the following one-dimensional time-fractional di usion-wave equation From the theoretical analysis, the numerical convergence order of WSGI-OSC (4.2) is expected to be O(τ + h ) when r = . 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 errors and the corresponding convergence orders of WSGI-OSC scheme for γ ∈ ( , ). 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 is xed. Numerical solution and global error for γ = . , h = / , τ = / are shown in Figure 1. .  .

Example 2 Consider the following one-dimensional fractional di usion-wave equation
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 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 , and γ = . , . , . , . , . , . . Table 4 presents the maximum L∞, L 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 = with γ = . , h = / , τ = / are shown in Figure 2.  .  .
where  Tables 5 and 6 list the maximum L∞, L 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 Tables 7 and 8 display L∞ and L errors and the corresponding convergence orders in time and space for some γ ∈ ( , ). Once again, the expected convergence rates with second-order accuracy in time direction . and fourth-order accuracy in spatial direction can be observed from two tables. Numerical solution and global error at T = with γ = . , h = / , τ = / are displayed in Figure 4. .  .

Conclusion
In this paper, we have constructed a Crank-Nicolson WSGI-OSC method for the two-dimensional timefractional di usion-wave equation. The original fractional di usion-wave equation is transformed into its equivalent partial integro-di erential equations, then Crank-Nicolson orthogonal spline collocation method with WSGI approximation is developed. The proposed method holds a higher convergence order than the convergence order O(τ −γ ) of general L approximation. The stability and convergence analysis are derived.
Some numerical examples are also given to con rm our theoretical analysis.