Meshless method based on RBFs for solving three-dimensional multi-term time fractional PDEs arising in engineering phenomenons

1 Introduction

Fractional PDE models are widely employed in many disciplines of physics, and variety of methods are used in the literature for their simulations (see Ain and He, 2019; Nadeem and Li, 2019; He and El-Dib, 2020a). Nowadays, it is in the limelight of active researchers with extensive research to advance numerical and analytical solutions for linear and nonlinear fractional PDEs (He and El-Dib, 2020b; He and Ain, 2020; He et al., 2012; Ahmad et al., 2018, 2020a,b,c; Srivastava et al., 2020; He, 2018; Inc et al., 2020; He, 2014). In the current work, we have considered the following two- and three-term time fractional Sobolev equation

Recent literature uses several meshless techniques to solve numerically different PDE models in almost all disciplines of physics and mathematics. In particular, the meshless techniques using radial basis function (RBF) are the main one of these processes. The meshless property makes it popular among researchers. These methods overcome the challenges of dimensionality using conventional methods. In contrast to mesh-based methods, meshless approaches do not need mesh in the domain. These techniques can compute the solution utilizing uniform or non-uniform nodes in both regular and irregular computational domains, which increases the applicability and usefulness of the methods. Meshless methods are in fact actionable and valuable that can be applied to solve physical problems (Wang et al., 2021; Hussain et al., 2020; Khan et al., 2020; Ahmad et al., 2019a, 2020d,e,f,g; Wang et al., 2019, 2020; Nawaz et al., 2021; Wang and Zheng, 2016).

Like other numerical methods, the meshless methods, have some drawbacks, which may be the most important one to choose the optimal value of the shape-parameter and dense ill-conditioned matrices. To circumvent these shortcomings, the local meshless technique is the best alternative proposed by researchers which is precise and stable in finding solutions for various integer and fractional PDE models (Ahmad et al., 2019b,c). These methods produce well conditioned sparse matrices and are less sensitive to the selection of shape-parameters than the global version. In addition, the local meshless method is extra valuable and effective than its global counterpart. Recently, these methods are explored in various form in different applications (Ahmad et al., 2017, 2020h; Shu, 2000).

In this paper, the local meshless technique is utilized for the solution of model Eqs. (1), (2). Inverse multiquadric (IMQ) radial basis functions (RBFs) are taken into account. Moreover, Two kinds of irregular domains are analyzed in numerical examinations.

2 Implementation of local meshless technique

According to the suggested technique, to approximated the derivatives of $\mathcal{U}(\bar{\mathbf{x}},t)$ at the centers ${\bar{\mathbf{x}}}_h$ by the neighborhood of ${\bar{\mathbf{x}}}_h,\{{\bar{\mathbf{x}}}_{h1},{\bar{\mathbf{x}}}_{h2},{\bar{\mathbf{x}}}_{h3},\dots ,{\bar{\mathbf{x}}}_{{\mathit{hn}}_h}\}\subset \{{\bar{\mathbf{x}}}_1,{\bar{\mathbf{x}}}_2,\dots ,{\bar{\mathbf{x}}}_{N^n}\},n_h\ll N^n$ , where $h=1,2,\dots ,N^n$ . We have used $\bar{\mathbf{x}}=x,\bar{\mathbf{x}}=(x,y)$ and $\bar{\mathbf{x}}=(x,y,z)$ for one-, two- and three-dimensional cases respectively.

Procedure for one-dimensional case is

Substituting the IMQ-RBF $\psi (\Vert y-y_p\Vert )=\frac{1}{\sqrt{1+{(c\Vert y_{\mathit{hk}}-y_p\Vert )}^2}}$ in (5)

Matrix form of (6)

From (8), we obtain

(5) and (9) implies $${\mathcal{U}}^{\left(m\right)}\left(x_h\right)={\left({\mathit{\lambda }}_{n_h}^{\left(m\right)}\right)}^T{\mathbf{V}}_{n_h},$$ where $${\mathbf{U}}_{n_h}={\left[\mathcal{U}\left(x_{h1}\right),\mathcal{U}\left(x_{h2}\right),\dots ,\mathcal{U}\left(x_{{\mathit{hn}}_h}\right)\right]}^T\text{.}$$

