Meshfree numerical approach for some time-space dependent order partial differential equations in porous media

1.   Introduction

The advection-dispersion equation has arisen in many chemical engineering and physical models ^[1]. In earth sciences, it is used as a solute transport water flow in surface and subsurface, stream, groundwater, ocean currents and deep river flows. Solute transport in streams, groundwater and rivers is controlled by heterogeneity or physical features in different levels. In the past, the advection-dispersion equation has been successfully used for mobile-immobile or transient storage models. In recent times, researchers highlighted the necessity of transport models from which heterogeneity and connectivity of spatial features inside solute transport can be described better. In porous media, the transport flow is controlled by advection and dispersion equations which will predict a breakthrough curve ^[2].

The study of numerical models of solute transport is a basic aspect in parameter reconnaissance at both field and micro scales. The characterization of motion of solute transport is measured by the application of the advection-dispersion equation in porous media. However, there is some evidence that the advection-dispersion model failed to explain transport in heterogeneous, fractured and even in homogeneous media. The advection-dispersion and mobile immobile have fitted in breakthrough curves ^[3]. In both fractured and porous media, the results show that the mobile-immobile model is better than the advection-dispersion model, especially in explaining long tails and peaks of breakthrough curves. Because of such reasons, in tough-walled fracture, the mobile-immobile model is more effective than in smooth-walled fracture. The single-rate mobile-immobile model appears as the advection-dispersion model transport in the mobile region, whereas the transport in the immobile region by diffusion, making physical non-equilibrium ^[3,4,5,6,7].

Most scientific phenomena are formed through the fractional PDEs, which are either linear or nonlinear. PDEs with fractional orders are now becoming a useful study in the fields of physics, acoustics, engineering, chemistry, viscoelasticity, and biology ^{[8,9,10,11,12]}.

For solute transport, the fractional mobile immobile advection-dispersion equation assumes power law waiting in the immobile region, yielding a fractional order in terms of time derivative. The fractional-order models are equipollent to non-fractional mobile-immobile transport with power law memory functions. There are some numerical methods used to study the mobile-immobile advection-dispersion model (see ^{[13,14,15,16,17]}). The advection-dispersion mobile-immobile model ^[13] is given by

with the following initial and boundary conditions:

where $r_1$, $r_2$, $r_3$ and $r_4$ are constants, and $R(\nu, \tau)$ is known function. Also, $0 < \zeta(\nu, \tau)\leq 1$, $1 < \xi_1(\nu, \tau)\leq 2$ and $0 < \xi_2(\nu, \tau)\leq 1$ are the time and space-dependent derivative orders of corresponding terms of Eq (1.1), $\mathbb{D}_\nu^{\xi(\nu, \tau)}$ w.r.t. $\nu$ in domain $D$.

There are two types of numerical methods to deal with the PDE models. One is meshgrid method and the other one is meshfree method (see ^{[18,19,20,21,22,23]}). Currently, the meshfree methods are getting more attention of researchers to find numerical solutions to the PDEs, especially in fractional orders. Because meshfree methods do not require meshing, these methods only need uniform distributed nodes in the domain. The accuracy of meshfree methods is higher and less time-consuming than meshgrid methods. The homotopy analysis ^[24,25], Variational iteration method ^[26], Spectral method ^[27], Adomian decomposition method ^[28,29] and Radial basis function method ^[30,31] are some well known meshfree methods.

In this research article, our focus is on the implementation of the radial basis function (RBF) method. Its implementation is easy in initial and boundary value problems and the choice of generating functions is also open. The meshfree radial basis function (RBF) method is used in both local and global forms. The RBF method is implemented through some generating functions such as multiquadric, inverse multiquadric, inverse quadric and Gaussian. The main issue of these functions is the choice of shape parameter, which is an open choice and necessary for the best numerical results. This issue sometimes creates problems to find the numerical results, especially in higher dimensions. However, in Gaussian, the choice of shape parameter is not very challenging. The shape parameter $c = \frac{1}{\sqrt{2}\sigma^2}$, where $\sigma$ is the variance between scattered data $\{\nu_i, i = 1, 2, \cdots, M\}$ and mean or center $\stackrel{\text{---}}{\nu}$ of the data (see ^[32]). Many researchers implemented RBF method on various PDE models. For example, Hosseini et al. ^[33] used RBF method to find the approximate solution of the fractional telegraph equation. Samad et al. ^[34] used RBF method to find the numerical solution of multi-term fractional-order PDEs. Similarly, Wang et al. ^[35,36], Khan et al. ^[37], Ali et al. ^[38] and Nikan et al. ^[19,20] have used meshfree RBF method to test the numerical results of various PDE models.

