Approximate Solutions of Time Fractional Diffusion Wave Models

: In this paper, a wavelet based collocation method is formulated for an approximate solution of (1 + 1)- and (1 + 2)-dimensional time fractional diffusion wave equations. The main objective of this study is to combine the ﬁnite difference method with Haar wavelets. One and two dimensional Haar wavelets are used for the discretization of a spatial operator while time fractional derivative is approximated using second order ﬁnite difference and quadrature rule. The scheme has an excellent feature that converts a time fractional partial differential equation to a system of algebraic equations which can be solved easily. The suggested technique is applied to solve some test problems. The obtained results have been compared with existing results in the literature. Also, the accuracy of the scheme has been checked by computing L 2 and L ∞ error norms. Computations validate that the proposed method produces good results, which are comparable with exact solutions and those presented before.


Introduction
The theory of fractional calculus is an ancient topic that has many applications. However, practical work in this direction has been recently started (see References [1][2][3]). Most of the physical phenomena in chemistry, physics, engineering and other fields of science can be modeled using parameters of fractional calculus [4,5], means fractional derivative and integral operators. Amongst these are electrolyte polarization [6], viscoelastic systems [7], dielectric polarization [8] and so forth. Fractional models in different circumstances lead towards more accurate behaviour than those of integer order models.
The time fractional diffusion wave equation (TFDWE) is such an important model which has extensive uses. The TFDWE is actually a wave equation [9] with a fractional time derivative which describes universal acoustic, electromagnetic and mechanical responses [10,11] with an enhanced method. Over the past few decades, extensive attention has been paid to the closed form solution of time fractional diffusion wave equations (TFDWEs) and is still an open area of research. The closed form solution of such problems is not an easy job and needs herculean efforts. Owing to the fact several authors proposed numerical methods for the solution of fractional models, Tadjeran et al. [12] used second order accurate approximation for fractional diffusion equations. Zhuang et al. [13] applied an implicit numerical method for the anomalous sub-diffusion equation. Yuste and Acedo [14] studied fractional diffusion equations via an explicit finite difference method. Chen et al. [15] proposed the Fourier method for fractional diffusion equations. Hosseini et al. [16] solved the fractional telegraph equation with the help of radial basis functions. Zhou and Xu [17] applied the Chebyshev wavelets collocation method for the solution of time fractional diffusion wave equations. Bhrawya [18] used the spectral Tau algorithm based on the Jacobi operational matrix for the numerical solution of time fractional diffusion-wave equations. Yaseen et al. [19] solved fractional diffusion wave equations with reaction terms using finite differences and a trigonometric B-splines technique. Khader [20] and his co-author applied the finite difference method coupled with the Hermite formula for solutions of fractional diffusion wave equations. Kanwal et al. [21] implemented two-dimensional Genocchi Polynomials combined with the Ritz-Galerkin Method for solutions of fractional diffusion wave and Klein-Gordon equations. Datsko et al. [22] studied time-fractional diffusion-wave equation with mass absorption in a sphere under harmonic impact.
Recently, numerical methods using wavelets have been given more emphasis because of their simple applicability. These methods also have some other interesting properties such as the ability to detect singularities and express the function in different resolution levels, which improves the accuracy. Amongst different classes of wavelets, Haar wavelets deserve special consideration. Haar wavelets consist of piece wise constant functions. The integration of these wavelets in different times is one of the best features. Also, Haar wavelets have orthogonality and normalization properties with compact support. For more discussion on Haar wavelets one can see References [23,24].
In the present study, we propose a hybrid numerical scheme, based on Haar wavelets and finite differences, to solve (1 + 1)-and (1 + 2)-dimensional TFDWEs. The stability of the proposed method is discussed with the matrix method which is an essential part of the manuscript. The models which will be under consideration are characterized in the following types: (1 + 2)-Dimensional Equation: In Equations (1)-(4), ∆ is two-dimensional Laplacian; A, B, f , g, α, χ, κ, χ 1 are known functions and w is unknown function. Equations (2) and (4) are the corresponding initial and boundary conditions. The symbols, Ω and ∂Ω, Φ and ∂Φ represent the domain and boundary of the domain respectively for (1 + 1)-and (1 + 2)-dimensional problems. Also c D δ t w denotes the time fractional derivative of w with respect to t in the Caputo sense which is given by