The derivatives of $\mathcal{U}(x,y,t)$ w.r.t. x and y can be found as $${\mathcal{U}}_x^{\left(m\right)}(x_h,y_h)\approx {\displaystyle \sum _{k=1}^{n_h}}{\gamma }_k^{\left(m\right)}\mathcal{U}(x_{\mathit{hk}},y_{\mathit{hk}}),h=1,2,\dots ,N^2,$$ $${\mathcal{U}}_y^{\left(m\right)}(x_h,y_h)\approx {\displaystyle \sum _{k=1}^{n_h}}{\eta }_k^{\left(m\right)}\mathcal{U}(x_{\mathit{hk}},y_{\mathit{hk}}),h=1,2,\dots ,N^2\text{.}$$

For ${\gamma }_k^{\left(m\right)}$ and ${\eta }_k^{\left(m\right)}$ ( $k=1,2,\dots ,n_h$ ), we continue as $${\mathit{\gamma }}_{n_h}^{\left(m\right)}={\mathbf{A}}_{n_h}^{-1}{\mathbf{\Phi }}_{n_h}^{\left(m\right)},$$ $${\mathit{\eta }}_{n_h}^{\left(m\right)}={\mathbf{A}}_{n_h}^{-1}{\mathbf{\Phi }}_{n_h}^{\left(m\right)}\text{.}$$