The current article covers the study of numerical solution of some time-space-dependent order PDE models using the Gaussian radial basis function (GA-RBF). The space derivatives are dealt with by Riesz fractional derivative and Grünwald-Letnikov derivative definitions. The time derivatives are iterated by the Caputo fractional derivative. The stability of the scheme is discussed in Section 3, whereas the numerical results are given in Section 4.

2.   Materials and methods

2.1.   Time derivatives

Let $\delta\tau$ be the time interval for $\tau\in[0, T]$, then $\tau_n = n\delta\tau$ for $n = 0, 1, 2, \cdots, N$.

Definition 2.1. For any variable order $\zeta(\nu, \tau)$, the Coimbra operator $\mathbb{D}$ is defined ^[39] by

where $0 < \zeta(\nu, \tau) < 1$. The Caputo fractional derivative can be expressed from Coimbra fractional derivative operator if $w(\nu, \tau_{0_{+}}) = w(\nu, \tau_{0_{-}})$, i.e.,

Using finite difference method to discretize the time fractional term of Eq (1.1),

Thus, the time discretization for Eq (1.1) is

where $\delta_k = (k+1)^{1-\zeta(\nu, \tau)}-k^{1-\zeta(\nu, \tau)}$ and the truncation error $\rho_{\delta\tau}^{n+1}$ is bounded (see ^[34]).

Lemma 2.1. $\delta_k = (k+1)^{1-\zeta(\nu, \tau)} - k^{1-\zeta(\nu, \tau)}$ has the following properties for $k = 0, 1, 2, \cdots$.

(i) $\delta_0 = 1$, $\delta_k > 0$,

(ii) $\delta_k \geq \delta_{k+1}$.

Let $R_1 = \frac{r_1}{\delta\tau}$ and $R_2 = \frac{r_2\delta\tau^{-\zeta(\nu, \tau)}}{\Gamma(2-\zeta(\nu, \tau))}$, then Eq (2.1) is written as

2.2.   Space discretization

2.2.1.   RBF method

Let the solution of Eq (1.1) at $n$-time level be

where $r_m = \|\nu-\nu_m\|$, $\nu_m$ are the centers of $\nu$, and $\|\cdot\|$ is the Euclidean norm. Then, for the collocation points $\{\nu_i: 1\leq i\leq M\}$, Eq (2.3) can be written as

where $r_{im} = \|\nu_i-\nu_m\|$. The Gaussian function ^[32] for the proposed method is given by

where $c$ is the shape parameter. The value of the shape parameter can be adjusted ^[32] by the equation $c = \frac{1}{\sqrt{2}\sigma^2}$, where $\sigma^2$ is the variance between $\nu_i$ and the centers $\nu_m$. Thus, the matrix form of Eq (2.4) is obtained by

Let $A = \{\exp(-cr^2_{im}), 1\leq i, m\leq M\}$, $w^n = \{w(\nu_1, \tau_n), w(\nu_2, \tau_n), \cdots, w(\nu_M, \tau_n)\}^T$ and $\mu^n = \{\mu_1^n, \mu_2^n, \cdots, \mu_M^n\}^T$. Then, equation (2.5) is obtained by

Since $A$ is invertible (see ^{[33,34,35,36]}),

Also, $A$ can be split into $A_b$ and $A_d$ such that $A_b = \{\varphi(r_{im}), i = 1, M, 1\leq m\leq M, 0$ elsewhere} and $A_d = \{\varphi(r_{im}), 2\leq i\leq M-1, 1\leq m\leq M, 0$ elsewhere}.

2.2.2.   Space derivatives

Applying $\theta$-weighted method to Eq (1.1) at time levels $n$ and $n+1$,

where $0 \leq \theta \leq 1$.

Lemma 2.2. If $w(\nu, \tau)\in L_1(\Omega)$ and $\mathbb{D}_\nu^{\xi(\nu, \tau)}w(\nu, \tau)\in C(\Omega),$ where $0\leq p-1 < \xi(\nu, \tau) < p$, then the normalized Grünwald-Letnikov approximation is obtained by

where $\psi_j(\xi(\nu_i, \tau_n)) = \frac{\Gamma(j-\xi(\nu_i, \tau_n))}{\Gamma(j+1)\Gamma(-\xi(\nu_i, \tau_n))}$, $i = 1, 2, \cdots, M$ are the number of collocation nodes, $j = 0, 1, 2, \cdots$.

Lemma 2.3. For $1 \leq i \leq M$, $0 \leq n \leq N$, the coefficients $\psi_j(\xi_i^n)$ satisfy the following properties for $1 < \xi(\nu, \tau) \leq 2$.

