An effective approach based on Hybrid B-spline to solve Riesz space fractional partial differential equations

1.   Introduction

In many areas of the modern life, mathematical modeling and numerical outcomes are integrated. Most of the mathematical methods, which can be used in engineering and natural sciences are based on ordinary differential equations (ODEs), partial differential equations (PDEs) and integral equations (IEs). Many fields of engineering and sciences including probability theory, fluid flow and electromagnetic theory represent the advantages of fractional equations over the ordinary equations ^[1,2,3].

Fractional differential equations have rich physical background in different fields of sciences. Fractional derivative provides an adequate approach to describe the physical models, which preserve the memory and hereditary properties of physical substances. This versatility of fractional derivative makes it more applicable in the field of physics ^[4], biology ^[5], hydrology ^[6] and chemistry ^[7]. More expressly, they are used in dynamical systems with quasi-chaotic behaviour, chaotic dynamical behavior, porous media or dynamics of complex material. Now-a-days, the interest of some of the scholars have been demonstrated in the form of research on the study of mathematical theory of fractional differential equations ^[8,9,10].

Among various definitions of fractional derivative, the Riesz fractional derivative involves a sum of the Riemann-Liouville right- and left- sided derivatives. This induction facilitates to describe the impact of fluid flow on both sides of the domain. The Riesz space fractional PDE was derived by Saichev and Zaslavsky ^[11] from the kinetics of chaotic dynamics. The kinetic equation with the fractional derivatives are considered for the non-negative values of $K_\alpha, K_\beta$ in order to observe the flow from left to right. The discretized version of the Riesz fractional derivative in space can facilitate the solution of inconsistency problem. One of the general form of the Riesz space fractional diffusion equation (RSF-DE) ^[12] for $1 < \alpha\leqslant2$ is

with boundary conditions $G(a, t) = G(b, t) = 0$, initial condition $G(x, 0) = \psi(x)$ and here $K_\alpha$ is diffusion coefficient.

Let us consider the Riesz space fractional advection dispersion equation (RSF-ADE) ^[12] for $1 < \alpha \leqslant 2$ and $0 < \beta < 1$ as

with boundary conditions $G(a, t) = G(b, t) = 0$ and initial condition $G(x, 0) = \psi(x)$, where $\alpha, \beta$ are fractional parameters, $K_{\alpha}, K_{\beta}$ represent the dispersion coefficient and average fluid velocity, $G(x, t)$ is an unknown function of space and time, $\psi(x)$ is a continuous functions and $f(x, t)$ is a source term.

The Riesz fractional operator on a finite interval $0\leqslant x\leqslant L$ for $n-1 < \alpha\leqslant n$ ($n$ is positive integer and $L$ is a real number) is described as

where

Here $x_i = a+ih, i = 0, 1, 2, \ldots, N, h = \frac{b-a}{N}$, $\; \; G_{m+p-1}^{n+1} = G(x_{m+p-1}, t_{n+1})$, $w_{0}^{\nu} = 1, w_{p}^{\nu} = \frac{(-1)^\nu \nu(\nu-1)\cdots (\nu-p+1)}{p!}$, $p = 1, 2, \ldots$, $\nu\in\{\alpha, \beta\}.$ Moreover ${}_{0}D_{x}^{\nu}$ and ${}_{x}D_{L}^{\nu}$ for $\nu \in \{\alpha, \beta\}$ are Riemann-Liouville left and right sided derivatives ^[12].

Several numerical techniques have been used to solve RSF-PDEs, Shen et al. ^[13] have worked over the fundamental solution of RSF-ADE using Laplace and Fourier transforms methods. Zhang et al. ^[14] have proposed a backward difference scheme and a finite element method to solve a symmetric Space fractional PDE with the Riesz fractional operator and Riemann-Liouville discretization. Ding et al. ^[15] have introduced a new numerical method, in which they proved its effectiveness for Reisz problem. Liao et al. ^[16] has used second order BDF approximation for the solution of RSF-DE, a series of fast parallel preconditioning methods based on implicit technique for the solution of said equation has been utilized by Gu et al. ^[17,18]. Jian et al. ^[19,20,21] used as fast iterative methods for the solution of fractional reaction diffusion and convection reaction diffusion equation.

