**1. 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–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 + 1)-Dimensional Equation:*

$$\mathbf{u}^{\varepsilon}D\_{t}^{\delta}w(\mathbf{x},t) = -w\_{l}(\mathbf{x},t) + w\_{xx}(\mathbf{x},t) + \mathcal{A}(\mathbf{x},t), \quad \mathbf{x} \in \Omega, \quad t \in [0,T], \quad 1 < \delta \le 2,\tag{1}$$

$$\begin{cases} w(\mathbf{x},0) = f(\mathbf{x}), \quad w\_t(\mathbf{x},0) = \mathbf{g}(\mathbf{x}) \quad \mathbf{x} \in \tilde{\Omega} = \Omega \cup \partial\Omega, \\ w(\mathbf{x},t) = u(t), \quad \mathbf{x} \in \partial\Omega \qquad t \in [0,T]. \end{cases} \tag{2}$$

*(1 + 2)-Dimensional Equation:*

$$\mathbf{u}^{\varepsilon} \cdot \mathbf{D}\_t^{\delta} w(\mathbf{x}, y, t) = \Delta w(\mathbf{x}, y, t) + \mathcal{B}(\mathbf{x}, y, t), \quad (\mathbf{x}, y) \in \Phi, \quad t \in [0, T], \quad 1 < \delta \le 2,\tag{3}$$

$$\begin{cases} w(\mathbf{x}, y, 0) = \chi(\mathbf{x}, y), \quad w\_t(\mathbf{x}, y, 0) = \kappa(\mathbf{x}, y), \quad (\mathbf{x}, y) \in \breve{\Phi} = \Phi \cup \partial \Phi, \\ w(\mathbf{x}, y, t) = \chi\_1(\mathbf{x}, y, t), \quad (\mathbf{x}, y) \in \partial \Phi, \quad t \in [0, T]. \end{cases} \tag{4}$$

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 *cDδt w* denotes the time fractional derivative of *w* with respect to *t* in the Caputo sense which is given by

$${}^{c}D\_{t}^{\delta}w = \begin{cases} \frac{1}{\Gamma(2-\delta)} \int\_{0}^{t} \frac{w\_{\zeta\zeta}(\mathbf{x},\tilde{\zeta})}{(t-\zeta)^{\delta-1}} d\zeta, & 1 < \delta < 2, \\\frac{\partial^{2}w(\mathbf{x},t)}{\partial t^{2}}, & \delta = 2. \end{cases} \tag{5}$$

#### **2. 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 2*M* 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:

$$\mathcal{H}\_1(\mathbf{x}) = \begin{cases} 1, & \mathbf{x} \in [a, b] \\ 0, & \text{otherwise.} \end{cases} \tag{6}$$

$$\mathcal{H}\_i(\mathbf{x}) = \begin{cases} 1, & \mathbf{x} \in [\mathfrak{z}\_1(i), \mathfrak{z}\_2(i)) \\ -1, & \mathbf{x} \in [\mathfrak{z}\_2(i), \mathfrak{z}\_3(i)) \\ 0, & \text{otherwise.} \end{cases} \tag{7}$$

where

$$\mathbb{P}\_1\mathbb{P}\_1(i) = a + 2k\nu\delta\mathbf{x},\ \mathbb{P}\_2(i) = a + (2k+1)\nu\delta\mathbf{x},\ \mathbb{P}\_3(i) = a + 2(k+1)\nu\delta\mathbf{x},\ \nu = \frac{M}{m},\ \delta\mathbf{x} = \frac{b-a}{2M}.$$

To solve nth order time fractional PDEs the following repeated integrals are needed:

$$\mathcal{P}\_{i,\emptyset}(\mathbf{x}) = \int\_{a}^{\mathbf{x}} \int\_{a}^{\mathbf{x}} \cdots \int\_{a}^{\mathbf{x}} \mathcal{H}\_{i}(z) dz^{\emptyset} = \frac{1}{(\beta - 1)!} \int\_{a}^{\mathbf{x}} (\mathbf{x} - z)^{\beta - 1} \mathcal{H}\_{i}(z) dz,\tag{8}$$

where