(i) $\psi_0(\xi_i^n) = 1$, $\psi_1(\xi_i^n) < 0$,

(ii) $\psi_j(\xi_i^n) > 0$ for $j \geq 2$,

(iii) ${\sum\limits_{j = 0}^{\infty}\psi_j(\xi_i^n) = 0}$, ${\sum\limits_{j = 0}^{\lambda}\psi_j(\xi_i^n) < 0}$ for $\lambda \geq 1$.

Proof. Since for $j \geq 0$, $\psi_j(\xi_i^n) = \frac{\Gamma(j-\xi_i^n)}{\Gamma(j+1)\Gamma(-\xi_i^n)}$, it is easy to verify (i) by putting $j = 0, 1$. Also,

for $1 < \xi_i^n \leq 2$ it is clear from above that $\psi_j(\xi_i^n) > 0$ for $j \geq 0$. Hence, (ii) is satisfied. Also, consider ${(1-y)^{\xi_i^n} = \sum\limits_{j = 0}^{\infty}\psi_j(\xi_i^n)y^j}$ and let $y = 1$, we obtain ${\sum\limits_{j = 0}^{\infty}\psi_j(\xi_i^n) = 0}$. Further, ${\sum\limits_{j = 0}^{\lambda}\psi_j(\xi_i^n) < 0}$ for $\lambda \geq 1$.

Lemma 2.4. For $0 < \xi(\nu, \tau) \leq 1$, the coefficients $\{\psi_j(\xi_i^n), \ \ 1 \leq i \leq M, \ \ 0 \leq n \leq N\}$ satisfy the following properties:

(i) $\psi_0(\xi_i^n) = 1$, $\psi_j(\xi_i^n) < 0$ for $j\geq 1$,

(ii) ${\sum\limits_{j = 0}^{\infty}\psi_j(\xi_i^n) = 0}$, ${\sum\limits_{j = 0}^{\lambda}\psi_j(\xi_i^n) > 0}$ for $\lambda\geq1$.

Proof. Similar to Lemma 2.3.

Now, from Eqs (2.9) and (2.10), and according to the Riesz fractional derivative operator ^[40], we have

where $C_{\xi(\nu_i, \tau_n)} = \frac{1}{-2\cos\frac{\pi\xi(\nu\;_i\;,\; \tau\;_n\;)}{2}}$ . For $i = 1, 2, \cdots, M$ and $m = j+1$, using Eq (2.9), we have

and from Eq (2.10),

where $'\ast'$ shows the product of the $i$-th component of $h_i$ with corresponding row of either $G(i, m)$ or $G_T(i, m)$ at $n$ time level. $h_i^n = C_{\xi(\nu_i, \tau_n)}\ast(\delta\nu)^{-\xi(\nu_i, \tau_n)}$.

and $G_T$ is the transpose of order $M$ matrix $G$. Thus,

Now, using Eq (2.13), we obtain from Eq (2.8) as

Let $W_1^n = G_{\xi_1^n}+{G_{\xi_1^n}}_T$ and $W_2^n = G_{\xi_2^n}+{G_{\xi_2^n}}_T$, Eq (2.14) becomes

2.3.   Numerical scheme

Comparing Eq (2.2) with Eq (2.15), we obtain

After equation of the terms at $n$ and $n+1$ time levels and using Eq (2.6), we obtain

Thus,

where

Let

then

Again using Eqs (2.6) and (2.7), Eq (2.19) can be written as

Thus, Eq (2.20) is the required numerical scheme.

3.   Stability and convergence

The developed scheme for Eq (1.1) is the recurrence relation $w(\tau_n)$ and $w(\tau_{n+1})$ for $n = 0, 1, 2, ....., N$. The elements of the amplification matrix $\mathbb{F} = AP^{-1}QA^{-1}$ are dependent on the constant number $c_i = \frac{\delta \tau^{\zeta}}{(\delta \nu)^{(\xi_1, \xi_2)}}$, where $\zeta$ is the time order and $\xi_1, \xi_2$ are space order differential operators respectively. $\delta \tau$ is the time step size and $\delta \nu$ is spatial step size at successive nodes. Now let us consider the exact solution of Eq (1.1) be $\mathscr{W}^n$ at $\tau^n$.

There are some important theorems from G. E. Fasshauer ^[32], whose statements are given below.

Theorem 3.1.

provided that ${\delta \nu}_{\chi, \Omega} \leq {\delta \nu}_0$, where