Caputo derivative (Caputo, 1967) for ${\gamma }_1\in (0,1)$ is utilized for time derivative $\frac{{\partial }^{{\gamma }_1}\mathcal{U}(\bar{\mathbf{x}},t)}{\partial t^{{\gamma }_1}}$ as follows $$\frac{{\partial }^{{\gamma }_1}\mathcal{U}(\bar{\mathbf{x}},t)}{\partial t^{{\gamma }_1}}=\left\{\begin{array}{cc}\frac{1}{\mathrm{\Gamma }\left(1-{\gamma }_1\right)}{\int }_0^t\frac{\partial \mathcal{U}(\bar{\mathbf{x}},\vartheta )}{\partial \vartheta }{\left(t-\vartheta \right)}^{-{\gamma }_1}d\vartheta ,\hfill & 0<{\gamma }_1<1\hfill \\ \hfill & \hfill \\ \frac{\partial \mathcal{U}(\bar{\mathbf{x}},t)}{\partial t},\hfill & {\gamma }_1=1\text{.}\hfill \end{array}\right.$$

To compute the derivative term we can proceed as follows, where $t_q=q\tau ,q=0,1,2,\dots ,Q$ and time step size $\mathrm{\Delta }\tau$ in $[0,t]$ . $$\begin{array}{ccc}\hfill & \frac{{\partial }^{{\gamma }_1}\mathcal{U}(\bar{\mathbf{x}},t_{q+1})}{\partial t^{{\gamma }_1}}\hfill & =\frac{1}{\mathrm{\Gamma }(1-{\gamma }_1)}{\int }_0^{t_{q+1}}\frac{\partial \mathcal{U}(\bar{\mathbf{x}},\vartheta )}{\partial \vartheta }{\left(t_{q+1}-\vartheta \right)}^{-{\gamma }_1}d\vartheta ,\hfill \\ \hfill & =\frac{1}{\mathrm{\Gamma }(1-{\gamma }_1)}{\displaystyle \sum _{r=0}^q}{\int }_{r\mathrm{\Delta }\tau }^{(r+1)\mathrm{\Delta }\tau }\frac{\partial \mathcal{U}(\bar{\mathbf{x}},\vartheta )}{\partial \vartheta }{\left(t_{r+1}-\vartheta \right)}^{-{\gamma }_1}d\vartheta ,\hfill \\ \hfill & \approx \frac{1}{\mathrm{\Gamma }(1-{\gamma }_1)}{\displaystyle \sum _{r=0}^q}{\int }_{r\mathrm{\Delta }\tau }^{(r+1)\mathrm{\Delta }\tau }\frac{\partial \mathcal{U}(\bar{\mathbf{x}},{\vartheta }_r)}{\partial \vartheta }{\left(t_{r+1}-\vartheta \right)}^{-{\gamma }_1}d\vartheta \text{.}\hfill \end{array}$$

The term $\frac{\partial \mathcal{U}(\bar{\mathbf{x}},{\vartheta }_r)}{\partial \vartheta }$ can be approximated as $$\frac{\partial \mathcal{U}(\bar{\mathbf{x}},{\vartheta }_r)}{\partial \vartheta }=\frac{\mathcal{U}(\bar{\mathbf{x}},{\vartheta }_{r+1})-\mathcal{U}(\bar{\mathbf{x}},{\vartheta }_r)}{\vartheta }+\mathcal{O}\left(\mathrm{\Delta }\tau \right)\text{.}$$

Then $$\begin{array}{cc}\hfill \frac{{\partial }^{{\gamma }_1}\mathcal{U}(\bar{\mathbf{x}},t_{q+1})}{\partial t^{{\gamma }_1}}& \approx \frac{1}{\mathrm{\Gamma }(1-{\gamma }_1)}{\displaystyle \sum _{r=0}^q}\frac{\mathcal{U}(\bar{\mathbf{x}},t_{r+1})-\mathcal{U}(\bar{\mathbf{x}},t_r)}{\mathrm{\Delta }\tau }{\int }_{r\mathrm{\Delta }\tau }^{(r+1)\mathrm{\Delta }\tau }{\left(t_{r+1}-\vartheta \right)}^{-{\gamma }_1}d\vartheta ,\hfill \\ \hfill & =\frac{1}{\mathrm{\Gamma }(1-{\gamma }_1)}{\displaystyle \sum _{r=0}^q}\frac{\mathcal{U}(\bar{\mathbf{x}},t_{q+1-r})-\mathcal{U}(\bar{\mathbf{x}},t_{q-r})}{\mathrm{\Delta }\tau }{\int }_{r\mathrm{\Delta }\tau }^{(r+1)\mathrm{\Delta }\tau }{\left(t_{r+1}-\vartheta \right)}^{-{\gamma }_1}d\vartheta ,\hfill \\ \hfill & =\left\{\begin{array}{cc}\frac{\mathrm{\Delta }{\tau }^{-{\gamma }_1}}{\mathrm{\Gamma }(2-{\gamma }_1)}({\mathcal{U}}^{q+1}-{\mathcal{U}}^q)+\frac{\mathrm{\Delta }{\tau }^{-{\gamma }_1}}{\mathrm{\Gamma }(2-{\gamma }_1)}{\displaystyle \sum _{r=1}^q}({\mathcal{U}}^{q+1-r}-{\mathcal{U}}^{q-r})[{(r+1)}^{1-{\gamma }_1}-r^{1-{\gamma }_1}],\hfill & q⩾1\hfill \\ \frac{\mathrm{\Delta }{\tau }^{-{\gamma }_1}}{\mathrm{\Gamma }(2-{\gamma }_1)}({\mathcal{U}}^1-{\mathcal{U}}^0).\hfill & q=0,\hfill \end{array}\right.\hfill \end{array}$$

Letting $a_0=\frac{\mathrm{\Delta }{\tau }^{-{\gamma }_1}}{\mathrm{\Gamma }(2-{\gamma }_1)}$ and $b_r={(r+1)}^{1-{\gamma }_1}-r^{1-{\gamma }_1}$ ,   $r=0,1,\dots ,q$ , we have

The fractional derivative of order ${\gamma }_2$ and ${\gamma }_3$ can be found as above.

3 Results

The local meshless technique is tested for applicability, accuracy and efficiency to approximate the solution of model Eqs. (1), (2). In this paper, different computational domains are utilized with uniform and scatted nodes. Throughout the paper, we have employed the Crank-Nicolson scheme, IMQ RBF with $c=100,\mathrm{\Delta }\tau =0.005$ , and spatial domain $[0,1]$ unless specifically stated. The accuracy is measure as follows

Close

Fig. 1

Problem 1, c versus error norms (left) two-term model equation, (right) three-term model equation.

Export to PPT

Close

Fig. 2

Domain-1 (left) and Domain-2 (right).

Export to PPT

Close

Fig. 3

Problem 1, error profile of two-term model equation (left) Domain-1 (right) Domain-2.

Export to PPT

Close

Fig. 4

Problem 1, error profile of three-term model equation (left) Domain-1 (right) Domain-2.

Export to PPT

Close

Fig. 5

Problem 3, two-term model equation at $z=0.5$ (left) numerical solution (right) exact solution.

Export to PPT

Close

Fig. 6

Problem 3, three-term model equation at $z=1$ (left) numerical solution (right) exact solution.

Export to PPT

4 Conclusion

In this paper, we have examined the effectiveness and applicability of the suggested local meshless technique to compute the numerical solution of time fractional Sobolev model equations. The obtained results prove that the suggested technique works amazingly in fractional PDE models. Numerical experiments show that the algorithm gives good accuracy and in light of these analyses, we suggest that the local meshless technique can be implemented to such types of fractional PDE models which appear in physical problems.