$$
\beta = 1, 2, \dots, n, \quad i = 1, 2, \dots, 2M.
$$

Keeping in view Equations (6) and (7) the close form expressions of these integrals are given by

$$\mathcal{P}\_{1,\beta}(\mathbf{x}) = \frac{(\mathbf{x} - \mathbf{a})^{\beta}}{\beta!}. \tag{9}$$

$$\mathcal{P}\_{i,\beta}(\mathbf{x}) = \begin{cases} 0, & \mathbf{x} < \mathbb{S}\_{1}(\mathbf{x}) \\ \frac{1}{\mathbb{S}!} \left[ \mathbf{x} - \mathbb{S}\_{1}(i) \right]^{\beta} & \mathbf{x} \in \left[ \mathbb{S}\_{1}(i), \mathbb{S}\_{2}(i) \right) \\ \frac{1}{\mathbb{S}!} \left[ (\mathbf{x} - \mathbb{S}\_{1}(i))^{\beta} - 2((\mathbf{x} - \mathbb{S}\_{2}(i))^{\beta}) \right] & \mathbf{x} \in \left[ \mathbb{S}\_{2}(i), \mathbb{S}\_{3}(i) \right) \\ \frac{1}{\mathbb{S}!} \left[ (\mathbf{x} - \mathbb{S}\_{1}(i))^{\beta} - 2((\mathbf{x} - \mathbb{S}\_{2}(i))^{\beta} + (\mathbf{x} - \mathbb{S}\_{3}(i))^{\beta}) \right] & \mathbf{x} \ge \mathbf{x}\_{3}(i). \end{cases} \tag{10}$$

#### **3. 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]

$$\begin{split} \,^c\mathrm{D}\_t^\delta w(x, t^{j+1}) &= \frac{1}{\Gamma(2-\delta)} \int\_0^{j+1} w^2(\xi) \left(x, \xi\right) \left(t^{j+1} - \xi\right)^{1-\delta} d\xi \\ &= \frac{1}{\Gamma(2-\delta)} \sum\_{k=0}^j \int\_{i^j}^{j+1} \left[\frac{w^{k+1} - 2w^k + w^{k-1}}{\tau^2}\right] \left(t^{j+1} - \xi\right)^{1-\delta} d\xi \\ &= \frac{1}{\Gamma(2-\delta)} \sum\_{k=0}^j \left[\frac{w^{k+1} - 2w^k + w^{k-1}}{\tau^2}\right] \int\_{i^j}^{j+1} \left[(j+1)\tau - \xi\right]^{1-\delta} d\xi \\ &= \frac{1}{\Gamma(2-\delta)} \sum\_{k=0}^j \left[\frac{w^{k+1} - 2w^k + w^{k-1}}{\tau^2}\right] \frac{(j-k+1)^{2-\delta} - (j-k)^{2-\delta}}{(2-\delta)(\tau^{\delta-2})} \\ &= \frac{\tau^{-\delta}}{\Gamma(3-\delta)} \sum\_{k=0}^j \left[w^{j-k+1} - 2w^{j-k} + w^{j-k-1}\right] \left[(k+1)^{2-\delta} - (k)^{2-\delta}\right] \\ &= A\_\delta \left[w^{j+1} - 2w^j + w^{j-1}\right] + A\_\delta \sum\_{k=1}^j \left[w^{j-k+1} - 2w^{j-k} + w^{j-k-1}\right] B(k), \end{split} \tag{11}$$

where *Aδ* = *τ* −*δ* <sup>Γ</sup>(<sup>3</sup>−*<sup>δ</sup>*), *τ* is time step size and *B*(*k*) = (*k* + 1)<sup>2</sup>−*<sup>δ</sup>* − (*k*)<sup>2</sup>−*<sup>δ</sup> Casei:*

Using Equation (11) and *θ*−weighted scheme (0 ≤ *θ* ≤ 1) in Equation (1), we obtain

$$\begin{split} A\_{\delta} \left[ w^{j+1} - 2w^j + w^{j-1} \right] + A\_{\delta} \sum\_{k=1}^{j} \left[ w^{j-k+1} - 2w^{j-k} + w^{j-k-1} \right] B(k) + \frac{1}{\tau} \left\{ w^{j+1} - w^j \right\} \\ = \theta w\_{\text{xx}}^{j+1} + (1 - \theta) w\_{\text{xx}}^{j} + \mathcal{A}(\mathbf{x}, t^{j+1}). \end{split} \tag{12}$$

