A numerical investigation of Caputo time fractional Allen–Cahn equation using redefined cubic B-spline functions

We present a collocation approach based on redefined cubic B-spline (RCBS) functions and finite difference formulation to study the approximate solution of time fractional Allen–Cahn equation (ACE). We discretize the time fractional derivative of order α ∈ (0, 1] by using finite forward difference formula and bring RCBS functions into action for spatial discretization. We find that the numerical scheme is of order O(h2 + t2–α) and unconditionally stable. We test the computational efficiency of the proposed method through some numerical examples subject to homogeneous/nonhomogeneous boundary constraints. The simulation results show a superior agreement with the exact solution as compared to those found in the literature.

In this paper, we consider the following time fractional ACE, which arises in mathematical modeling of phase separation in alloys of iron [1]: subject to the initial conditions u(z, 0) = φ 0 (z) ( 2 ) and boundary conditions where φ 0 (z) and ψ i (z) are supposed to be smooth functions with continuous first-order derivatives. There are many interpretations of fractional-order derivatives [16][17][18][19]. However, we have considered the fractional derivative in Caputo sense as [20] ∂ α u(z, t) where Γ is the usual gamma function. A lot of work on analytical solutions of time fractional ACE is available in the literature. Esen et al. [21] employed homotopy analysis method (HAM) to find the exact solution of time fractional ACE. The exact solution of space-time fractional ACE was discussed in [22] by means of a fractional subdiffusion method. The authors in [23] studied the coarsening of domains in a subdiffusive ACE in the context of the Seki-Lindenberg subdiffusion-reaction model. Yasar and Giresunlu [24] examined the exact solution of nonlinear space-time fractional ACE using (Ǵ G . 1 G )expansion algorithm. Guner et al. [25] solved the time fractional nonlinear ACE by means of the first integral, the exp-function and (Ǵ G )-expansion methods. Zhai et al. [26] applied a robust explicit operator splitting spectral approach to solve the nonlocal fractional Allen-Cahn model. The existence and uniqueness of weak solutions for the fractional Allen-Cahn, Cahn-Hilliard, and porous medium equations have been discussed by Akagi [27].
Tariq and Akram [1] proposed an analytical approach to determine the exact solutions for time fractional ACE. The authors applied the fractional complex transform method to reduce the equation into nonlinear ordinary differential equation (ODE). To explore numerical solution of space fractional ACEs, Hou et al. [28] employed the finite difference and Crank-Nicolson schemes for temporal and spatial discretization, respectively. Li et al. [29] explored the numerical solution of space-time fractional Allen-Cahn phase-field equation arising in the study of fluid mixture having two immiscible fluid phases. Some new exact solutions for time fractional ACE were presented by Hosseini et al. [30] using a new Kudryashov method in perspective of the conformable fractional derivative.
Most recently, Sakar et al. [31] proposed a numerical scheme based on the iterative reproducing kernel method (IRKM) to investigate the approximate solution of time fractional ACE. Liu et al. [32] developed a numerical scheme based on finite difference formulation and the Fourier spectral method to solve time fractional ACE in one and two space dimensions. Smooth and nonsmooth solutions for nonlinear space fractional ACE were studied by Yin et al. [33] using a fast algorithm based on time two-mesh finite difference method. Inc et al. [34] reduced the time fractional ACE and time fractional Klein-Gordon equations into the corresponding nonlinear fractional ODEs and employed an explicit power series method to solve these fractional ODEs. Khalid et al. [35] proposed a numerical approach based on cubic modified extended basis spline functions for time fractional diffusion wave equations. The authors in [36] employed nonpolynomial quintic spline functions for numerical investigation of fourth-order fractional boundary value problems involving product terms.
In this paper, we present a redefined cubic B-spline (RCBS) algorithm for numerical investigation of time fractional ACEs. The Caputo time fractional derivative has been discretized by finite difference formula, whereas spatial derivatives are discretized by RCBS functions. This approach is novel for numerical solution of fractional-order ACEs, and to the best of our knowledge, spline solution of time fractional ACE has never been studied yet. Moreover, this scheme is equally effective for homogeneous and nonhomogeneous boundary conditions. This paper is organized as follows. Section 2 evolves a brief description of temporal discretization, cubic B-spline functions, and spatial discretization. In Sect. 3, we discuss the stability of the proposed algorithm. The theoretical convergence analysis is presented in Sect. 4. The discussion on numerical results of three test problems has been reported in Sect. 5.