From the existing results, it has been seen that spline functions have been used to solve DEs. Zahra and Elkholy ^[22] presented a numerical solution of FDEs by the use of cubic splines. Non classical diffusion problems by Abbas et al. ^[23], advection-diffusion problems by Nazir et al. ^[24] and fractional sub-diffusion equation by Zhu et al. ^[25] are solved by using B-spline collocation method. Yaseen et al. ^[26] presented a finite difference scheme based on cubic trigonometric B-splines for a time fractional diffusion-wave equation. Hashmi et al. ^[27] has solved Hunter-Saxton equation by using cubic trigonometric B-spline collocation method. Most recently, B-spline based collocation method ^[28,29] and differential quadrature method have been applied by Hashmi et al. ^[30,31] for the solution of fractional PDEs.

From the existing results, it is evident that no attempt has been made to solve RSF-PDEs using HBCM. In present research, formulation of Riesz space fractional derivatives using the Riemann Liouville right and left sided derivatives, indebted with HBCM has been used. Von Neumann stability analysis approach is used to achieve the convergence, which shows that our proposed method is unconditionally stable. Numerical examples with their tabular and graphical representations are included to illustrate the efficiency of proposed method.

The layout of this paper is as follows: HBCM is given in Section $2$. In Section $3$, temporal discretization, proposed methodology and convergence analysis of RSF-DE is presented. In Section $4$ temporal discretization of RSF-ADE along with proposed method is presented and demonstrated with examples. The conclusion of the proposed technique for both models is discussed in last section.

2.   Hybrid B-spline collocation method

Let $\tau$ be the time step such that $t_n = n \tau$, $n = 0, 1, 2, \ldots, M,$ i.e. $0 \leqslant t_n \leqslant T$, and $h$ is the step size with $x_m = a+mh, a \leqslant x_m \leqslant b, m = 0, 1, 2, \ldots, N$. Let $G_m^n$ be the numerical estimate of $G(x, t)$ at the mesh point $(x_m, t_n)$. Now HBCM is used to solve the FPDEs. The approximate solution of $G(x, t)$ is considered as

where $C_j(t)$ are unknown parameters and $H_{3, j}(x)$ is hybrid cubic B-spline basis function of order three and described as

where $\gamma = 1$ and $\gamma = 0$, represent the cubic B-spline basis $B_{3, j}(x)$ and trigonometric B-spline basis $T_{3, j}(x)$ function respectively. For Hybrid bases, Table 1 represents the functional values of the hybrid B-spline: Here

and

The B-spline representation of $G(x, t)$ with its derivatives at $(x_m, t_{n+1})$ have the form

3.   Riesz space fractional diffusion equation (RSF-DE)

FPDEs are undergoing rapid development due to their wide application in engineering, economics and physics. In the present research, RSF-DE is solved by using HBCM. The proposed methodology is as follows:

3.1.   Proposed methodology

The time derivative and the Riesz fractional derivatives are discretized by implicit Euler scheme and the shifted Grunwald formulae respectively, in order to apply HCBM on (1.1) as follows:

simplifying (3.1) with the help of (2.3), the resultant equation becomes

where $\eta = \frac{K_\alpha \tau}{2\cos(\frac{\pi\alpha}{2}){h^\alpha}}\; {\rm{and}} \; w_0^\alpha = 0$. For $(m = 0, \ldots, N)$, the above equation provides a rectangular system of unknowns. So the BCs are used to obtained two more additional constraints for finding a unique solution of the system described, which are given below

Now the obtained system can be solved by a variant of Thomas algorithm.

3.2.   Initial state

First, we will find an initial vector to start the iterative process i.e. $C^0 = [C_{-1}^{0}, C_{0}^{0}, \ldots, C_{N+1}^{0}]^T$. This can be determined by utilizing the IC on the approximate solution

where $C_{m}(0)'s$ are unknown parameters. Here $G(x, 0)$ is required, which satisfies the following conditions:

This can be written as $(N+3)\times (N+3)$ matrix system:

where $C^0 = [C_{-1}^{0}, C_{0}^{0}, \ldots, C_{N+1}^{0}]^t$,

and $B = [{\psi}^{'}(x_{0}), {\psi}(x_{0}), \ldots, {\psi}(x_{N}), {\psi}^{'}(x_{N})]^t.$ This will provide us a starting solution to the iterative process given in (3.2)–(3.4).

3.3.   Stability analysis

The stability analysis is performed using Von Neumann stability analysis in terms of evaluation of amplitude of the error. For this, we considered the homogenous problem, it is as follows:

where $m$ and $\phi$ are real numbers. By substituting $(3.8)$ in $(3.1)$ and simplifying, we have

The Eq (3.9) can be re-written as

where

Now the relation has the form

It concludes that the system is unconditionally stable.

3.4.   Convergence of the scheme