Ground Work
In this section, some basic definitions of fractional calculus and Haar wavelets are presented, which will be required for the demonstration of our results. For a basic definition of Haar wavelets and its integrals we refer to Reference [23]. Let us consider x ∈ [a, b] where a and b are the limits of the interval. Next, the interval is subdivided into 2M intervals where M = 2 J and J denote the maximal level of resolution. Further, the two parameters j = 0, . · · · , J and k = 0, . · · · , 2 j − 1 are introduced. These parameters show the integer decomposition of wavelet number i = m + k + 1, where m = 2 j . The first and ith wavelets are defined as follows: where To solve nth order time fractional PDEs the following repeated integrals are needed: where β = 1, 2, . . . n, i = 1, 2, . . . 2M.
Keeping in view Equations (6) and (7) the close form expressions of these integrals are given by

Description of the Method
This section is devoted to discussing the scheme for Equations (1) and (3) separately. In both cases, the fractional order time derivative has been approximated by the quadrature formula [16] Case i: Using Equation (11) and θ−weighted scheme (0 ≤ θ ≤ 1) in Equation (1), we obtain After simplification, the above equation transforms to In our analysis we take θ = 1/2. Now approximating the highest order derivative by a truncated Haar wavelets series as: Integrating Equation (14) Integrating Equation (15) from 0 to 1, we get Substituting Equation (16) in Equation (15), the resultant equation reduces to Integration of Equation (17) from 0 to x yields Substituting values from Equations (14), (17) and (18) in Equation (13) and using collocation points x m = m−0.5 2M , m = 1, 2, . . . 2M, leads to the following system of algebraic equation where Equation (19) contains 2M equations. The unknown wavelet coefficients can be computed from this system. After determination of these unknown constants, the required solution at each time can be calculated from Equation (18).

Case ii:
Following a similar approach, as discussed earlier, Equation (3) gives Now we approximate w j+1 xxyy (x, y) with a two dimensional truncated Haar wavelets series as: where a j+1 i,l are unknowns to be determined. Integration of Equation (21) w.r.t. to y, between 0 and y, gives Integrating Equation (22) w.r.t y from 0 to 1, the unknown term w j+1 xxy (x, 0) is given by Substituting Equation (23) in Equation (22), the obtained result is Integrating Equation (24) from 0 to y, we get Repeating the same procedure one can easily derive the subsequent expressions Substitution of Equations (25), (26) and (29) in Equation (20) and using the collocation points, x m = m−0.5 2M , y n = n−0.5 2M , m, n = 1, 2, . . . 2M, produces the following system of equations where xx (x m , 1) Equation (30) represents 2M × 2M equations in so many unknowns which can be solved easily. After calculation of these unknowns, an approximate solution can be obtained from Equation (29).

Stability Analysis
Here we present the stability analysis of the proposed scheme for (1 + 2)-dimensional problems; a similar result can be proved for (1 + 1)-dimensional problems. In matrix form Equations (25), (26) and (29) can be written as w where α j+1 = α j+1 (i, l), U , V, Z andŨ j+1 ,Ṽ j+1 ,Z j+1 are interpolation matrices of w j+1 xx , w j+1 yy , w j+1 at collocation points and boundary terms, respectively. Now using Equations (31), (32) and (33) in Equation (20), we get Using Equation (33) in Equation (36) we have The above equation shows a recurrence relation of a full discretization scheme which allow us refinement in time. Ifw j+1 is numerical solution theñ Let e j+1 = w j+1 −w j+1 be the error at (j + 1) th time level. Subtracting Equation (37) from Equation (38) then e j+1 = Λe j , where Λ = Z C −1 T Z −1 is the amplification matrix. According to Lax-Richtmyer criterion, the scheme will be stable if Λ ≤ 1. It has been verified computationally that Λ ≤ 1. For J = 1 the spectral radius is 0.01025 which lies in the stability domain.

Convergence Analysis
The convergence analysis of scheme (18) and (29) is similar to the following theorems, therefore the proofs are omitted.

