Using Fractional Bernoulli Wavelets for Solving Fractional Diffusion Wave Equations with Initial and Boundary Conditions

In this paper, fractional-order Bernoulli wavelets based on the Bernoulli polynomials are constructed and applied to evaluate the numerical solution of the general form of Caputo fractional order diffusion wave equations. The operational matrices of ordinary and fractional derivatives for Bernoulli wavelets are set via fractional Riemann–Liouville integral operator. Then, these wavelets and their operational matrices are utilized to reduce the nonlinear fractional problem to a set of algebraic equations. For solving the obtained system of equations, Galerkin and collocation spectral methods are employed. To demonstrate the validity and applicability of the presented method, we offer five significant examples, including generalized Cattaneo diffusion wave and Klein–Gordon equations. The implementation of algorithms exposes high accuracy of the presented numerical method. The advantage of having compact support and orthogonality of these family of wavelets trigger having sparse operational matrices, which reduces the computational time and CPU requirements.


Introduction
Fractional partial differential equations have appeared in variety of applications [1][2][3][4][5][6]. Diffusion wave equations, as an important branch of these equations, have attracted considerable attention amongst researchers due to their extensive applications in science and engineering. These equations have been widely utilized in some significant physical materials, such as percolation clusters, amorphous, colloid, glassy and porous materials through fractals, dielectrics, semiconductors, polymers and biological systems [7]. In this paper, we aim to solve the general form of Time Fractional Diffusion Wave Equation (TFDWE), which is defined as: subject to the initial and boundary conditions and u(0, t) = g 0 (t), u(L, t) = g 1 (t), (3) where T, L, λ and µ are positive real constants, the fractional derivative, ∂ ξ ∂t ξ , is considered in Caputo sense for ξ (1 < ξ ≤ 2) and the positive constant θ depends on the temperature, the friction coefficient, the universal gas constant and finally on the Avagadro number. The functions f 0 , f 1 , g 0 , g 1 and F are given and represent sufficiently smooth functions. The unknown function u represents the probability density of finding a particle at x place in time t.
The TFDWE defined in (1) is given in the general form that includes particular cases: (I) For ϑ = 1, µ = 0 and F = 0, Equation (1)  for ϑ = 0. We observe that linear and nonlinear Klein-Gordon equations model many problems in classical and quantum mechanics, solitons and condensed matter physics. For example, nonlinear sine Klein-Gordon equation models a Josephson junction, the motion of rigid pendula attached to a stretched wire, dislocations in crystals and DNA dynamics [11][12][13][14].
Since the aforementioned fractional problems does not possess analytical solution, the numerical methods gain much importance in obtaining the solutions to these models. For solving TFDWE, a spectral Tau algorithm based on Jacobi operational matrix is presented in [15] and Legendre wavelets and their fractional operational matrix of integration have been applied in [16]. Authors of [17] employed Euler wavelets and matrix inversion of M × M collocation points via Gaussian quadrature to reduce the fractional problem into an algebraic system.
Haar wavelets collocation and finite difference methods are applied for solving (1 + 1) and (1 + 2) dimensional TFDWEs in [18]. The solutions of multi-term time fractional wave diffusion equation by using transform based local meshless method and quadrature is discussed in [19]. The Spectral collocation method and the finite difference scheme based on cubic trigonometric B-spline were introduced in [20,21], respectively. On the other hand, the family of wavelets are employed extensively for solving the fractional order of differential and integro-differential equations [22][23][24].
It has been proven that the wavelets family are appropriate tools for time-frequency analysis and specially for approximating the signals (time-variable functions), and thus many numerical approaches based on the wavelets have been purposed in the literature for solving a variety of problems; such as [25][26][27][28][29][30]. However, under the extension of research objects and scope, the wavelets theory was discovered to have shortcomings. For example, some signals, such as chirp-like signals [31] (which are ubiquitous in nature) and man-made systems, cannot be processed by classic wavelet theory.
In order to fix wavelet shortcomings, some novel tools based on wavelets have been introduced by researchers, such as fractional wavelet transform [32], a generalization of the wavelet transform. On the other hand, in the computational methods, the operational matrices of the fractional order derivative have attracted the attention of researchers for solving fractional order equations. The determination of the operational matrices of fractional order derivatives has computational complexity by non-fractional basis rather than the fractional basis.
To make a comparison, suppose we have two fractional B f and non-fractional B bases, and t 2 ∈ B and t 3/2 ∈ B f . It is clear that the representing D 1/2 {t 3/2 } in B f is simpler than showing D 1/2 {t 2 } in B. Indeed, we should approximate the operational matrix of derivative on a non-fractional basis, but on a fractional basis, we can evaluate the operational matrix exactly. Therefore, the fractional basis and fractional wavelets family came to the attention of researchers in the numerically solving fractional order differential equations.
Motivated by this, in the present study, the fractional Bernoulli wavelets, which are constructed via operational matrix of fractional Caputo derivative are applied for solving Equations (1)-(3). These fractional wavelets are constructed by Bernoulli polynomials [33][34][35] and have significant properties, such as having compact support and orthogonality. The important advantage of the fractional Bernoulli wavelets is having exact operational matrices, which improves the accuracy of the employed method.
We organize the rest of the paper as follows: Section 2 introduces some essential preliminaries and definitions. The definition of Bernoulli wavelets, fractional Bernoulli wavelets, function approximation , and operational matrix of derivatives are given in Section 3. Section 4 is devoted to the numerical implementation. Section 5 exposes the accuracy and efficiency of the presented method as well as presents some illustrative examples. Our conclusions regarding the proposed method are discussed in Section 6.