Theorem 3.2. Let $\Gamma\subseteq\mathscr{R}^s$ be open and bounded which satisfies the interior cone condition, and let $\Theta\in\mathbb{C}^{2k}(\Gamma\times\Gamma)$ be symmetric and strictly conditionally positive definite of order $n$ on $\mathscr{R}^s$. Assign the interpolant to $g\in\aleph_{\phi}(\Gamma)$ on the $(n-1)$ unisolvent set $\chi$ by $P_g$. Fix $\gamma\in \mathbb{N}_0^s$ taking $|\alpha|\leq k$. Then, there exist positive constants $\delta \nu_0$ and $C$ (independent of $\nu$, $g$ and $\Theta$) such that

where $\aleph_\phi(\Gamma)$ represents the native space of RBF and $g\in\aleph_\phi(\Gamma)$. Since the Gaussian is an infinitely smooth function, application of the above theorem yields high algebraic convergence rates. So it is concluded that for every $k\in\mathbb{N}$ and $|\alpha|\leq k$, we have

where $\mathscr{W}$ and $w$ are the respective exact and numerical solutions. Let us assume that Eq (2.20) is accurate with respect to space orders $1 < \xi_1\leq2$ and $0 < \xi_2\leq 1$, then

Let us define the residual by $e^p = \mathscr{W}^p-w^p$, then

Now, by Lax-Richtmyer definition of stability ^[18], the numerical model (2.20) is stable if

If $\mathscr{R}$ is normal, then $\|\mathscr{R}\| = \rho(\mathscr{R})$; otherwise, $\rho(\mathscr{R}) \leq \|\mathscr{R}\|$ is always satisfied. Let us assume that $\delta \nu$ and $\delta \tau$ are small enough such that the constant $c = \frac{\delta \tau^{\zeta}}{\delta \nu^({\xi_1, \xi_2})}$. Therefore, there exists a constant $\wp$ such that

Now, $e^{0} = 0$, as $e^n$ satisfies the initial condition and boundary condition. Thus, using mathematical induction, we obtain from Eq (3.4)

Using geometric series property, we obtain

Thus, the numerical model given in Eq (2.20) is convergent.

4.   Numerical results and discussion

In this part, the numerical results are discussed for numerical model (Eq (2.19)) obtained from Eq (1.1) along with initial and boundary conditions. The accuracy of all examples is analyzed by the following equations:

where $\mathscr{W}$ and $w$ are the exact and numerical solutions. The space and time convergence rates are measured by the following formulae:

Example 4.1. Set $r_1 = r_2 = r_3 = 1$ and $r_4 = -1$ in Eq (1.1) and consider $\zeta(\nu, \tau) = 1-0.5e^{-\nu\tau}$, $\xi_1(\nu, \tau) = 1.7+0.5e^{-\frac{\nu^2}{1000}-\frac{\tau}{50}-1}$ and $\xi_2(\nu, \tau) = 0.7+0.5e^{-\frac{\nu^2}{1000}-\frac{\tau}{50}-1}$, where $(\nu, \tau) = [0, 1]\times[0, 1].$ The initial and boundary conditions are

where $\mathscr{R}_1(\nu, \tau) = \nu(1-\nu)\Big(2\tau+\frac{2\tau^{1-\zeta(\nu, \tau)}}{\Gamma(3-\zeta(\nu, \tau))}\Big)$, $\mathscr{R}_2(\nu, \tau) = \tau^2\Big(\frac{\nu^{1-\xi_1(\nu, \tau)}}{\Gamma(2-\xi_1(\nu, \tau))}-\frac{2\nu^{2-\xi_1(\nu, \tau)}}{\Gamma(3-\xi_1(\nu, \tau))}\Big),$ $\mathscr{R}_3(\nu, \tau) = \tau^2\Big(\frac{\nu^{1-\xi_2(\nu, \tau)}}{\Gamma(2-\xi_2(\nu, \tau))}-\frac{2\nu^{2-\xi_2(\nu, \tau)}}{\Gamma(3-\xi_2(\nu, \tau))}\Big),$ $\mathscr{W}(\nu, \tau) = \tau^2\nu(1-\nu)$ is the exact solution. The numerical results vs exact solutions along with errors are shown in Table 1. The computed $\max(\|L\|)$ values are shown in Table 2 along with time and space convergence. The results are also shown in Figures 1–3. In Figure 1, $\delta\nu = 0.001$ and $\delta\tau = 0.001$ for 3D surface plot, while $\delta\nu = 0.1$ and $\delta\tau = 0.1$ for plot3 graph. In Figure 2, the error graphs are shown on 3D surfaces at $N = 50,500, 1000$ and $\delta\tau = 0.001$ where $\tau = 1$ while Figure 3 shows convergence rates w.r.t. space and time respectively.