Description of the method 2.1 Temporal discretization
Let the time interval [0, T] be divided into M subintervals of equal length t = T M such as 0 = t 0 < t 1 < · · · < t M = T, where t n = n t, n = 0, 1, . . . , M. The Caputo definition given in Eq. (4) for time fractional derivative can be rewritten as Using forward difference formulation, Eq. (5) can be modified as where δ k = (k + 1) 1-α -(k) 1-α , ρ = t n+1τ , and the truncation error η n+1 t is bounded, where is a finite constant.

Cubic B-spline functions
In typical cubic B-spline collocation method, the approximation U * (z, t) to the exact solution u(z, t) is where γ j (t) are time-dependent unknowns to be computed. The cubic B-spline functions are twice differentiable at the nodal points z j over the domain [a, b] and preserve identical properties like nonnegativity, local support, convex, full property, symmetry, partition of unity, and geometric invariability. The blending function B j (z) of cubic B-spline can be represented in the following form [37]: where B -1 , B 0 , . . . , B N+1 can serve as spline basis over the domain [a, b]. The values of B j (z), B j (z), and B j (z) are given in Table 1 [38]. The approximate solution (U * ) n j = U * (z j , t n ) with its first-and second-order derivatives at the nth time level in terms of γ j can be written as

Redefined cubic B-spline functions
Usually, in collocation techniques, the Dirichlet type end conditions are imposed where the basis of spline functions vanish, but the CBS functions B -1 , B 0 , . . . , B N+1 do not vanish at the boundaries. So we redefine these bases so that they vanish at the boundaries when Dirichlet-type end conditions are prescribed [38,39]. The spline solution U(z, t) for the exact solution u(z, t) is obtained by eliminating γ n -1 and γ n N+1 from Eq. (10) as follows: where the weight function W (z, t) and RCBS B j (z) functions are given by Using Eq. (13) in Eq. (9), we obtain Let ( W * ) n+1 be the resultant term obtained from the weight functions. Then the last equation becomes Letting g n+1 = f n+1 -F n -( W * ) n+1 and using the values given in Table 1 with Eqs. (14) and (15) in Eq. (16), we form a recurrence relation after some simplifications: where ρ 0 = 6 h 2 , ρ 1 = ( r-1 6 -1 h 2 ), and ρ 2 = 2( r-1 3 + 1 h 2 ). The system of equations (17) in matrix form is given by where Using the initial conditions given in Eq. (2), we obtain the initial vector γ 0 = [γ 0 0 , γ 0 1 , . . . , γ 0 N ] T required to initiate the iteration process. We apply the initial conditions as follows: This gives N + 1 linear equations and can be written as where

Stability analysis
A numerical algorithm is considered to be stable if the errors do not propagate during the execution [40]. Here we employ the Fourier method to carry out the stability analysis of the proposed scheme. Let Φ n represent the growth factor in Fourier mode, and letΦ n be its approximate value. The error λ n i at the nth time level is given by where For simplicity, we study the stability in Eq. (17) for linear case (g = 0) only. Using Eqs. (21) and (17), we get the round-off error in the form The initial condition is satisfied by the error equation and similarly the boundary conditions become Define a grid function in the Fourier form: Now, in the Fourier series form, λ n (z) can be expressed as where Applying the norm, we obtain Hence Let λ n j = n e iνjh be the solution in the Fourier series form, where i = √ -1 and ν = 2π m b-a .

Theorem 1 The fully implicit scheme given in
Proof First of all, we have to prove by induction that | n | ≤ | 0 | for n = 0, . . . , T × M. For n = 0, Eq. (30) takes the form Let us assume that the result | n | ≤ | 0 | is true for n = 1, 2, . . . T × M -1. Now expression (30) can be written as Now, using Eqs. (28) and (31), we get This proves that the implicit scheme (17) is unconditionally stable.

Convergence analysis
We follow the technique used in [41] to investigate the uniform convergence of the proposed algorithm. Firstly, we state the following theorem [42,43].  N . LetŨ(z, t) be the unique spline approximation to the given problem at the spatial grid points z j ∈ , j = 0, . . . , N , then for all t ≥ 0, there exist μ j , independent of h, such that Proof By the triangle inequality we can write For grid points z j , we get Further, for a point z ∈ [z j , z j+1 ], we have