The convergence estimates for the proposed scheme is presented here for the homogenous case:

Theorem 3.1. Let $\{G(x, t^n)\}_{n = 0}^{M}$ be the exact solution of the problem and let $\{G^n\}$ be the approximate solution of proposed problem then

where $e^{n+1} = \{|G(x_{m}, t^{n+1})-G_{m}^{n+1}|; \; \; m = 0, 1, \ldots, N\}$ and $D = \left(1+\eta \sum\limits_{p = 0}^{m+1}w_p^{\alpha} +\eta \sum\limits_{p = 0}^{N-m+1}w_p^{\alpha} \right).$

Proof. For $h > 0$, we know that $R_1 = \frac{\gamma}{6}+(1-\gamma)r_1, R_{2} = \frac{{2\gamma}}{3}+(1-\gamma)r_2,$ then:

● For $\gamma = 1$ implies $R_1 = \frac{1}{6}, R_2 = \frac{2}{3}$ implies $R_2 > R_1$

● For $\gamma = 0$ implies $R_1 = \frac{\sin^{2}(\frac{h}{2})}{\sin(h)\sin(\frac{3h}{2})}, R_2 = \frac{2}{1+2\cos(h)}$ implies $R_2 > R_1$ for all $h > 0$ in the domain of the problem. This results that the coefficient matrix $A$ given in (3.7) is strictly diagonally dominant and will produce a unique solution.

For the proposed problem, let the exact solution satisfies the numerical scheme (3.1), then we have

or

and then the approximate solution will have the form

subtracting above two equations gives

where $e^{n+1} = \max \{|G(x_{m}, t^{n+1})-G_{m}^{n+1}|; \; \; m = 0, 1, \ldots, N\}$. Taking the inner product with $e^{n+1}$ on both sides and using the fact that $\langle y, y\rangle = \|y\|^2$, we have

Dividing it by $\|e^{n+1}\|$ gives

which means

It is clear that $D\geqslant 1$, hence the proposed scheme is convergent.

3.5.   Example 1

In this example, HBCM is applied on RSF-DE to showcase the efficiency of proposed method. Consider an example of the form ^[12]

with BCs $G(0, t) = 0, G(\pi, t) = 0$ and the IC $G(x, 0) = x^2(\pi-x)$. The exact solution is $G(x, t) = \sum\limits_{n = 1}^{\infty}b_n sin(n x)e^{(-(n^2)^\frac{\alpha}{2}K_\alpha t)}$ where $b_n = \frac{8}{n^3}(-1)^{n+1}-\frac{4}{n^3}$.

Table 2 shows the absolute errors at different values of $\gamma$, which indicates that the results are acceptable for all choices of $\gamma$. Table 3 at $N = 2^4$ and Table 4 at $N = 2^6$ present the numerical display of absolute errors of the RSF-DE for different values of $\alpha$. Similarly Table 5 describes the numerical display of errors for $N = 2^8$ at various values of $\alpha$. From these taables, we concludes that proposed scheme has more convergent solutions for increasing values of N. Table 6 is providing the comparison between the spline collocation method and the solution obtained by ^[12]. In Figure 1(a–d), the geometrical display of absolute errors at different value of $\alpha$ has been illustrated, which indicates that the higher values of $\alpha$ provides the improved solution.

4.   Riesz space fractional advection-dispersion equation (RSF-ADE)

RSF-ADE is used in groundwater hydrology to model the transport of passive tracers carried by fluid flow in a porous medium. The space RSF-ADE have been used to describe important phenomena such as colloid, comb structures, bio-logical systems, glassy, disordered and porous media. In current section, the proposed methodology, stability analysis and numerical application of HBCM on RSF-ADE is performed.

4.1.   Proposed methodology

Time derivative and RSF are discretized by implicit Euler scheme and shifted Grunwald formulae respectively in order to apply HBCM on (1.2) as follows:

re-arranging and simplifying we have

By substituting (2.3) in the above equation, we get

where

For $(m = 0, \ldots, N)$, the above equation provides a rectangular system of $\left\{C^{n+1}_{-1}, C^{n+1}_{0}, \right.$ $\left.C^{n+1}_{1}, \ldots, C^{n+1}_{N+1}\right\}$. The BCs are used to obtained two more constraints for finding a unique solution of the system, which are given below:

The initial vector can be obtained by the procedure as described in pervious section.

4.2.   Stability analysis

From the methodology of previous section, we can proceed by substituting

in homogenous case of (4.3), which will have the form

Simplifying the above relation, we get

which can also be written as

where

Hence