Preliminaries on Fractional Calculus
In this section, we present some basic definitions and concepts on fractional calculus that are essential for th subsequent discussion. There are various definitions for fractional integration and derivatives. However, we use fractional Riemann-Liouville integration and fractional Caputo derivatives due to their advantages over other fractional operators. Definition 1. The Riemann-Liouville fractional integral operator of nonnegative order ι is defined as [36] where we set J 0 f (x) = f (x). Some of these are operator features as under and, for the real constant λ, we have Definition 2. The Caputo fractional derivative operator of nonnegative order ι is defined as [36] For the Caputo derivative, we have [37] D ι x ϑ = 0, ϑ ∈ N 0 , ϑ < ι , and where ι denotes the smallest integer greater than or equal to ι, and ι denotes the largest integer less than or equal to ι and N 0 = {0, 1, 2, . . . } and R is the set of real numbers. Similar to the Riemann-Liouville fractional integral operator, the Caputo fractional derivative is linear-that is, for the real constant λ, we have The Laplace transform of Caputo fractional derivative is as follows The relations between the Reimann-Liouville fractional integral and Caputo fractional derivative operators can be addressed by the following identities [38]: and

Main Results
In this section, definitions of Bernoulli wavelets (BWs), fractional Bernoulli wavelets (FBWs), operational matrix of fractional Caputo derivative, and function approximation by these wavelets are presented.

BWs and FBWs
BWs of order m, ψ nm (t) = ψ(k,n, m, t) consist of four arguments,n = n − 1, n = 1, 2, . . . , 2 k−1 , k is a positive integer and t is the normalized time. These wavelets on the interval [0, 1) are defined by [39] whereB 0 (t) = 1 andB and Λ m = ϑ 2m is the normality coefficient. B m , m = 0, 1, . . . , M − 1 are known as Bernoulli polynomials and defined by where ϑ j := B j (0) are Bernoulli numbers. Therefore, Bernoulli wavelets for m > 0 can be rewritten as where and constructed by change of variable t to x ι , (ι > 0) on the BWs [39], that is,  Remark 1. The Bernoulli polynomials are not orthogonal, and they are satisfied in the following relation [39]: Thus, the FBWs, which are constructed by Bernoulli polynomials, are not orthogonal either.

Function Approximation by FBWs
A function f ∈ L 2 [0, 1] can be approximated by FBWs as If the infinite series in equation (17) is truncated in some suitable m and n, we find where C M and Ψ ι k,M are 2 k−1 × M-dimensional column vectors and defined as In order to determine the coefficients in (18), we place and Now, substituting (18) in (21), we find and thus placing The two variable function v(x, t) could be approximated by two dimensional FBWs as where It is clear that for