Theorem 3 Let U be the numerical approximation to the analytical exact solution u for
where Ω > 0 is a constant independent from h, and h is sufficiently small.
Proof Let u be the exact solution, and letŨ = N j=0 c j (t)B j be the spline approximation for U.
At the nth time stage, the present BVP can be written in the form of the difference equation L(Ũ(z j , t) -U(z j , t)): Also, the boundary constraints are and σ n j = h 2 g n jg j n , j = 0, . . . , N.
It is evident from (32) that σ n j = h 2 g n jg j n ≤ μh 4 .
Using the boundary conditions, we conclude that where μ 1 does not depend on h.
The proved theorem and relation (7) yield that the presented numerical scheme is convergent. Hence where and Ω are real constants, and α ∈ (0, 1]. This implies that the experimental order of convergence (EOC) of the proposed method is O(h 2 + t 2-α ).

Numerical experiments and discussion
In this section, we present three numerical experiments to check the efficiency of the proposed scheme for time fractional ACE. We test the computational results through error norms L 2 and L ∞ and the experimental order of convergence (EOC) as All computations are done by using MATHEMATICA 9.0 software.

Problem 1 Consider the time fractional ACE [31]
: The initial and boundary constraints can be extracted from the exact solution (z 2z) × t 1+α . A comparison of the maximum absolute error with IRKM [31] is presented in Tables 2, 3, 4, 5. The computational results are reported for α = 0.7, 0.9, t = 1, and t = 0.001 for 0.1 ≤ z ≤ 0.9. It is clear that we have achieved self-explanatory results as compared to the outcomes obtained from IRKM [31]. Table 6 reports the present experimental results at t = 1 with t = h 2 and α = 0.25, 0.5, 0.75. In Table 7, the error norms L 2 and L ∞ are tabulated for α = 0.2, 0.5, 0.8 with the variation of time t. Figure 1 depicts the physical behavior of exact and numerical solutions at different time levels for N = 80, α = 0.6, and t = 0.001 in the domain z ∈ [0, 1]. The three-dimensional plots given in Fig. 2 show the accuracy of the proposed scheme for α = 0.6, N = 80, t = 0.001, and t = 0.2. In Fig. 3 the spline solution at different time stages is displayed in a single frame when -1 ≤ z ≤ 2. The 3D plots of analytical exact and numerical solutions for N = 100, α = 0.4, t = 0.2, and t = 0.001 are shown in Fig. 4. Figure 5 shows the numerical results with variation of α for N = 20, t = 0.4, and 0 ≤ z ≤ 1. Problem 2 Consider the following time fractional ACE [29]: The source term f (z, t) on the right-hand side is given by    Table 4 Absolute errors for Example 1 by IRKM [31] at t = 1 when α = 0.9, N = 100, and t = 0.001     where E β,γ (ζ ), the Mittag-Leffler function, is defined as [29] E β,γ (ζ ) = ∞ k=0 ζ k Γ (βk + γ ) .
The initial and boundary conditions can be derived from the given exact solution z(1z 2 ) 3 t 2 E 1,3 (t). To affirm the impact and applicability of the proposed approach, the absolute error appeared in computations is given in Table 8 corresponding to different values of z using α = 0.5, n = 10, t = 0.1, t = 0.001, and 0 ≤ z ≤ 1. In Table 9 the L 2 norm is presented against different values of α and t for h -3 in the domain z ∈ [-1, 1]. The graphic representation of exact and approximate solutions is shown in Fig. 6 subject to N = 100, t = 0.001 α = 0.6, and -1 ≤ z ≤ 1. The 3D plots of exact and approximate solution for Example 2 are shown in Figs. 7-8.
Problem 3 Consider the following time fractional ACE [32]: where the forcing term on the right-hand side is given by   The initial and boundary conditions are u(z, 0) = φ 0 (z), u(a, t) = ψ 1 (t), u(b, t) = ψ 2 (t).
The exact solution of the given problem is t 2 sin(z). In Table 10 the exact and numerical solutions have been showcased for α = 0.9, t = 0.1, t = 0.001, and z ∈ [0, π]. Moreover, the absolute computational error has also been reported using the same choice of parameters.
Here we have intentionally taken n = 10 to demonstrate the efficiency of our method. In Figs. 9, 10, 11 the 2D and 3D plots of exact and numerical solutions are displayed. We can

Concluding remarks
We conclude this work by the following remarks: • An efficient algorithm based on redefined cubic B-spline collocation approach has been presented for approximate solution of time fractional Allen-Cahn equation. • The proposed scheme engages usual finite forward difference formulation and a redefined set of cubic B-spline functions for temporal and spatial discretizations, respectively. • The unconditional stability of the proposed method has been proved rigorously.
• The computational order of convergence is conformable with the theoretical estimations. • The numerical simulation has been run for three test examples, which show that the proposed scheme can efficiently be employed for numerical treatment of time fractional problems.