So, we get $\xi_{n+1} < |\xi_{n}|$ because $\mu = (1+C_1+C_2+C_3+C_4) > 1$.

This concludes that the system is unconditionally convergent.

4.3.   Convergence of the scheme

The convergence estimates for the proposed scheme is presented here for the homogenous case:

Theorem 4.1. Let $\{G(x, t^n)\}_{n = 0}^{M}$ be the exact solution of the problem and let $\{G^n\}$ be the approximate solution of proposed problem, then

where $e^{n+1} = \{|G(x_{m}, t^{n+1})-G_{m}^{n+1}|; \; \; m = 0, 1, \ldots, N\}$ and

Proof. For the proposed problem, let the exact solution satisfies the numerical scheme (4.2), then we have

or

and

Subtracting above two equations gives

where $e^{n+1} = \max \{|G(x_{m}, t^{n+1})-G_{m}^{n+1}|; \; \; m = 0, 1, \ldots, N\}$. Taking the inner product with $e^{n+1}$ on both sides and using the fact that $\langle y, y\rangle = \|y\|^2$, we have

Dividing by $\|e^{n+1}\|$ gives

which means

It is clear that $D\geqslant 1$, hence the proposed scheme is convergent.

4.4.   Example 2

Consider the following example ^[12]

with BCs and IC

where

The exact solution is $G(x, t) = t^\alpha e^{\beta t}x^2(1-x)^2.$ Table 7 shows the absolute errors at different values of $\gamma$, which indicates that the results are nearly same for different values of $\gamma$. Table 8 for $N = 40$ shows the numerical display of absolute errors of the RSF-ADE for different values of $\alpha$ in example $2$. Similarly Table 9 describes the numerical display of absolute errors for $N = 40$ and $T = 0.001$ at different values of $\alpha$. From these tables, we see that the solution is better for smaller values of $T$. Table 10 is providing the comparison between the maximum error of the approximate solution and solution in ^[12].

Figure 2(a) shows the geometrical representation of error for $N = 2^6$, $T = 0.001$ at different values of $\alpha$, which indicates that by increasing the value of $\alpha$ the approximate solution move closer to exact one. Geometrical representation of both solutions for $N = 2^6$ is shown in Figure 2(b) which indicates that exact and approximate solution are very much coinciding.

Algorithm: Hybrid B-spline collocation method for RSF-DE and RSF-ADE

$(1)$ Discretize $x_i, t_j$ for $i = 0, 1, \ldots, N, j = 0, 1, \ldots, M$.

$(2)$ Evaluate $w_p^\alpha, \; \; p = 0, 1, \ldots, N$ for RSF-DE ($w_p^\alpha, w_p^\beta, \; \; p = 0, 1, \ldots, N$ for RSF-ADE).

$(3)$ Initialize the vector $C^0 = [C_{-1}^0, C_{0}^0, \ldots, C_{N+1}^0]$ using initial conditions and Thomas algorithm.

$(4)$ Discretize $G(x_m, t_{n+1})$ and replace by $R_{1}C_{m-1}^{n+1}+R_{2}C_{m}^{n+1}+R_{1}C_{m+1}^{n+1}$ for $k = 0, 1, \ldots, M-1$.

$(5)$ For $i = 1, 2, \ldots, N-1$

(a) Construction of square system of equations using (3.2) for RSF-DE and (4.2) for RSF-ADE,

(b) Using fsolve evaluate the solution at each time level till final time.

(6) Substitute in B-spline basis function (2.1).

5.   Conclusions and future work

The aim of this research was to solve RSF-DE and RSF-ADE using HBCM. We have used the left and right sided Riemann Liouville formulae to discretize the Riesz fractional derivative. The discretized equation is then updated with the spanning form of hybrid B-spline basis. The convergence of numerical scheme is ensured by using VN stability analysis. The methods are employed on numerical examples of both models to show the applicability of prescribed method. The results are displayed using different tables, which shows that solution becomes convergent by decreasing the values of $T$. Numerical results are also shown graphically, which indicate that the approximate solutions are very much coinciding with the exact solution. The comparison provided in Tables 6 and 10 with ^[12] indicate that approximate solution are comparable. The display of CPU time and Memory allocation is placed in Table 11. However, the proposed numerical algorithms can be extended for tempered fractional diffusion equation ^[32] with some adjustments.

Acknowledgments

This work is supported by the Higher Education Commission of Pakistan under the NRPU project grant number 9306/Punjab/NRPU/R & D/HEC/2017.

Conflict of interest

The authors declare no conflicts of interest.