After simplification, the above equation transforms to

$$\begin{aligned} (\tau A\_\delta + 1)w^{j+1} - \tau \theta w\_{xx}^{j+1} &= 2\tau A\_\delta w^j - \tau A\_\delta w^{j-1} - \tau A\_\delta \sum\_{k=1}^j \left[ w^{j-k+1} - 2w^{j-k} + w^{j-k-1} \right] B(k) \\ &+ w^j + \tau (1 - \theta) w\_{xx}^j + \tau \mathcal{A}(x, t^{j+1}). \end{aligned} \tag{13}$$

In our analysis we take *θ* = 1/2. Now approximating the highest order derivative by a truncated Haar wavelets series as:

$$w\_{xx}^{j+1}(\mathbf{x}) = \sum\_{i=1}^{2M} a\_i^{j+1} \mathcal{H}\_i(\mathbf{x}). \tag{14}$$

.

Integrating Equation (14) from 0 to *x*

$$w\_x^{j+1}(\mathbf{x}) = \sum\_{i=1}^{2M} a\_i^{j+1} \mathcal{P}\_{i,1}(\mathbf{x}) + w\_x^{j+1}(\mathbf{0}).\tag{15}$$

Integrating Equation (15) from 0 to 1, we ge<sup>t</sup>

$$w\_x^{j+1}(0) = w^{j+1}(1) - w^{j+1}(0) - \sum\_{i=1}^{2M} a\_i^{j+1} \mathcal{P}\_{i,2}(1). \tag{16}$$

Substituting Equation (16) in Equation (15), the resultant equation reduces to

$$w\_{\mathbf{x}}^{j+1}(\mathbf{x}) = \sum\_{i=1}^{2M} a\_i^{j+1} \left[ \mathcal{P}\_{i,1}(\mathbf{x}) - \mathcal{P}\_{i,2}(1) \right] + w^{j+1}(1) - w^{j+1}(0). \tag{17}$$

Integration of Equation (17) from 0 to *x* yields

$$w^{j+1}(\mathbf{x}) = \sum\_{i=1}^{2M} a\_i^{j+1} \left[ \mathcal{P}\_{i,2}(\mathbf{x}) - \mathbf{x} \mathcal{P}\_{i,2}(\mathbf{1}) \right] + \mathbf{x} \left[ w^{j+1}(\mathbf{1}) - w^{j+1}(\mathbf{0}) \right] + w^{j+1}(\mathbf{0}).\tag{18}$$

Substituting values from Equations (14), (17) and (18) in Equation (13) and using collocation points *xm* = *<sup>m</sup>*−0.5 2*M* , *m* = 1, 2, . . . 2*M*, leads to the following system of algebraic equation

$$\sum\_{i=1}^{2M} a\_i^{j+1} \left[ \left( \tau A\_\delta + 1 \right) \left\{ \mathcal{P}\_{i,2}(\mathbf{x}) - \mathbf{x} \mathcal{P}\_{i,2}(1) \right\} - \tau \theta \mathcal{H}\_i(\mathbf{x}) \right]\_{\mathbf{x} = \mathbf{x}\_\mathbf{m}} = \mathcal{R}(m), \tag{19}$$

where

$$\begin{split} \mathcal{R}(m) &= 2\tau A\_{\delta}w^{j} - \tau A\_{\delta}w^{j-1} - \tau A\_{\delta} \sum\_{k=1}^{j} \left[ w^{j-k+1} - 2w^{j-k} + w^{j-k-1} \right] B(k) + w^{j} \\ &+ \tau (1 - \theta) w\_{xx}^{j} + \tau \mathcal{A}(x\_{m}, t^{j+1}) - (\tau A\_{\delta} + 1) \left\{ x\_{m} \left( w^{j+1}(1) - w^{j+1}(0) \right) + w^{j+1}(0) \right\}. \end{split}$$

Equation (19) contains 2*M* 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).
