Approximate Solutions of Time Fractional Diffusion Wave Models

Abstract

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 finite 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 finite 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_{\infty }$ error norms. Computations validate that the proposed method produces good results, which are comparable with exact solutions and those presented before.

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,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 + 1)-Dimensional Equation:

  $$\begin{array}{c}{}^c{D}_t^{\delta }w(x,t)=-w_t(x,t)+w_{xx}(x,t)+\mathcal{A}(x,t),x\in \mathsf{\Omega },t\in [0,T],1<\delta \le 2,\end{array}$$

  (1)

  $$\begin{array}{c}\left\{\begin{array}{c}w(x,0)=f\left(x\right),w_t(x,0)=g\left(x\right)x\in \tilde{\mathsf{\Omega }}=\mathsf{\Omega }\cup \partial \mathsf{\Omega },\hfill \\ w(x,t)=\alpha \left(t\right),x\in \partial \mathsf{\Omega }t\in [0,T].\hfill \end{array}\right.\end{array}$$

  (2)
• (1 + 2)-Dimensional Equation:

  $$\begin{array}{c}{}^c{D}_t^{\delta }w(x,y,t)=\Delta w(x,y,t)+\mathcal{B}(x,y,t),(x,y)\in \mathsf{\Phi },t\in [0,T],1<\delta \le 2,\end{array}$$

  (3)

  $$\begin{array}{c}\left\{\begin{array}{c}w(x,y,0)=\chi (x,y),w_t(x,y,0)=\kappa (x,y),(x,y)\in \tilde{\mathsf{\Phi }}=\mathsf{\Phi }\cup \partial \mathsf{\Phi },\hfill \\ w(x,y,t)={\chi }_1(x,y,t),(x,y)\in \partial \mathsf{\Phi },t\in [0,T].\hfill \end{array}\right.\end{array}$$

  (4)

In Equations (1)–(4), $\Delta$ is two-dimensional Laplacian; $\mathcal{A},\mathcal{B},f,g,\alpha ,\chi ,\kappa ,{\chi }_1$ are known functions and w is unknown function. Equations (2) and (4) are the corresponding initial and boundary conditions. The symbols, $\mathsf{\Omega }$ and $\partial \mathsf{\Omega }$, $\mathsf{\Phi }$ and $\partial \mathsf{\Phi }$ represent the domain and boundary of the domain respectively for (1 + 1)- and (1 + 2)-dimensional problems. Also ${}^c{D}_t^{\delta }w$ denotes the time fractional derivative of w with respect to t in the Caputo sense which is given by

$$\begin{array}{c}{}^c{D}_t^{\delta }w=\left\{\begin{array}{c}\frac{1}{\mathsf{\Gamma }\left(2-\delta \right)}{\int }_0^t\frac{w_{\zeta \zeta }(x,\zeta )}{{\left(t-\zeta \right)}^{\delta -1}}d\zeta ,1<\delta <2,\hfill \\ \frac{{\partial }^{2}w(x,t)}{\partial t^2},\delta =2.\hfill \end{array}\right.\hfill \end{array}$$

(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\in [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,.\cdots ,J$ and $k=0,.\cdots ,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:

$$\begin{array}{c}{\mathcal{H}}_1\left(x\right)=\left\{\begin{array}{cc}1,\hfill & x\in \left[a,b\right]\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

(6)

$$\begin{array}{c}{\mathcal{H}}_i\left(x\right)=\left\{\begin{array}{cc}1,\hfill & x\in \left[{\xi }_1\left(i\right),{\xi }_2\left(i\right)\right)\hfill \\ -1,\hfill & x\in \left[{\xi }_2\left(i\right),{\xi }_3\left(i\right)\right)\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\hfill \end{array}$$

(7)

where

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

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

$$\begin{array}{c}\hfill {\mathcal{P}}_{i,\beta }\left(x\right)={\int }_a^x{\int }_a^x\cdots {\int }_a^x{\mathcal{H}}_i\left(z\right)d{z}^{\beta }=\frac{1}{\left(\beta -1\right)!}{\int }_a^x{\left(x-z\right)}^{\beta -1}{\mathcal{H}}_i\left(z\right)dz,\end{array}$$

(8)

where

$$\beta =1,2,\cdots n,i=1,2,\cdots 2M.$$

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

$$\begin{array}{c}{\mathcal{P}}_{1,\beta }\left(x\right)=\frac{{\left(x-a\right)}^{\beta }}{\beta !}.\hfill \end{array}$$

(9)

$$\begin{array}{c}{\mathcal{P}}_{i,\beta }\left(x\right)=\left\{\begin{array}{cc}0,\hfill & x<{\xi }_1\left(x\right)\hfill \\ \frac{1}{\beta !}{\left[x-{\xi }_1\left(i\right)\right)}^{\beta }\hfill & x\in \left[{\xi }_1\left(i\right),{\xi }_2\left(i\right)\right)\hfill \\ \frac{1}{\beta !}\left[{(x-{\xi }_1\left(i\right))}^{\beta }-2({(x-{\xi }_2\left(i\right))}^{\beta }\right)\hfill & x\in [{\xi }_2\left(i\right),{\xi }_3\left(i\right))\hfill \\ \frac{1}{\beta !}\left[{(x-{\xi }_1\left(i\right))}^{\beta }-2({(x-{\xi }_2\left(i\right))}^{\beta }+{(x-{\xi }_3\left(i\right))}^{\beta }\right)\hfill & x\ge {\xi }_3\left(i\right).\hfill \end{array}\right.\hfill \end{array}$$

(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{array}{cc}\hfill {}^c{\mathbb{D}}_t^{\delta }w(x,t^{j+1})& =\frac{1}{\mathsf{\Gamma }\left(2-\delta \right)}{\int }_0^{t^{j+1}}{w}^{\left(2\right)}(x,\zeta ){\left(t^{j+1}-\zeta \right)}^{1-\delta }d\zeta \hfill \\ \hfill & =\frac{1}{\mathsf{\Gamma }\left(2-\delta \right)}\sum _{k=0}^j{\int }_{t^j}^{t^{j+1}}\left[\frac{w^{k+1}-2w^k+w^{k-1}}{{\tau }^2}\right]{\left(t^{j+1}-\zeta \right)}^{1-\delta }d\zeta \hfill \\ \hfill & =\frac{1}{\mathsf{\Gamma }\left(2-\delta \right)}\sum _{k=0}^j\left[\frac{w^{k+1}-2w^k+w^{k-1}}{{\tau }^2}\right]{\int }_{t^j}^{t^{j+1}}{\left[(j+1)\tau -\zeta \right]}^{1-\delta }d\zeta \hfill \\ \hfill & =\frac{1}{\mathsf{\Gamma }\left(2-\delta \right)}\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 )\left({\tau }^{\delta -2}\right)}\hfill \\ \hfill & =\frac{{\tau }^{-\delta }}{\mathsf{\Gamma }\left(3-\delta \right)}\sum _{k=0}^j\left[w^{j-k+1}-2w^{j-k}+w^{j-k-1}\right]\left[{\left(k+1\right)}^{2-\delta }-{\left(k\right)}^{2-\delta }\right]\hfill \\ \hfill & =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\left(k\right),\hfill \end{array}$$

(11)

where $A_{\delta }=\frac{{\tau }^{-\delta }}{\mathsf{\Gamma }\left(3-\delta \right)}$, $\tau$ is time step size and $B\left(k\right)={\left(k+1\right)}^{2-\delta }-{\left(k\right)}^{2-\delta }.$

• Case i:

Using Equation (11) and $\theta -$ weighted scheme $(0\le \theta \le 1)$ in Equation (1), we obtain

$$\begin{array}{ll}{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\left(k\right)+\frac{1}{\tau }\left\{w^{j+1}-w^j\right\}\\ & =\theta w_{xx}^{j+1}+(1-\theta )w_{xx}^j+\mathcal{A}(x,t^{j+1}).\end{array}$$

(12)

After simplification, the above equation transforms to

$$\begin{array}{cc}\hfill (\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\left(k\right)\hfill \\ \hfill & +w^j+\tau (1-\theta )w_{xx}^j+\tau \mathcal{A}(x,t^{j+1}).\hfill \end{array}$$

(13)

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

$$w_{xx}^{j+1}\left(x\right)=\sum _{i=1}^{2M}{a}_i^{j+1}{\mathcal{H}}_i\left(x\right).$$

(14)

Integrating Equation (14) from 0 to x

$$w_x^{j+1}\left(x\right)=\sum _{i=1}^{2M}{a}_i^{j+1}{\mathcal{P}}_{i,1}\left(x\right)+w_x^{j+1}\left(0\right).$$

(15)

Integrating Equation (15) from 0 to 1, we get

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

(16)

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

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

(17)

Integration of Equation (17) from 0 to x yields

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

(18)

Substituting values from Equations (14), (17) and (18) in Equation (13) and using collocation points $x_m=\frac{m-0.5}{2M},m=1,2,\cdots 2M$, leads to the following system of algebraic equation

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

(19)

where

$$\begin{array}{cc}\hfill \mathcal{R}\left(m\right)& =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\left(k\right)+w^j\hfill \\ \hfill & +\tau (1-\theta )w_{xx}^j+\tau \mathcal{A}(x_m,t^{j+1})-\left(\tau A_{\delta }+1\right)\left\{x_m\left(w^{j+1}\left(1\right)-w^{j+1}\left(0\right)\right)+w^{j+1}\left(0\right)\right\}.\hfill \end{array}$$

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

$$\begin{array}{cc}\hfill A_{\delta }{w}^{j+1}-\theta \left[w_{xx}^{j+1}+w_{yy}^{j+1}\right]& =(1-\theta )\left[w_{xx}^j+w_{yy}^j\right]+\mathcal{B}(x,y,t^{j+1})+2A_{\delta }{w}^j-A_{\delta }{w}^{j-1}\hfill \\ \hfill & -A_{\delta }\sum _{k=1}^j\left[w^{j-k+1}-2w^{j-k}+w^{j-k-1}\right]B\left(k\right).\hfill \end{array}$$

(20)

Now we approximate $w_{xxyy}^{j+1}(x,y)$ with a two dimensional truncated Haar wavelets series as:

$$w_{xxyy}^{j+1}(x,y)=\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}{\mathcal{H}}_i\left(x\right){\mathcal{H}}_i\left(y\right),$$

(21)

where $a_{i,l}^{j+1}$ are unknowns to be determined. Integration of Equation (21) w.r.t. to y, between 0 and y, gives

$$w_{xxy}^{j+1}(x,y)=\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}{\mathcal{H}}_i\left(x\right){\mathcal{P}}_{l,1}\left(y\right)+w_{xxy}^{j+1}(x,0).$$

(22)

Integrating Equation (22) w.r.t y from 0 to 1, the unknown term $w_{xxy}^{j+1}(x,0)$ is given by

$$w_{xxy}^{j+1}(x,0)=w_{xx}^{j+1}(x,1)-w_{xx}^{j+1}(x,0)-\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}{\mathcal{H}}_i\left(x\right){\mathcal{P}}_{l,2}\left(1\right).$$

(23)

Substituting Equation (23) in Equation (22), the obtained result is

$$w_{xxy}^{j+1}(x,y)=\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}{\mathcal{H}}_i\left(x\right)\left[{\mathcal{P}}_{l,1}\left(y\right)-{\mathcal{P}}_{l,2}\left(1\right)\right]+w_{xx}^{j+1}(x,1)-w_{xx}^{j+1}(x,0).$$

(24)

Integrating Equation (24) from 0 to y, we get

$$w_{xx}^{j+1}(x,y)=\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}{\mathcal{H}}_i\left(x\right)\left[{\mathcal{P}}_{l,2}\left(y\right)-y{\mathcal{P}}_{l,2}\left(1\right)\right]+y{w}_{xx}^{j+1}(x,1)+(1-y)w_{xx}^{j+1}(x,0).$$

(25)

Repeating the same procedure one can easily derive the subsequent expressions

$$w_{yy}^{j+1}(x,y)=\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}\left[{\mathcal{P}}_{i,2}\left(x\right)-x{\mathcal{P}}_{i,2}\left(1\right)\right]{\mathcal{H}}_l\left(y\right)+x{w}_{yy}^{j+1}(1,y)+\left(1-x\right)w_{yy}^{j+1}(0,y).$$

(26)

$$\begin{array}{cc}\hfill w_x^{j+1}(x,y)& =\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}\left[{\mathcal{P}}_{i,1}\left(x\right)-{\mathcal{P}}_{i,2}\left(1\right)\right]\left[{\mathcal{P}}_{l,2}\left(y\right)-y{\mathcal{P}}_{l,2}\left(1\right)\right]+y{w}_x^{j+1}(x,1)\hfill \\ \hfill & +(1-y)w_x^{j+1}(x,0)+w^{j+1}(1,y)-w^{j+1}(0,y)-y{w}^{j+1}(1,1)+y{w}^{j+1}(0,1)\hfill \\ \hfill & +(y-1)w^{j+1}(1,0)+(1-y)w^{j+1}(0,0).\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill w_y^{j+1}(x,y)& =\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}\left[{\mathcal{P}}_{i,2}\left(x\right)-x{\mathcal{P}}_{i,2}\left(1\right)\right]\left[{\mathcal{P}}_{l,1}\left(y\right)-{\mathcal{P}}_{l,2}\left(1\right)\right]+x{w}_y^{j+1}(1,y)\hfill \\ \hfill & +(1-x)w_y^{j+1}(0,y)+w^{j+1}(x,1)-w^{j+1}(x,0)-x{w}^j(1,1)+x{w}^{j+1}(1,0)\hfill \\ \hfill & +(x-1)w^{j+1}(0,1)+(1-x)w^{j+1}(0,0).\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill w^{j+1}(x,y)& =\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}\left[{\mathcal{P}}_{i,2}\left(x\right)-x{\mathcal{P}}_{i,2}\left(1\right)\right]\left[{\mathcal{P}}_{l,2}\left(y\right)-y{\mathcal{P}}_{l,2}\left(1\right)\right]+y{w}^{j+1}(x,1)\hfill \\ \hfill & -y{w}^{j+1}(0,1)+(1-y)\left[w^{j+1}(x,0)-w^{j+1}(0,0)\right]+x{w}^{j+1}(1,y)\hfill \\ \hfill & -x{w}^{j+1}(0,y)-xy\left[w^{j+1}(1,1)-w^{j+1}(0,1)\right]+x\left(y-1\right)w^{j+1}(1,0)\hfill \\ \hfill & +x(1-y)w^{j+1}(0,0)+w^{j+1}(0,y).\hfill \end{array}$$

(29)

Substitution of Equations (25), (26) and (29) in Equation (20) and using the collocation points, $x_m=\frac{m-0.5}{2M},y_n=\frac{n-0.5}{2M},m,n=1,2,\cdots 2M,$ produces the following system of equations

$$\sum _{i=1}^{2M}\sum _{l=1}^{2M}{a}_{i,l}^{j+1}\left[A_{\delta }\mathcal{D}(i,l,m,n)-\theta \mathcal{E}(i,l,m,n)-\theta \mathcal{F}(i,l,m,n)\right]=\mathcal{L}(m,n)+\mathcal{M}(m,n),$$

(30)

where

$$\begin{array}{cc}\hfill \mathcal{D}(i,l,m,n)& =\left[{\mathcal{P}}_{i,2}\left(x_m\right)-x_m{\mathcal{P}}_{i,2}\left(1\right)\right]\left[{\mathcal{P}}_{l,2}\left(y_n\right)-y_n{\mathcal{P}}_{l,2}\left(1\right)\right],\hfill \\ \hfill \mathcal{E}(i,l,m,n)& ={\mathcal{H}}_i\left(x_m\right)\left[{\mathcal{P}}_{l,2}\left(y_n\right)-y_n{\mathcal{P}}_{l,2}\left(1\right)\right],\hfill \\ \hfill \mathcal{F}(i,l,m,n)& =\left[{\mathcal{P}}_{i,2}\left(x_m\right)-x_m{\mathcal{P}}_{i,2}\left(1\right)\right]{\mathcal{H}}_l\left(y_n\right),\hfill \\ \hfill \mathcal{L}(m,n)& =(1-\theta )\left[w_{xx}^j+w_{yy}^j\right],+\mathcal{B}(x_m,y_n,t^{j+1})+2A_{\delta }{w}^j-A_{\delta }{w}^{j-1}\hfill \\ \hfill & -A_{\delta }\sum _{k=1}^j\left[w^{j-k+1}-2w^{j-k}+w^{j-k-1}\right]B\left(k\right),\hfill \\ \hfill \mathcal{M}(m,n)& =-A_{\delta }[y_n{w}_x^{j+1}(x_m,1)-y_n{w}^{j+1}(0,1)+(1-y_n)\left\{w^{j+1}(x_m,0)-w^{j+1}(0,0)\right\}\hfill \\ \hfill & +x_m{w}^{j+1}(1,y_n)-x_m{w}^{j+1}(0,y_n)-x_m{y}_n\left\{w^{j+1}(1,1)-w^{j+1}(0,1)\right\}\hfill \\ \hfill & +x_m\left(y_n-1\right)w^{j+1}(1,0)+x_m(1-y_n)w^{j+1}(0,0)+w^{j+1}(0,y_n)]+\theta [y_n{w}_{xx}^{j+1}(x_m,1)\hfill \\ \hfill & +(1-y_n)w_{xx}^{j+1}(x_m,0)+x_m{w}_{yy}^{j+1}(1,y_n)+\left(1-x_m\right)w_{yy}^{j+1}(0,y_n)].\hfill \end{array}$$

Equation (30) represents $2M\times 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).

4. 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_{xx}^{j+1}=\mathcal{U}{\alpha }^{j+1}+{\tilde{\mathcal{U}}}^{j+1},$$

(31)

$$w_{yy}^{j+1}=\mathcal{V}{\alpha }^{j+1}+{\tilde{\mathcal{V}}}^{j+1},$$

(32)

$$w^{j+1}=\mathcal{Z}{\alpha }^{j+1}+{\tilde{\mathcal{Z}}}^{j+1},$$

(33)

where ${\alpha }^{j+1}={\alpha }^{j+1}(i,l)$, $\mathcal{U},\mathcal{V},\mathcal{Z}$ and ${\tilde{\mathcal{U}}}^{j+1},{\tilde{\mathcal{V}}}^{j+1},{\tilde{\mathcal{Z}}}^{j+1}$ are interpolation matrices of $w_{xx}^{j+1},w_{yy}^{j+1},w^{j+1}$ at collocation points and boundary terms, respectively. Now using Equations (31), (32) and (33) in Equation (20), we get

$$\left[A_{\delta }\mathcal{Z}-\theta \left(\mathcal{U}+\mathcal{V}\right)\right]{\alpha }^{j+1}=\left[2A_{\delta }\mathcal{Z}+(1-\theta )\left(\mathcal{U}+\mathcal{V}\right)\right]{\alpha }^j+{\mathcal{G}}^{j+1},$$

(34)

where ${\mathcal{G}}^{j+1}=-A_{\delta }{\tilde{\mathcal{Z}}}^{j+1}+\theta \left({\tilde{\mathcal{U}}}^{j+1}+{\tilde{\mathcal{V}}}^{j+1}\right)+2A_{\delta }\widetilde{{\mathcal{Z}}^j}+(1-\theta )(\widetilde{{\mathcal{U}}^j}+\widetilde{{\mathcal{V}}^j})+{\mathbf{B}}^{j+1}-A_{\delta }{w}^{j-1}-A_{\delta }{\sum }_{k=1}^j\left[w^{j-k+1}-2w^{j-k}+w^{j-k-1}\right]B\left(k\right)$.

Now From Equation (34) one can write

$${\alpha }^{j+1}={\mathcal{C}}^{-1}\mathcal{T}{\alpha }^j+{\mathcal{C}}^{-1}{\mathcal{G}}^{j+1},$$

(35)

where $\mathcal{C}=\left[A_{\delta }\mathcal{Z}-\theta \left(\mathcal{U}+\mathcal{V}\right)\right],\mathcal{T}=2A_{\delta }\mathcal{Z}+(1-\theta )\left[\mathcal{U}+\mathcal{V}\right).$ Putting Equation (35) in Equation (33) we get

$$w^{j+1}=\mathcal{Z}{\mathcal{C}}^{-1}\mathcal{T}{\alpha }^j+\mathcal{Z}{\mathcal{C}}^{-1}{\mathcal{G}}^{j+1}+{\tilde{\mathcal{Z}}}^{j+1}.$$

(36)

Using Equation (33) in Equation (36) we have

$$w^{j+1}=\mathcal{Z}{\mathcal{C}}^{-1}\mathcal{T}{\mathcal{Z}}^{-1}{w}^j-\mathcal{Z}{\mathcal{C}}^{-1}\mathcal{T}{\mathcal{Z}}^{-1}{\tilde{\mathcal{Z}}}^j+\mathcal{Z}{\mathcal{C}}^{-1}{\mathcal{G}}^{j+1}+{\tilde{\mathcal{Z}}}^{j+1}.$$

(37)

The above equation shows a recurrence relation of a full discretization scheme which allow us refinement in time. If ${\tilde{w}}^{j+1}$ is numerical solution then

$${\tilde{w}}^{j+1}=\mathcal{Z}{\mathcal{C}}^{-1}\mathcal{T}{\mathcal{Z}}^{-1}\widetilde{w^j}-\mathcal{Z}{\mathcal{C}}^{-1}\mathcal{T}{\mathcal{C}}^{-1}{\tilde{\mathcal{Z}}}^j+\mathcal{Z}{\mathcal{C}}^{-1}{\mathcal{G}}^{j+1}+{\tilde{\mathcal{Z}}}^{j+1}.$$

(38)

Let $e^{j+1}=w^{j+1}-{\tilde{w}}^{j+1}$ be the error at ${(j+1)}^{th}$ time level. Subtracting Equation (37) from Equation (38) then

$$e^{j+1}=\mathsf{\Lambda }{e}^j,$$

where $\mathsf{\Lambda }=\mathcal{Z}{\mathcal{C}}^{-1}\mathcal{T}{\mathcal{Z}}^{-1}$ is the amplification matrix. According to Lax-Richtmyer criterion, the scheme will be stable if $\parallel \mathsf{\Lambda }\parallel \le 1$. It has been verified computationally that $\parallel \mathsf{\Lambda }\parallel \le 1$. For $J=1$ the spectral radius is $0.01025$ which lies in the stability domain.

5. Convergence Analysis

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

Lemma 1

(see [24]). If $w\left(x\right)\in L^2\left(R\right)$ with $\left|w^{\prime }\left(x\right)\right|\le \rho$, for all $x\in (0,1)$, $\rho >0$ and $w\left(x\right)={\sum }_{i=0}^{\infty }{a}_i{\mathcal{H}}_i\left(x\right)$ then $\mid a_i\mid \le \frac{\rho }{2^{j+1}}.$

Lemma 2

(see [25]). If $f(x,y)$ satisfies a Lipschitz condition on $[0,1]\times [0,1]$, that is, there exists a positive L such that for all $(x_1,y),(x_2,y)\in [0,1]\times [0,1]$ we have $\mid f(x_1,y)-f(x_2,y)\mid \le L\mid x_1-x_2\mid$ then

$$a_{i,l}^2\le \frac{L^2}{2^{4j+4}{m}^2}$$

Theorem 1.

If $w\left(x\right)$ and $w_{2M}\left(x\right)$ are the exact and approximate solution of Equation (1), then the error norm $\Vert E_J\Vert$ at $J^{th}$ resolution level is

$$\Vert E_J\Vert \le \frac{4\rho }{3}{\left(\frac{1}{2^{J+1}}\right)}^2.$$

(39)

Proof. 

See [26]. ☐

Theorem 2.

Assume $w(x,y)$ and $w_{2M}(x,y)$ be the exact and approximate solution of Equation (3), then

$$\Vert E_J\Vert \le \frac{L}{4\sqrt{255}}\frac{1}{2^{4J}}.$$

(40)

Proof. 

See [27]. ☐

6. 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_{\infty }$ and $L_2$ error norm are figured out which are defined as follows:

$$L_2=\sqrt{\sum _{i=1}^{2M}{\left(w^{ext}-w^{app}\right)}^2},L_{\infty }=\underset{1\le i\le 2M}{max}|w^{ext}-w^{app}|,$$

(41)

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

$$\begin{array}{cc}\hfill {}^c{D}_t^{\delta }w(x,t)& =-w_t(x,t)+w_{xx}(x,t)+\mathcal{A}(x,t),x\in [0,1],t\in [0,1],1<\delta \le 2,\end{array}$$

(42)

with $\mathcal{A}(x,t)=\frac{2x(1-x)t^{2-\delta }}{\mathsf{\Gamma }(3-\delta )}+2tx(1-x)+2t^2$. 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,\delta =1.1,1.3,1.5,1.7,1.9$. The obtained error norms are shown in

Table 1. Error norms of problem 5.1 for at $J=4$.

$\mathit{\delta }$	$\mathit{t}=0.01,\mathit{\tau }=0.0001$		$\mathit{t}=0.1,\mathit{\tau }=0.001$		$\mathit{t}=1,\mathit{\tau }=0.01$
	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$
1.1	7.0694 × 10${}^{-8}$	2.9496 × 10${}^{-7}$	1.0799 × 10${}^{-5}$	4.5921 × 10${}^{-5}$	2.1556 × 10${}^{-3}$	8.9537 × 10${}^{-3}$
1.3	3.1776 × 10${}^{-8}$	1.3294 × 10${}^{-7}$	7.6592 × 10${}^{-6}$	3.2979 × 10${}^{-5}$	2.1082 × 10${}^{-3}$	8.7615 × 10${}^{-3}$
1.5	1.1646 × 10${}^{-8}$	4.8890 × 10${}^{-8}$	4.8457 × 10${}^{-6}$	2.1318 × 10${}^{-5}$	2.0653 × 10${}^{-3}$	8.5871 × 10${}^{-3}$
1.7	5.2296 × 10${}^{-9}$	2.1899 × 10${}^{-8}$	3.3635 × 10${}^{-6}$	1.4944 × 10${}^{-5}$	2.1431 × 10${}^{-3}$	8.8989 × 10${}^{-3}$
1.9	2.2087 × 10${}^{-9}$	9.2078 × 10${}^{-9}$	2.1386 × 10${}^{-6}$	9.3949 × 10${}^{-5}$	2.4382 × 10${}^{-3}$	1.0094 × 10${}^{-2}$

. From table it is obvious that results of the present scheme match well with exact solution. Also in

Table 2. Comparison of maximum error of problem 5.1 with previous work at $t=1$ and $\delta =1.7$.

[28]				Present Method
$\mathit{h}$	$\mathit{\tau }$	$\mathit{Error}$		$\mathit{dx}$	$\mathit{\tau }$	$\mathit{Error}$
0.05	0.05	4.4333 × 10${}^{-3}$		1/4	0.05	8.2306 × 10${}^{-4}$
0.025	0.0125	7.7368 × 10${}^{-4}$		1/8	…	5.4184 × 10${}^{-4}$
0.0125	0.00625	3.1827 × 10${}^{-4}$		1/16	…	4.9195 × 10${}^{-4}$

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. Error norms of problem 5.1 for different values of $\tau$ and $\delta$.

$\mathit{\delta }=1.5$				$\mathit{\delta }=1.7$
$\mathit{\tau }$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$		${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$
1/4	5.4216 × 10${}^{-2}$	2.2515 × 10${}^{-1}$		5.4689 × 10${}^{-2}$	2.2708 × 10${}^{-1}$
1/8	2.7891 × 10${}^{-2}$	1.1571 × 10${}^{-1}$		2.8861 × 10${}^{-2}$	1.1961 × 10${}^{-1}$
1/16	1.3645 × 10${}^{-2}$	5.6647 × 10${}^{-2}$		1.4343 × 10${}^{-2}$	5.9443 × 10${}^{-2}$
1/32	6.6683 × 10${}^{-3}$	2.7699 × 10${}^{-2}$		7.0022 × 10${}^{-3}$	2.9034 × 10${}^{-2}$
1/64	3.2674 × 10${}^{-3}$	1.3580 × 10${}^{-2}$		3.4061 × 10${}^{-3}$	1.4135 × 10${}^{-2}$

shows convergence in time for fixed $dx=1/32$. The convergence rate of the proposed scheme has been addressed in

Table 4. Convergence rate of maximum error of problem 5.1 at $t=1$ and $\delta =1.7$.

J	$\mathit{\tau }$	$\mathit{Error}$	$\mathit{Rate}$
1	1/10	2.6194 × 10${}^{-3}$	-
2	1/20	1.3086 × 10${}^{-3}$	1.0012
3	1/40	5.5920 × 10${}^{-4}$	1.2265
4	1/80	2.2921 × 10${}^{-4}$	1.2866
5	1/160	9.3261 × 10${}^{-5}$	1.2973
6	1/320	3.7377 × 10${}^{-5}$	1.3191

. 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

$$\begin{array}{cc}\hfill {}^c{D}_t^{\delta }w(x,t)& =-w_t(x,t)+w_{xx}(x,t)+\mathcal{A}(x,t),x\in [0,1],t\in [0,1],1<\delta \le 2,\end{array}$$

(43)

coupled with initial and boundary conditions

$$\begin{array}{c}\left\{\begin{array}{c}w(x,0)=0,w_t(x,0)=0x\in (0,1)\hfill \\ w(0,t)=t^3,w(1,t)=e{t}^3,t\in [0,1].\hfill \end{array}\right.\hfill \end{array}$$

(44)

The exact solution and source term are given by $w(x,t)=e^x{t}^3$ and $\mathcal{A}(x,t)=\frac{6t^{3-\delta }{e}^x}{\mathsf{\Gamma }(4-\delta )}+3t^2{e}^x-t^3{e}^x$. In

Table 5. Error norms of problem 5.2 at $J=4$.

$\mathit{\delta }$	$\mathit{t}=0.01,\mathit{\tau }=0.0001$		$\mathit{t}=0.1,\mathit{\tau }=0.001$
	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$
1.1	1.7079 × 10${}^{-7}$	6.8446 × 10${}^{-7}$	1.2504 × 10${}^{-4}$	5.3397 × 10${}^{-4}$
1.3	6.5331 × 10${}^{-7}$	2.5683 × 10${}^{-6}$	4.4278 × 10${}^{-4}$	1.8777 × 10${}^{-3}$
1.5	1.2494 × 10${}^{-6}$	4.7989 × 10${}^{-6}$	8.7071 × 10${}^{-4}$	3.6354 × 10${}^{-3}$
1.7	1.3386 × 10${}^{-6}$	5.0827 × 10${}^{-6}$	1.0489 × 10${}^{-3}$	4.2541 × 10${}^{-3}$
1.9	5.6739 × 10${}^{-7}$	2.1936 × 10${}^{-6}$	5.1085 × 10${}^{-4}$	2.0046 × 10${}^{-3}$

the obtained error norms are shown for parameters $t=0.01,0.1,$ $\delta =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]

$$\begin{array}{cc}\hfill {}^c{D}_t^{\delta }w(x,y,t)& =\Delta w(x,y,t)+\mathcal{B}(x,y,t),(x,y)\in [0,1]\times [0,1],t\in [0,1],1<\delta \le 2,\hfill \end{array}$$

(45)

with exact solution $w(x,y,t)=sin\left(\pi x\right)sin\left(\pi y\right)t^{\delta +3}$, and source term

$$\mathcal{B}(x,y,t)=sin\left(\pi x\right)sin\left(\pi y\right)\left[\frac{\mathsf{\Gamma }(\delta +3)t^2}{2}-2t^{\delta +2}\right].$$

We solved this problem for resolution level $J=4$ and the obtained results are recorded in

Table 6. Comparison of problem 5.4 at $t=1$ and $\delta$ with previous results.

$\mathit{\delta }$	$\mathit{t}=0.1,\mathit{\tau }=0.001$		$\mathit{t}=0.2,\mathit{\tau }=0.01$		$\mathit{t}=0.5,\mathit{\tau }=0.05$
	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$
1.5	1.6049 × 10${}^{-4}$	8.0439 × 10${}^{-5}$	4.3534 × 10${}^{-4}$	2.1819 × 10${}^{-4}$	2.5390 × 10${}^{-2}$	1.2725 × 10${}^{-2}$
1.7	1.1635 × 10${}^{-4}$	5.8320 × 10${}^{-5}$	5.9673 × 10${}^{-4}$	2.9908 × 10${}^{-4}$	7.5824 × 10${}^{-3}$	3.8003 × 10${}^{-3}$
1.9	3.0965 × 10${}^{-5}$	1.5519 × 10${}^{-5}$	3.9390 × 10${}^{-4}$	1.9742 × 10${}^{-4}$	4.8842 × 10${}^{-3}$	2.4480 × 10${}^{-3}$

for different values of time and $\tau$. From Table 6 it is clear that the proposed scheme works well for the solution of two dimensional problems.

Table 7. Error norms of problem 5.3 for different values of $\tau$ and $\delta$.

		${\mathit{L}}_{\infty }$
$\mathit{\delta }$	$\mathit{\tau }$	$\mathit{Present}$	[29]
1.25	1/10	8.1748 × 10${}^{-3}$	8.1577 × 10${}^{-2}$
	1/20	6.5092 × 10${}^{-3}$	3.4379 × 10${}^{-2}$
	1/40	5.7150 × 10${}^{-3}$	1.4484 × 10${}^{-2}$
1.5	1/10	6.7087 × 10${}^{-3}$	2.9942 × 10${}^{-2}$
	1/20	4.8922 × 10${}^{-3}$	1.0749 × 10${}^{-2}$
	1/40	4.1390 × 10${}^{-3}$	3.8291 × 10${}^{-3}$
1.75	1/10	6.7087 × 10${}^{-3}$	8.4482 × 10${}^{-3}$
	1/20	4.8922 × 10${}^{-3}$	2.5877 × 10${}^{-3}$
	1/40	4.1390 × 10${}^{-3}$	7.8500 × 10${}^{-4}$

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]

$$\begin{array}{cc}\hfill {}^c{D}_t^{\delta }w(x,t)+w(x,t)& =w_{xx}(x,t)+\mathcal{A}(x,t),x\in [0,1],t\in [0,1],1<\delta \le 2,\end{array}$$

(46)

coupled with initial and boundary conditions

$$\begin{array}{c}\left\{\begin{array}{c}w(x,0)=0,w_t(x,0)=0x\in (0,1)\hfill \\ w(0,t)=o,w(1,t)=0,t\in [0,1],\hfill \end{array}\right.\hfill \end{array}$$

(47)

where the forcing terms are $\mathcal{A}(x,t)=\frac{2t^{2-\delta }x(1-x)}{\mathsf{\Gamma }(3-\delta )}+t^{2}x(1-x)-2t^2.$ This problem has been solved with the help of the proposed scheme. In

Table 8. Absolute error at different points of example 5.4 at $\tau =0.001$.

$(\mathit{x},\mathit{t})$	$\mathit{\delta }=1.1$		$\mathit{\delta }=1.3$		$\mathit{\delta }=1.5$		$\mathit{\delta }=1.9$
	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\infty }$ [19]	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\infty }$ [19]	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\infty }$ [19]	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\infty }$ [19]