Example 4.2. Put $r_1 = r_2 = 1$, $r_3 = 3$ and $r_4 = -2$ and consider $\zeta(\nu, \tau) = 0.75-0.5e^{\nu\tau}$, $\xi_1(\nu, \tau) = 1.55+0.35\cos(\pi\nu\tau)$ and $\xi_2(\nu, \tau) = 0.65+0.25\sin(\pi\nu\tau)$, where $(\nu, \tau) = [0, 1]\times[0, 1].$

The initial and boundary conditions are

where $\mathscr{R}_1(\nu, \tau) = \nu^2(1-\nu)\Big(1+\frac{\tau^{1-\zeta(\nu, \tau)}}{\Gamma(2-\zeta(\nu, \tau))}\Big)$, $\mathscr{R}_2(\nu, \tau) = \Big(\frac{2\nu^{2-\xi_1(\nu, \tau)}}{\Gamma(3-\xi_1(\nu, \tau))}-\frac{6\nu^{3-\xi_1(\nu, \tau)}}{\Gamma(4-\xi_1(\nu, \tau))}\Big)\tau$, $\mathscr{R}_3(\nu, \tau) = \Big(\frac{2\nu^{2-\xi_2(\nu, \tau)}}{\Gamma(3-\xi_2(\nu, \tau))}-\frac{6\nu^{3-\xi_2(\nu, \tau)}}{\Gamma(4-\xi_2(\nu, \tau))}\Big)\tau.$ Exact solution: $\mathscr{W}(\nu, \tau) = \nu^2(1-\nu)\tau$. The numerical approximation, exact solution and errors between exact solution and numerical approximation are shown in Table 3, while the maximum absolute errors are shown in Table 4 at different time and space intervals along with convergence rate. The results are also shown in Figures 4–6. The 3D surface plot is constructed on $\delta\nu = \delta\tau = 0.001$, while plot3 is constructed on $\delta\nu = \delta\tau = 0.1$ as shown in Figure 4. In Figure 5, the error plots are shown for $N = 100,500, 1000$ and $\delta\tau = 0.001$ at $\tau = 1$ on 3D surface, while Figure 6 shows the convergence rates w.r.t. space and time.

Example 4.3. Let us take $\zeta(\nu, \tau) = 0.55-0.25e^{-\nu\tau}$, $\xi_1(\nu, \tau) = 1.55+0.35\cos(3\pi\nu\tau)$, $\xi_2(\nu, \tau) = 0.65+0.25\sin(3\pi\nu\tau)$, $r_1 = r_2 = r_3 = 1$ and $r_4 = -1$. $w(\nu, 0) = 0$ is the initial condition while $w(0, \tau) = 0 = w(L, \tau)$ is the boundary condition. $R(\nu, \tau) = \mathscr{R}_1(\nu, \tau)-\tau\mathscr{R}_2(\nu, \tau)$, where $\mathscr{R}_1(\nu, \tau) = \Big(1+\frac{\tau^{1-\zeta(\nu, \tau)}}{\Gamma(2-\zeta(\nu, \tau))}\Big)\sin(2\pi\nu)$, $\mathscr{R}_2(\nu, \tau) = \sin(2\pi\nu+\pi^2\xi_1(\nu, \tau))-\sin(2\pi\nu+\pi^2\xi_2(\nu, \tau))$. The exact solution is $\mathscr{W}(\nu, \tau) = \tau \sin(2\pi\nu)$. All results are shown in Tables 5, 6 and Figures 7–9.

Table 7 represents the $\|L\|_\infty$ at $N = 50,100,500, 1000$ and $\delta\tau = 0.001$ vs CPU machine time. These results are computed by MATLAB 2018b software while the hardware machine is HP core i7, 11th generation processor, 16GB RAM, 2GB GPU.

5.   Conclusions

Meshfree RBF method based on the Gaussian function is used for the numerical solution of the time-space dependent order advection-dispersion mobile-immobile equation. Time derivatives are dealt with by Caputo derivative operator, while space derivatives are dealt with by Riesz, Grüwald-Letnikov derivative operators. The stability and convergence of the method are discussed and the efficiency is analyzed by the maximum norm. The method is tested through some examples and the results are shown in tables and figures. In the future, based on these results, it will be more interesting to see how the proposed method can work on such PDE models defined on irregular domains, complex function orders, and in case when the exact solution of the PDE model is unknown.

Conflict of interest

The authors have no conflicts of interest.

Acknowledgments

The work was supported by Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2023R8), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.