Wavelet based algorithm for numerical study of (1+2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(1+2)$\end{document}-dimensional time fractional diffusion problems

An effective and robust scheme is developed for solutions of two-dimensional time fractional heat flow problems. The proposed scheme is based on two-dimensional Haar wavelets coupled with finite differences. The time fractional derivative is approximated by an L1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$L_{1}$\end{document}-formula while spatial part is approximated by two-dimensional Haar wavelets. The proposed methodology first converts the problem to a discrete form and then with collocation approach to a system of linear equations which is easily solvable. To check the efficiency of the scheme, two error norms, E∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$E_{\infty }$\end{document} an Erms\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$E_{\mathrm{rms}}$\end{document}, have been computed. The stability of the scheme has been discussed which is an important part of the manuscript. It is also observed that the spectral radius of the amplification matrix satisfies a stability condition. From computation it is clear that computed results are comparable with the exact solution.


Introduction
Fractional differential equations (FDEs) play a crucial role in modeling of different physical phenomena such as heat and mass transfer [1,2], fluid mechanics [3], and financial theory [4,5]. Usually these equations arise in the form of fractional partial differential equations (FPDEs). The well known application of FPDEs is their nonlocal property which means that the next state of the dynamical system not only depends on the present state but also on all previous states. That's the reason why fractional calculus has been given more emphasis in the last few decades. One of the most important models of FPDEs is time fractional diffusion equation (TFDE) which is used to describe anomalous diffusion phenomena in transport processes. This model can be obtained from the classical diffusion model by replacing the first order time derivative with a fractional one. Analytical methods have limited ability to solve these problems, therefore, numerical techniques are essential.
In the literature several approaches have been proposed for the solution of TFDEs, for example, finite difference methods [6,7], spectral method [8], finite element method [9], etc. In contrast to finite difference methods, compact finite difference methods have good accuracy [10,11], which have been applied for the solution of 1D TFDEs in [12]. For 2D problems, explicit and implicit central difference formulas coupled with extrapolation and kernel based approximation with adaptive technique have been studied in [13,14]. For multidimensional problems, another method is operator-splitting technique or ADI which transforms a given problem to progressions of 1D linear problems having low computational cost. For complete details of ADI, the interested readers are referred to the papers [15,16]. The idea of the ADI method was extended to FPDEs [17,18].
Among diverse numerical methods, wavelet based methods achieved much attention for solution of FPDEs recently because of ease implementation. Also an interested feature of wavelet based methods is that they are capable of noticing singularities and approximating a function with different resolution levels. There are several families of wavelets, like Daubechies, Symlet, Coiflet, etc. All of them are equally interesting, but Daubechies wavelets of order one, known as Haar wavelets (HW), deserve special attention. HW are composed of piecewise constant functions which attain the values -1 and 1. These wavelets are good for approximations but have the drawback of discontinuity at the endpoints of an interval. This problem has been discussed by Cattani [19,20] who used interpolating splines to regularize these wavelets. Another approach was used by Chen and Hasio [21], in their work the highest order derivative has been approximated by truncated HW series instead of solution. Further, this approach was used by Lepik [22,23] to study different problems.
The above discussion about Haar wavelets is limited to integer order differential equations, however, this approach also works for solution of fractional differential equations. Lepik [24] solved fractional integral equations with the help of Haar wavelets. Li [8] applied Haar wavelets for solution of fractional ordinary differential equations (FODEs). Chen et al. [25] analyzed errors of the numerical method for FODEs by Haar wavelets. Ray and Patra [26] solved nonlinear oscillatory van der Pol system using Haar wavelet operational matrices. Yi and Huang [27] implemented an operational matrix method for variable coefficient fractional differential equations. Saeed et al. [28] used Haar wavelet Picard technique for numerical solution of fractional order initial and boundary value problems. Saeed and Rehman [29] applied the same technique for numerical solution of fractional partial differential equations. Majak et al. [30] studied higher order Haar wavelet method for FGM structures. Recently, Alderremya et al. [31] implemented the spectral collocation method together with third order Chebyshev polynomials to certain new models of multispace-fractional Gardner equation. Zhang et al. [32] computed the general solution for impulsive differential equations with Riemann-Liouville fractional-order q ∈ (1, 2). Agarwal and Singh [33] studied mathematical model of Nipah virus with fractional order approach. Morales-Delgado et al. [34] discussed fractional dynamics of oxygen diffusion through capillary to tissues under the influence of external forces considering the fractional operators of Liouville-Caputo and Caputo-Fabrizio. Choi and Agarwal [35] established a fractional integration formula using generalized multiindex Mittag-Leffler function. Agarwal et al. [36,37] studied fixed point results, some new differential equation formulas for extended Mittag-Leffler-type functions and an extension of the Caputo fractional derivative operator involving generalized hypergeometric-type functions. Amiri et al. [38] developed a numerical method based on cosine-trigonometric basis functions to solve Fredholm integral equations of the second kind. Moghadam et al. [39] applied Bernoulli wavelet method for numerical solutions of advection dispersion equations. Farooq et al. [40] used Chebyshev wavelet method for solution of delay differential equations. Khalil et al. [41] introduced some new operational matrices of integration for fractional Poisson equations with integral type boundary conditions.

Motivation
Developing methods for solutions of TFDEs either analytical, semianalytical, or numerical is an important task. The motivation of this work is to compute numerical solutions of 2D TFDEs using two-dimensional Haar wavelets and finite differences. The proposed scheme is hybrid and gives good results which are comparable with exact solutions. Also we will examine the computational stability of the proposed algorithm which is related to the spectral radius of the amplification matrix. The problem which will be consider in this study is defined as follows: together with the initial condition and boundary conditions where X = X(x, y), Ψ and ∂Ψ denote solution domain and boundary, B(X, t), β(X), η(X, t) are known functions, while C(X, t) is an unknown function to be computed. In Eq. (1.1), is two-dimensional Laplacian and ∂ γ C(X,t) ∂t γ is the fractional derivative of C with respect to t in Caputo sense which is defined as: The rest of the paper is organized as follows: Basic definitions of HW and its integrals are reported in Sect. 2. The proposed methodology and convergence with stability analysis have been discussed in Sects. 3 and 4, respectively. Numerical results are given in Sect. 5, and finally conclusion is summarized in Sect. 6.

Convergence and stability
First we address the convergence theorem of the proposed scheme. The following lemma is necessary for convergence theorem. The above result indicates that the resolution level and error norm are inversely proportional. It has been verified that when resolution increases the error norm decreases.

Lemma 1 ([42]) If f (x, y) satisfies a Lipschitz condition on
Proof To prove the theorem, take the asymptotic expansion of Eq. (3.10) defined as follows: From Eq. (4.2) at resolution level J, we can write The L 2 -norm is given by To evaluate the complex integral, it is sufficient to estimate the maximum bounds of Haar wavelets integrals which can be obtained from the expression [30]: From Eq. (4.5), one can deduce that (4.7) Using Eq. (4.7) and the above lemma in Eq. (4.4), we can derive which gives This completes the proof of the theorem.

Illustrative test problems
In this section of the manuscript, we solve some test problems to validate the proposed method. To check goodness, we compute the error norms, E ∞ and E rms , which are defined as follows:   3)) have been extracted from the exact solution C(x, y, t) = t 2 sin(x) sin(y). The problem has been solved for τ = 0.01, J = 4, γ = 0.1, 0.3, 0.5, 0.7, 0.9 and the obtained results tabulated in Table 1. In Table 2 calculated error norms with decreasing time step size have been given, which show that the scheme is convergent in time. In Tables 1 and 2, the spectral radius of the amplification matrix has been mentioned. It is obvious from the tables that computed solutions agree with exact, and the spectral radius lies in the stability domain. Solution profile and absolute errors are plotted in Fig. 1. From the figure it is clear that the proposed solutions agree well with the exact solution.
Example 5.2 Consider Eq. (1.1) with the exact solution C(x, y, t) = exp(x + y)t 1+γ . Initial and boundary conditions have been used from the exact solution. The source term B(x, y, t) is easy to derive corresponding to the exact solution. The problem has been solved for different times, and the obtained results have been reported in Table 3. In this problem, we fix t = 0.2, γ = 0.9, while varying resolution level J. The achieved results are presented in Table 4. One can see clearly that when the resolution level increases the accuracy also increases which guarantees spatial convergence. The same table reveals that the spectral radius remains the same when the resolution level varies. To clarify the spatial convergence    more, the absolute errors E abs at different mesh points have been addressed in Table 5.
Graphical solutions together with absolute errors are shown in Fig. 2. From the figure it is obvious that the obtained results and exact solutions have strong agreement.
Example 5.3 Consider Eq. (1.1) with the exact solution C(x, y, t) = t 2 (xx 2 ) 2 (yy 2 ) 2 . The source term is easy to compute from the exact solution. In Table 6, we recorded the computed error norms for parameters t = 0.5, 0.1, γ = 0.2, 0.4, 0.6, 0.8, τ = 0.01. In Fig. 3, exact versus approximate solutions and absolute errors have been plotted. From tabulated data and graphical solution, it is clear that computed solutions match well with the exact one.
where B(X, t) = sin(x) sin(y) 4 The associated initial and boundary conditions have been extracted from the exact solution C(x, y, t) = sin(x) sin(y)(t 1.4 + t 1.2 + t + 1) 4 .
In Table 7, we presented computed error norms for the nonlinear problem at different values of γ . From the table it is obvious that the current scheme also works for twodimensional nonlinear time fractional problems. Graphical solutions together with absolute errors have been shown Fig. 4. From the figure one can observe that computed solutions match well with the exact solution.  Example 5.5 Finally, we consider a nonlinear time fractional differential equation of the form ∂ γ C(X, t) ∂t γ = C(X, t) + C(X, t) -C 2 (X, t) + B(X, t), t > 0, X ∈ Ψ , t γ + t 2 x 2 (1x) 2 y 2 (1y) 2 1t γ + t 2 x 2 (1x) 2 y 2 (1y) 2 .
The exact solution of this problem is C(x, y, t) = (t γ + t 2 )x 2 (1x) 2 y 2 (1y) 2 . It is obvious that the first derivative of the exact solution approaches infinity when t → 0, which shows that there is an initial singularity in time. Computed results have been addressed in Table 8 for J = 3. From the table it is clear that the suggested scheme works if an initial singularity exists in time. Also the same table shows time convergence of the suggested scheme for a singularity problem.

Concluding remarks
In this work, we proposed a numerical scheme based on two-dimensional Haar wavelets combined with finite differences for solution of time fractional diffusion equations. The scheme has been applied for solutions of five test problems and noted the achieved results in tabulated as well as in graphical form. It has been observed that two-dimensional Haar wavelets work well for solution of two-dimensional time fractional heat problems.