$(0.1,0.1)$	2.1684 × 10${}^{-19}$	9.5133 × 10${}^{-9}$	1.0842 × 10${}^{-19}$	6.6004 × 10${}^{-9}$	3.2526 × 10${}^{-19}$	4.4920 × 10${}^{-9}$	5.8546 × 10${}^{-18}$	1.9326 × 10${}^{-9}$
$(0.2,0.2)$	9.5409 × 10${}^{-18}$	1.0530 × 10${}^{-7}$	2.1684 × 10${}^{-17}$	7.9127 × 10${}^{-8}$	3.4694 × 10${}^{-18}$	5.7844 × 10${}^{-8}$	1.8735 × 10${}^{-16}$	2.8903 × 10${}^{-8}$
$(0.3,0.3)$	4.1633 × 10${}^{-17}$	9.6665 × 10${}^{-7}$	4.1633 × 10${}^{-17}$	3.3461 × 10${}^{-7}$	1.0755 × 10${}^{-16}$	2.5678 × 10${}^{-7}$	1.3634 × 10${}^{-15}$	1.4105 × 10${}^{-7}$
$(0.4,0.4)$	1.3877 × 10${}^{-17}$	1.0813 × 10${}^{-6}$	1.8735 × 10${}^{-16}$	9.1574 × 10${}^{-7}$	7.0776 × 10${}^{-16}$	7.3594 × 10${}^{-7}$	4.5033 × 10${}^{-15}$	4.3402 × 10${}^{-7}$
$(0.5,0.5)$	1.3877 × 10${}^{-16}$	2.2190 × 10${}^{-6}$	3.9551 × 10${}^{-16}$	1.6516 × 10${}^{-6}$	2.0261 × 10${}^{-15}$	1.6516 × 10${}^{-6}$	1.4231 × 10${}^{-14}$	1.0367 × 10${}^{-6}$

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.

• Problem 5.5:

Now we consider the following equation

$$\begin{array}{c}{}^c{D}_t^{\delta }w(x,y,t)=a_1\Delta w(x,y,t)-b_{1}sin\left(w(x,y,t)\right),(x,y)\in \mathsf{\Phi },t\in [0,T],1<\delta \le 2,\end{array}$$

(48)

where $a_1$ and $b_1$ are constants and the initial and boundary conditions are

$$\begin{array}{c}\left\{\begin{array}{c}w(x,y,0)=arctan\left(exp(\frac{1}{2}-\sqrt{15x^2+15y^2)}\right),w_t(x,y,0)=0,(x,y)\in \tilde{\mathsf{\Phi }}=\mathsf{\Phi }\cup \partial \mathsf{\Phi },\hfill \\ w(x,y,t)=0,(x,y)\in \partial \mathsf{\Phi },t\in [0,T].\hfill \end{array}\right.\end{array}$$

(49)

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.

7. 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.