Theorem 2.
Assume w(x, y) and w 2M (x, y) be the exact and approximate solution of Equation (3), then Proof. See [27].

Illustrative Test Problems
In this part, we chose some test problem to confirm the reliability and efficiency of the present scheme. For validation of our results L ∞ and L 2 error norm are figured out which are defined as follows: where w app and w ext are respectively approximate and exact solutions.

Problem 5.1
Let us take the following (1 + 1)-dimensional TFDWE with damping Initial and boundary conditions are derived from the exact solution w(x, t) = t 2 x(1 − x). This problem has been solved for parameters J = 4, t = 0.01, 0.1, 1, δ = 1.1, 1.3, 1.5, 1.7, 1.9. The obtained error norms are shown in Table 1. From table it is obvious that results of the present scheme match well with exact solution. Also in Table 2 it has been observed that accuracy increases with increasing resolution level which shows the convergence in the spatial direction. In the same table, the results have been matched with existing results in the literature which clarify that computed solutions are in good agreement with the work of Chen et al. [28]. Table 3 shows convergence in time for fixed dx = 1/32. The convergence rate of the proposed scheme has been addressed in Table 4. the graphical solution and error plot are given in Figure 1. From this Figure  it is clear that approximate solutions are matchable with exact.

Problem 5.2:
Consider the following TFDWE with damping coupled with initial and boundary conditions The exact solution and source term are given by w(x, t) = e x t 3 and A(x, t) = 6t 3−δ e x Γ(4−δ) + 3t 2 e x − t 3 e x . In Table 5 the obtained error norms are shown for parameters t = 0.01, 0.1, δ = 1.1, 1.3, 1.5, 1.7, 1.9, J = 4. Table 5 shows that exact and approximate solutions agree with each other. The solution profile and absolute error are displayed Figure 2. From the Figure, the coincidence of both solutions are visible.

Problem 5.3:
Now we consider (1+2)-dimensional TFDWE [29] with exact solution w(x, y, t) = sin(πx)sin(πy)t δ+3 , and source term B(x, y, t) = sin(πx)sin(πy) We solved this problem for resolution level J = 4 and the obtained results are recorded in Table 6 for different values of time and τ. From Table 6 it is clear that the proposed scheme works well for the solution of two dimensional problems. Table 7 shows the comparison of the computed results with the previous work of Zhang [29]. One can see that our results are matchable with existing results. The same table presents convergence in time for (1 + 2)-dimensional problems. The graphical solution and absolute error of the problem are shown in Figure 3. It is obvious from Figure 3 that the exact and approximate solutions have strong agreement.

Problem 5.4:
Consider the following TFDWE with reaction term [19] coupled with initial and boundary conditions where the forcing terms are A( This problem has been solved with the help of the proposed scheme. In Table 8 we presented the solutions at different points. Also the obtained results have been compared with the work presented in Reference [19]. It is clear from table that our results are more accurate. From the table it is also obvious that the exact and numerical solutions are in good agreement. Exact verses numerical solutions are plotted in Figure 4. Graphical solutions also indicate that the proposed scheme works in the case where the reaction term exists. Exact and approximate Absolute error

Problem 5.5:
Now we consider the following equation where a 1 and b 1 are constants and the initial and boundary conditions are Here, we examine the behaviour of circular ring soliton numerically. Due to pulsating behaviour, such waves are also known as pulsons. We choose different values of parameters a 1 , b 1 to present surface plots to study the time evolution of the circular ring soliton. We observe the effect of a 1 and b 1 on solutions. In Figure 5, numerical solutions for different values of a 1 and b 1 have been plotted. Figure 6 shows the numerical solution for a 1 = 0.05 while varying b 1 . In Figure 7 the results are plotted for b 1 = 10, in which the wave peak value at the centre becomes lower as a 1 increases. This reveals that the solitary wave moves in a stable way up to a large time under finite initial condition.

Conclusions
In this paper, we proposed a hybrid method based on finite difference and Haar wavelets approximations. The scheme is applied for the numerical solution of (1 + 1)-and (1 + 2)-dimensional time fraction partial differential equations. The accuracy and applicability of the scheme is validated through some test problems. The tabulated data and graphical solution show that the scheme works very well for time fractional problems.