Operational Matrix of Riemann-Liouville Fractional Integration for FBWs
The Riemann-Liouville fractional integration of Ψ ι can be obtained as where F ξ,ι is the relative operational square matrix of dimension 2 k−1 × M and could be evaluated by using Equations (5) and (15) as follows On the other hand, Thus, using Equation (28) and (29), we can write where Now, we expand x ιi+ξ in terms of FBWs, as Equation (18) x ιi+ξ by (30) and (31), we find where , and 0 is a 1 × M-dimensional row matrix where all its entries are zero.

Operational Matrix of Derivative for FBWs
The derivative of Ψ ι can be obtained as where D is a relative operational square matrix of dimension 2 k−1 × M and can be evaluated as follows On the other hand, Therefore, by using Equations (34) and (35), we can write where Now, we expand x ιi−1 in terms of FBWs, as Equation (18) x ιi−1 by (36) and (37), we find Therefore, we have where D l = µ l,0,0,0 , µ l,1,0,0 , . . . , µ l,M−1,j,i 1×M .

Numerical Implementation
By using the properties of FBWs and their operational matrices of derivative and Riemann-Liouville fractional integration and spectral collocation and Galerkin methods, a new approach is introduced in this section for solving general form of TFDWEs. For solving this kind of equations, we first expand ∂t ξ by FBWs of order ι as where U is 2 k−1 × M-dimensional square unknown coefficients matrix. By using Riemann-Liouville fractional integral of order ξ with respect to variable t, we find by Equation (40) and derivative operational matrix for FBWs, we find and By substituting Equations (39)- (42) in Equation (1), we have where For solving the current equation first, we apply the Galerkin method via FBWs of order ι, that is, the operator 1 0 (.)Ψ ι (t)t ι−1 dt is utilized for Equation (43). Thus, we find where λ n,m n,m , and λ n,m n,m was defined in (22). Therefore, we obtain an algebraic system of equations with 2 k−1 × M equation and (2 k−1 × M) × (2 k−1 × M) unknowns. For solving the obtained system, the collocation method is employed. First, we define residual error function as: and now we consider the collocation points as On substituting the collocation points (in (45)), we find Now, we use thus, we obtain an algebraic linear system with 2 k−1 × M − 2 equations and 2 k−1 × M unknowns. Further, using the boundary conditions (3) we have by collocating the modified boundary conditions in points we also find 2 × 2 k−1 × M equations. Thus, an algebraic linear system is obtained with the same equations and unknowns, which can be easily solved. In this study, the achieved system was solved using the symbolic calculus software Wolframe Mathematica 10. By specifying the unknown matrix U, we find the approximate solution of TFDWE as (40).

Illustrative Examples
In this section, some illustrative examples are presented for the sake of demonstrating the accuracy and efficiency of the proposed method. Let υ and υ * be the exact and approximate solutions of the problems, respectively. Symbols used in the tables are defined as under where T is the final time of problem. We start by the following example where the fractional Bernoulli wavelets are defined for k = 2.

Example 1 ([40]). Consider the following TFDWE
subject to the initial and boundary conditions where for 1 < ξ ≤ 2 the exact solution of this problem is υ(x, t) = (t 2 − t) sin(πx). This problem is solved by the presented method in Section 4 for M = 2, 3, 4 and ι = 0.5, 1.5, 2. Figures 2 and 3 show the space-time curves of exact and approximate solutions and error curves of the approximate solutions of Example 1 for ξ = 1.5 and T = 2, respectively for M = 3 and M = 4. In Table 1, the maximum absolute error ∞ obtained by the purposed method for T = 0.2, 0.4, N = 20 and ξ = 1.5, 1.7 are reported at the final time T and compared with the result of [40], which applied the fractional finite difference method based on the Hermite formula, and the result [21], which used cubic trigonometric basis functions for solving Example 1.
For analysing the effect of M and ι on the numerical solution, the example (50) has been considered and solved for different values of M and ι. According to the reported results in Table 1, we can see that the maximum absolute errors were decreased by increasing the value of M-in other words, when we used a large basis of Bernoulli wavelets, we achieved more small absolute errors. This is valid for the fractional wavelets order ι as well-that is, for the larger ι, we found smaller absolute errors.
where the exact solution of this problem is υ(x, t) = t(t ξ+2 + 1)e x . This problem is solved by the presented method in Section 4 for M = 2, 3, 4 and ι = 0.75, 1.25, 2. Figure 4 shows the space-time curves of exact and approximate solutions of Example 2 for M = 3, ι = 1.5 and ξ = 1.5, and in Figure 5 the error curves of the approximate solutions of Example 2 for M = 4, ι = 1 and ι = 1.5 are shown. In Table 2, the maximum absolute error ∞ , for N = 20 and different values of M, ι and ξ = 1.3, 1.5, 1.8 are reported at the final time T = 1 and compared with the similar results of [15,41]. The authors of [15] applied the shifted Jacobi tau spectral procedure together with the Jacobi operational matrix for solving Example 2, where N and M in the second column of Table 2 are the degrees of shifted Jacobi polynomials in expanding the unknown two-variable function υ(x, t) with respect to the variables t and x, respectively. In [41], the finite difference scheme is presented for finding the solution of Example 2.

