Non-Polynomial Quintic Spline for Numerical Solution of Fourth--Order Time Fractional Partial Differential Equations

This paper presents a novel approach for numerical solution of a class of fourth order time fractional partial differential equations (PDE's). The finite difference formulation has been used for temporal discretization, whereas, the space discretization is achieved by means of non polynomial quintic spline method. The proposed algorithm is proved to be stable and convergent. In order to corroborate this work, some test problems have been considered and the computational outcomes are compared with those found in the exiting literature. It is revealed that the presented scheme is more accurate as compared to current variants on the topic.


Introduction
In the modern era, fractional order differential equations have gained a significant amount of research work due to their vide range of applications in various branches of science and engineering such as Physics, electrical networks, fluid mechanics, control theory, theory of viscoelasticity, neurology and theory of electromagnetic acoustics [1,2]. Wang [3] introduced the very first approximate solution of nonlinear fractional Korteweg-de Vries (KdV) Burger equation involving space and time fractional derivatives using Adomian Decomposition method. Zurigat at al. [4] examined the approximate solution of fractional order algebraic differential equations using Homotopy analysis method. Turut and Guzel [5] implemented Adomian decomposition method and multivariate Pade approximation method for solving fractional order nonlinear partial differential equations (PDE's). In [6], Liu and Hou applied the Generalized differential transform method to solve the coupled Burger equation with space and time fractional derivatives. Khan et al. [7] used Adomian decomposition method and Variational iteration method for numerical solution of fourth order time fractional PDE's with variable coefficients. Later on, Abbass et al. [8] employed a finite difference approach based on third degree trigonometric B-spline functions for approximate solution of one-dimensional wave equation. Javidi and Ahmad [9] developed a computational technique based on Homotopy perturbation method Laplace transform and Stehfests numerical inversion algorithm for solving fourth-order time-fractional PDE's with variable coefficients.The fractional differential transform method and modified fractional differential transform method were proposed by Kanth and Aruna [10] for series solution to higher dimensional third-order dispersive fractional PDE's. Pandey and Mishra [11] applied Sumudu transforms and Homotopy analysis approach for solving time-fractional third order dispersive type of PDE's. The fractional Variational iteration method was put into action by Prakash and Kumar in [12] for series solution to third-order fractional dispersive PDE's in higher dimensional space.
The spline approximation techniques have been applied extensively for numerical solution of ODE's and PDE's. The spline functions have a variety of significant gains over finite difference schemes. These functions provide a continuous differentiable estimation to solution over the whole spatial domain with great accuracy. The straightforward employment of spline functions provides a solid ground for applying them in the context of numerical approximations for initial/boundary problems. Khan and Aziz [13], solved third order boundary-value problems (BVP's) using a numerical method based on quintic spline functions. In [14], non polynomial quintic spline method was employed for numerical solution of fourth order two-point BVP's. Khan and Sultana [15] proposed non-polynomial quintic spline functions for numerical solution of third order BVP's associated with odd-order obstacle problems. In [16], Srivastava discussed numerical solution of differential equations using polynomial spline functions of different orders. Siddiqi and Arshed [17] brought the fifth degree basis spline collocation functions into use for approximate solution of fourth order time fractional PDE's. Rashidinia and Mohsenyzadeh [18] used non-polynomial quintic spline technique for one-dimensional heat and wave equations. Recently, in [19], fifth degree spline approximation technique has been utilized for approximate solution of fourth-Order time-fractional PDE's. In [20], the new fractional order spline functions were considered to obtain the approximate solution for fractional Bagely-Torvik Equation. Arshed [21] employed quintic B-spline collocation scheme for solving fourth order time-fractional super diffusion equation. More recently, parametric quintic spline approach and Grunwald-Letnikov approximation have been proposed in [22] for a distributed order fractional sub-diffusion problem.
In the field of modern science and engineering the fourth-order initial/boundary value problems are of great importance. For example, airplane wings, bridge slabs, floor systems and window glasses are being modeled as plates subject to different types end supports which are successfully described in terms of fourth-order PDE's [19]. In this work, we consider the following class of the fourth-order time-fractional with the following initial and boundary conditions where γ ∈ (0, 1), is the order of fractional time derivative, α represents the ratio of flexural-rigidity of beam to its mass per unit length, y(x, t) is the beam transverse displacement, u(x, t) describes the dynamic driving force per unit mass and the function v 0 (x) is known to be continuous on [0, L]. There are many descriptions to the concept of fractional differentiation but Caputo and Riemann-Liouville have been the most common definitions. Here, we shall use the Caputo's approach because it is more appropriate for real world problems and it permits initial and boundary conditions in terms of ordinary derivatives. The Caputo's definition of fractional derivative of order γ is given by This paper has been composed with the aim to develop a spline collocation method for approximate solution of fourth order time-fractional PDE's. The backward Euler's scheme has been utilized for temporal discretization, whereas, non polynomial quintic spline function, comprised of a trigonometric part and a polynomial part, has been used to interpolate the unknown function in spatial direction. The presented technique has also been proved to be stable and convergent.
This work is arranged as follows: In section 2, a brief explanation of quintic spline scheme has been presented and the consistency relations between the values of spline approximation and its derivatives at the nodal points are derived. Section 3 describes the use of L1 approximation in time direction to achieve a backward Euler technique. Non-polynomial quintic spline scheme for the spatial discretization has been discussed in section 4. The computational results and discussions are given in section 5.

Description of Non Polynomial Quintic Spline Function
Consider x i = ih, be the mesh points of uniform partition of [0, L] into sub-intervals [x i , x i−1 ], where h = L n and i = 0, 1, 2, · · · , n. Let y(x) be a sufficiently smooth function defined on [0, L]. We denote the non polynomial quintic spline approximation to y(x) by S(x). Each non polynomial spline segment R i (x) has the following form where a i , b i , c i , d i , e i and f i are the constants and the parameter ξ, the frequency of the trigonometric functions, will be used to enhance the accuracy of the technique. When ξ approaches to zero, Eq.(2.1) reduces to quintic polynomial spline function in [a, b]. The non polynomial quintic spline can be defined as (2.2) First of all, we establish the consistency relations for all the coefficients involved in (2.1) in terms of S i 's, The values of coefficients introduced in (2.1) can be calculated as The relation (2.7) provides (n − 3) linear equations with (n − 1) unknowns S i , i = 1(1)n − 1. Hence, we require two more equations for direct calculation of S i , one at each end of the range of integration, which can be formulated as setting i = 1, 2 in (2.4) we have Similarly, subtracting (2.11) from (2.9), we get Now, the first end condition is obtained by substituting (2.12), (2.13) into (2.8) for i = 1.
Similarly, the second end condition for i = n, is given by Lemma 2.1. The local truncation error t i , i = 1(1)n − 1 associated with the Eqs (2.7), (2.14) and (2.15) is given by Proof. We have to find local truncation error t i , i = 1, 2, ..., n − 1 for the present scheme. First of all, we write Eqs (2.7), (2.14), and (2.15) as for τ = 4, 5, 6, 7, we get The local truncation error given in Eq (2.16) takes the following form

Temporal Discretization
In order to discretize the time fractional derivative, backward Euler scheme is employed. We consider, t p = p∆t for p = 0(1)K with ∆t = T K as the step size in time direction. The computation of Caputo time-fractional derivative at t = t p+1 can be made as The above equation along with the definition of Caputo fractional derivative gives the following relation.
Now, we define a semi-discrete fractional differential operator G γ t as Then, Eq. (3.1) can be written as Here, l p+1 ∆t denotes the truncation error between ∂ γ ∂t γ y(x, t p+1 ) and G γ t y(x, t p+1 ). Let G γ t y(x, t p+1 ) be the approximation of Caputo time-fractional derivative at t = t p+1 , then Eq. (1.1) can be expressed as Using (3.1), the above equation can be written as where, β = Γ(2 − γ)∆t γ and y p+1 (x) = y(x, t p+1 ) with the initial and boundary conditions as follow Moreover, the coefficients b j involved in (3.1) have the following properties where the constant c is dependant on y. To apply this scheme, we need the values y 0 and y 1 .
For p = 0, (3.4) takes the following form Now (3.4) and (3.6) with initial and boundary conditions formulate a complete set of semi-discrete problem for (1.1).
The error term l p+1 can also be defined as [23] From Eqs.(3.2) and (3.5), the error term can be expressed as Now, we define some functional spaces and their standard norms as x , ∀r ≤ n} where L 2 (η) denotes the space of all measurable functions whose square is Lebesgue integrable in η . The inner product and norm in L 2 (η) are given by The inner product and norm in S 2 (η) are given by Also, the norm . in H n (η) is defined in the following way It is also preferred to define . 2 g 2 = g 2 0 + βα g (2) x 2 0 1 2 (3.9) Now, for the stability and convergence analysis, we are to find y p+1 ∈ H 2 0 (η) such that for all g ∈ H 2 0 (η), Eqs.(3.4) and (3.6) give the following two relations and < y 1 , g > +βα < y 1 xxxx , g >=< y 0 , g > +β < u 1 , g > (3.11) The theorem given below describes the unconditional stability of the semi-discrete problem.
Theorem 1. The discrete problem is unconditionally stable in such a way that ∀∆t > 0, it holds where . 2 is discussed in Eq.(3.9).
Proof. In order to prove this result, mathematical induction is used. For p = 0 and g = y 1 , Eq. (3.11) takes the following form < y 1 , y 1 > +βα < y 1 xxxx , y 1 >=< y 0 , y 1 > +β < u 1 , y 1 > Integrating by parts, the above result can be written as Due to the boundary conditions on g, all the boundary related contributions are disappeared. From Schwarz inequality and the inequality g 0 ≤ g 2 , Eq. (3.13) becomes Suppose that the result is true for g = y j i.e y j 2 ≤ y 0 0 + β j i=1 u i 0 , j = 2, 3, · · · , p. (3.14) Let g = y p+1 in Eq.(3.10) Integrating by parts, we get Again due to the boundary conditions on g all the boundary related contributions are disappeared. From Schwarz inequality and the inequality g 0 ≤ g 2 , the above expression changes to Using (3.14), the above relation takes the following form Using the properties of b j , we can write < e 1 , g > +βα < e 1 xx , g xx >=< e 0 , g > + < l 1 , g >, ∀g ∈ S 2 0 (η).
Let g = e 1 and e 0 = 0 gives the following relation For p = 1, Eq. (3.17) is satisfied.
Since, for p = 1, p −γ b −1 p−1 = 1. Therefore, it can be written in the following form Therefore, ∀ p such that p∆t ≤ T , The above discussion can be summed up in following theorem.
The operator Φ is defined as Now, Eq.(2.7) takes the following form Applying the operator Φ on Eq.(4.1), we get the following result , p = 1, 2, 3, · · · , K − 1. (4.4) After simplifying the system (4.4) takes the following form = Q i , i = 2, 3, · · · , n − 2, p = 1, 2, · · · , K − 1. (4.5) where System (4.6) provides (n− 3) equations involving S p+1 i , i = 1, 2, · · · , n− 1. Therefore, we further need two equations for complete solution of S p+1 i . The required two end conditions can be derived using simply supported boundary conditions as The proposed algorithm is a five point scheme. In order to implement it, the numerical values of S 2 = [S 2 1 , S 2 2 , S 2 3 , · · · , S 2 n−1 ] T and S 1 = [S 1 1 , S 1 2 , S 1 3 , · · · , S 1 n−1 ] T are needed. To calculate the values of S 2 , it is required to find S 1 . Solving Eq.(3.6) and using the non polynomial quintic spline technique, value of S 1 can be found as: where The system (4.9) consists of (n − 3) equations involving S 1 i , i = 1, 2, · · · , n − 1. Hence, to get a unique solution to this system, two additional end equations can be obtained from simply supported boundary conditions in the following way are column vectors with dimension (n − 1). The system in (4.9)-(4.11) can be expressed as where A,B and C are square matrices of order (n − 1), such that

Calculation of Truncation Error
The Eq.(4.4) can be written in the following form Expanding equation (4.4) with Taylor series in terms of S(x i , t p ) and its spatial derivatives , the truncation error is obtained,as Where From above discussion and Theorem 2, It is concluded that the scheme is of O(h 4 + ∆t 4−γ )

Numerical Results
In this section, we consider three test problems to check the validity and efficiency of the proposed numerical scheme. The approximate results are compared with quintic spline collocation method (QnSM) used in [19]. All the computations are executed in M athematica 9.0. The accuracy of presented technique is tested by error norms L ∞ , L 2 and order of convergence (χ), which are calculated as where y i , Y i represent the exact and approximate solution at i th knot respectively.
Problem 1. Consider the fourth order time-fractional PDE [19] ∂ γ y ∂t γ + α with initial condition y(x, 0) = sin(πx) and the boundary conditions The exact solution is y(x, t) = sin(πx)e t . The computational error norms L ∞ and L 2 corresponding to different values of γ are listed in Table 1 when α = 0.01 and n = 100. It is obvious that our proposed computational approach produces more accurate results with ∆t = 0.01 as compared to QnSM used in [19] with ∆t = 0.000001. A comparison of L ∞ , L 2 and order of convergence χ with QnSM [19] at t = 1 corresponding to γ = 0.5 and ∆t = h is reported in Table 2. It is observed that the order of convergence in numerical results exhibits a good agreement with the theoretical estimation. In Figure 1, three dimensional visuals of exact and approximate solutions are displayed for n = 100, ∆t = 0.01. The absolute numerical error at t = 1 corresponding to n = 100, ∆t = 0.01 and γ = 0.5 is portrayed in Figure   2.     1.00 1.0000 × 10 −6 7.0711 × 10 −8 9.8911 × 10 −9 9.0927 × 10 −9 Problem 2. Consider the following fourth order time-fractional PDE [19] ∂ γ y ∂t γ + 0.05 with initial condition y(x, 0) = 0 and the boundary conditions The close form solution is y(x, t) = t sin(πx). The computational error norms L ∞ and L 2 for n = 100 and different choices of γ are presented in Table 3. It can be observed that our presented approach yields more accurate results with ∆t = 0.01 as compared to QnSM employed in [19] with ∆t = 0.000001. Table 4 presents the error norms L ∞ , L 2 and the corresponding order of convergence at t = 1 for ∆t = handγ = 0.5. Figure 3 shows three dimensional plots of exact and approximate solutions for n = 100, ∆t = 0.01 and γ = 0.5. The absolute numerical error at t = 1 corresponding to n = 100, ∆t = 0.01 and γ = 0.25 is portrayed in Figure 4.    Problem 3. Consider the fourth order time-fractional PDE [19] ∂ γ y ∂t γ + 0.05 The analytical exact solution is y(x, t) = (t + 1) sin πx. Table 5 shows a comparison of computational error norms with QnSM [19] corresponding to different selections of γ. It is found that our computational outcomes are better than QnSM [19]. In Figure 5, three dimensional visuals of exact and approximate solutions are displayed for n = 40, γ = 0.5 and ∆t = 0.01. The absolute numerical error at t = 1 corresponding to n = 40, ∆t = 0.01 and γ = 0.5 is portrayed in Figure 6. It is obvious that approximate solution is highly consistent with the analytical exact solution, which proves the effectiveness of proposed scheme.
(a) Exact solution (b) Approximate solution

Conclusion
In this work, non polynomial quintic spline collocation method has been employed for approximate solution of fourth order time-fractional partial differential equations. The backward Euler's method has been used for temporal discretization, whereas, non polynomial quintic spline function composed of a trigonometric part and a polynomial part has been employed to interpolate the unknown function in spatial direction. The proposed numerical algorithm is proved to be convergent and unconditionally stable. The numerical outcomes are found to be more accurate as compared to QnSM [19].