Example 3. Consider the following TFKGE [42]
the exact solution of this problem is υ(x, t) = t 2 (e − e x ) sin(x). This problem is solved by the presented method in Section 4 for M = 3, 4 and different values of ι. Figure 6 shows  [42,43].
In [42], the authors used a high order difference method solved the Example 3 where τ = ∆t and h = ∆x. In [43], two-dimensional Genocchi polynomials and the Ritz-Galerkin method via a satisfier function were developed to investigate the numerical solution of Example 3. The parameter N in the second column of Table 3 investigates the degree of Genocchi polynomials, employed for expanding the unknown function in [43].
The exact solution of this problem is υ(x, t) = t 2 x sin(x − 1). Numerical results for this example are reported in Table 5 and Figure 8.

Example 5. Consider the TFGCE [42]
with For obtaining the solution of Equations (53) and (54), we use the Laplace transform and Riemann sum approximation for Laplace inversion. Let V(x, s) be the Laplace transform of υ(x, t) with respect to the variable t. Then, the Laplace transform of this problem is sV(x, s) + s 1+ξ V(x, s) = ∂ 2 ∂x 2 V(x, s), By solving this second order boundary value problem (with respect to variable x), we find Using the Riemann sum approximation for Laplace inversion, we find For obtaining the boundary conditions (54), we choose N = 500 in (55). The numerical results for this example are reported in Table 6 and Figures 9 and 10. Further, the approximated solutions for different values of ξ at t = 0.4 are shown in Figure 11.

Conclusions
In this paper, the fractional Bernoulli wavelets were defined in new settings and applied based on the wavelet spectral methods for solving an important family of time fractional partial differential equations, including the generalized Cattaneo diffusion wave and Klein-Gordon equations. First, the operational matrices of Caputo fractional and ordinary derivatives were constructed and then employed for reducing the time fractional partial equation to an algebraic system. For solving the proposed system, wavelet Galerkin and collocation methods were used. In the following sentences, we briefly describe the novelty of the study and the significance of the results.
Most presented spectral methods in the literature have used classic non-fractional base functions or wavelets, such as Chebyshev, Legendre, and Splines. We applied fractional Bernoulli wavelets as basic functions. The operational matrices construction for fractional order derivatives has computational complexity by a non-fractional basis rather than a fractional basis. Indeed, for most of the bases, the operational matrices of the fractional derivatives in a non-fractional basis have been approximated; however, in fractional bases, we can evaluate the operational matrix exactly. The fractional Bernoulli wavelets have been introduced by some researchers already, but the purposed scheme for the construction of the operational matrices is more applicable and simpler. In this study, the operational matrices were constructed upon the direct definitions of wavelets and Riemann-Liouville fractional integration operator.
Most presented spectral techniques for the solving the obtained system have been applied double collocation method (with respect to time and local variables, t and x), but in the introduced method for reducing the computational cost, after substituting the spectral solution in main fractional PDE, the Galerkin approach was performed with respect to the time variable on the achieved equation. Thus, the two-variable residual error function was reduced to one variable, and the collocation method was applied by the defined meshes for variable x.
From the tabulated maximum absolute errors in the tables, the results obtained by other methods are easily attained by the presented method for a low number of collocation points and, consequently, have less computational cost. According to the exact solutions and their achieved approximations for the illustrative examples, the importance and efficiency of the suggested methods is clearly visible in the figures and tables. The proposed approach has promising applications as it can be extended and applied to the study of the numerical solutions of some other types of fractional integro-differential equations.