Numerical study for a class of time fractional diffusion equations using operational matrices based on Hosoya polynomial

1.   Introduction

Non-linear partial differential equations are extensively used in science and engineering to model real-world phenomena ^[1,2,3,4]. Using fractional operators like the Riemann-Liouville (RL) and the Caputo operators which have local and singular kernels, it is difficult to express many non-local dynamics systems. Thus to describe complex physical problems, fractional operators with non-local and non-singular kernels ^[5,6] were defined. The Atangana-Baleanu (AB) fractional derivative operator is one of these type of fractional operators which is introduced by Atangana and Baleanu^[7].

The time fractional Kolmogorov equations (TF-KEs) are defined as

with the initial and boundary conditions

where $(s, t)\in [0, 1]\times [0, 1]$, ${^{AB}\! D}_{t}^{\gamma}$ denotes the Atangana-Baleanu (AB) derivative operator, $D_{s} g(s, t) = \dfrac{\partial}{\partial s}g(s, t)$ and $D_{ss} g(s, t) = \dfrac{\partial^{2}}{\partial s^{2}}g(s, t)$. If $\vartheta_{1}(s)$ and $\vartheta_{2}(s)$ are constants, then Eq (1.1) is presenting the time fractional advection-diffusion equations (TF-ADEs).

Many researchers are developing methods to find the solution of partial differential equations of fractional order. Analytical solutions or formal solutions of such type of equations are difficult; therefore, numerical simulations of these equations inspire a large amount of attentions. High accuracy methods can illustrate the anomalous diffusion phenomenon more precisely. Some of the efficient techniques are Adomian decomposition ^[8,9], a two-grid temporal second-order scheme ^[10], the Galerkin finite element method ^[11], finite difference ^[12], a differential transform ^[13], the orthogonal spline collocation method ^[14], the optimal homotopy asymptotic method ^[15], an operational matrix (OM) ^{[16,17,18,19,20,21,22,23,24]}, etc.

The OM is one of the numerical tools to find the solution of a variety of differential equations. OMs of fractional derivatives and integration were derived using polynomials like the Chebyshev ^[16], Legendre ^[17,18], Bernstein ^[19], clique ^[20], Genocchi ^[21], Bernoulli ^[22], etc. In this work, with the help of the Hosoya polynomial (HS) of simple paths and OMs, we reduce problem (1.1) to the solution of a system of nonlinear algebraic equations, which greatly simplifies the problem under study.

The sections are arranged as follows. In Section 2, we review some basic preliminaries in fractional calculus and interesting properties of the HP. Section 3 presents a new technique to solve the TF-KEs. The efficiency and simplicity of the proposed method using examples are discussed in Section 5. In Section 6, the conclusion is given.

2.   Preliminaries

In this section we discuss some basic preliminaries of fractional calculus and the main properties of the HP. We also compute an error bound for the numerical solution.

2.1.   Fractional calculus

Definition 2.1. (See ^[25]) Let $0 < \gamma\leq 1$. The RL integral of order $\gamma$ is defined as

One of the properties of the fractional order of RL integral is

Definition 2.2. (See ^[7]) Let $0 < \gamma \leq 1$, $g\in H^{1}(0, 1)$ and $\Phi(\gamma)$ be a normalization function such that $\Phi(0) = \Phi(1) = 1$ and $\Phi(\gamma) = 1-\gamma+\dfrac{\gamma}{\Gamma(\gamma)}$. Then, the following holds

1) The AB derivative is defined as

where $E_{\gamma}(s) = \sum\limits_{j = 0}^{\infty}\frac{s^{j}}{\Gamma(\gamma j+1)}$ is the Mittag-Leffler function.

2) The AB integral is given as

Let $v_{\gamma} = \dfrac{1-\gamma}{\Phi(\gamma)}$ and $w_{\gamma} = \dfrac{1}{\Phi(\gamma)\Gamma(\gamma)}$; then, we can rewrite (2.1) as

The AB integral satisfies the following property ^[26]:

2.2.   The HP and their properties

In 1988, Haruo Hosoya introduced the concept of the HP ^[27,28]. This polynomial is used to calculate distance between vertices of a graph ^[29]. In ^[30,31], the HP of path graphs is obtained. The HP of the path graphs is described as

where $d(G, l)$ denotes the distance between vertex pairs in the path graph ^[32,33]. Here we consider path graph with vertices $n$ where $n\in N$. Based on $n$ vertex values the Hosoya polynomials are calculated ^[34]. Let us consider the path $P_{n}$ with $n$ vertices; then the HP of the $P_{i}, \, i = 1, 2, \cdots, n$ are computed as

Consider any function $g(s)$ in $L^{2}(0, 1)$; we can approximate it using the HP as follows:

where

and

From (2.2), we have

where ${\mathcal{Q}} = \langle\mathrm{H}(s), \mathrm{H}(s)\rangle$ and $\langle \cdot, \cdot\rangle$ denotes the inner product of two arbitrary functions.

Now, consider the function $g(s, t) \in L^{2}\left([0, 1]\times [0, 1]\right)$; then, it can be expanded in terms of the HP by using the infinite series,

If we consider the first $(N+1)^{2}$ terms in (2.4), an approximation of the function $g(s, t)$ is obtained as

where

Theorem 2.1. The integral of the vector $\mathrm{H}(s)$ given by (2.3) can be approximated as

where $R$ is called the OM of integration for the HP.

Proof. Firstly, we express the basis vector of the HP, $\mathrm{H}(s)$, in terms of the Taylor basis functions,

where

and

with

Now, we can write

where $\mathcal{B} = [b_{q, r}], \; q, r = 1, 2, \cdots, N+1$ is an $(N+1)\times(N+1)$ matrix with the following elements

and

Now, by approximating $s^{k}, \; k = 1, 2, \cdots, N+1$ in terms of the HP and by (2.7), we have

where $\mathcal{A}^{-1}_{r}$, $r = 2, 3, \cdots, N+1$ is the $r$-th row of the matrix $\mathcal{A}^{-1}$ and $\mathcal{L} = {\mathcal{Q}}^{-1}\langle s^{N+1}, \mathrm{H}(s)\rangle$. Then, we get

where $\mathcal{E} = \begin{bmatrix}\mathcal{A}_{2}^{-1}, \mathcal{A}_{3}^{-1}, \cdots, \mathcal{A}_{N+1}^{-1}, \mathcal{L}^{T} \end{bmatrix}^{T}$. Therefore, by taking $R = \mathcal{A}\mathcal{B}\mathcal{E}$, the proof is completed.

Theorem 2.2. The OM of the product based on the HP is given by (2.3) can be approximated as

where $\widehat{C}$ is called the OM of product for the HP.

Proof. Multiplying the vector $C = [c_1, c_2, \cdots, c_{N+1}]^T$ by $\mathrm{H}(s)$ and $\mathrm{H}^{T}(s)$ gives

Taking $e_{k, i} = \left[e_{k, i}^{1}, e_{k, i}^{2}, \cdots, e_{k, i}^{N+1} \right]^{T}$ and expanding $s^{k-1}\, \tilde{ H}(P_{i}, s)\simeq e_{k, i}^{T}\, \mathrm{H}(s), \; i, k = 1, 2, \cdots, N+1$ using the HP, we can write

Therefore,

where $E_{k}$ is an $(N+1)\times(N+1)$ matrix and the vectors $e_{k, i}$ for $k = 1, 2, \cdots, N+1$ are the columns of $E_{k}$. Let $\overline{E}_{k} = E_{k}\, C, \quad k = 1, 2, \cdots, N+1$. Setting $\overline{C} = \left[\overline{E}_{1}, \overline{E}_{2}, \cdots, \overline{E}_{N+1} \right]$ as an $(N+1)\times(N+1)$ matrix and using (2.8) and (2.9), we have

where by taking $\widehat{C} = \overline{C}\, A^{T}$, the proof is completed.

Theorem 2.3. Consider the given vector $\mathrm{H}(s)$ in (2.3); the fractional RL integral of this vector is approximated as

where $\mathcal{P}^{\gamma}$ is named the OM based on the HP which is given by

with

Proof. First, we rewrite $\tilde{H}(P_{i}, s)$ in the following form:

Let us apply, the RL integral operator, ${^{R\!L}\!I^{\gamma}_s }$, on $\tilde{H}(P_{i}, s), \; i = 1, \cdots, N+1$; this yields

Now, using the HP, the function $s^{k+\gamma-1}$ is approximated as:

By substituting (2.11) into (2.10), we have,

Theorem 2.4. Suppose that $0 < \gamma\leq 1$ and $\tilde{H}(P_{i}, x)$ is the HP vector; then,

where $\mathcal{I}^{\gamma} = v_{\gamma} I+w_{\gamma}\Gamma(\gamma+1) \mathcal{P}^{\gamma}$ is called the OM of the AB-integral based on the HP and $I$ is an $(N+1)\times(N+1)$ identity matrix.

Proof. Applying the AB integral operator, ${^{AB}\! I_{s}^{\gamma}}$, on $\mathrm{H}(s)$ yields

According to Theorem 2.3, we have that ${^{R\!L}\! I_{s}^{\gamma}}\mathrm{H}(s)\simeq \mathcal{P}^{\gamma}\mathrm{H}(s)$. Therefore

Setting $\mathcal{I}^{\gamma} = v_{\gamma}I+w_{\gamma}\Gamma(\gamma+1)\mathcal{P}^{\gamma}$, the proof is complete.

3.   The proposed technique

The main aim of this section is to introduce a technique based on the HP of simple paths to find the solution of the TF-KEs. To do this, we first expand $D_{ss} g(s, t)$ as

Integrating (3.1) with respect to $s$ gives

Again integrating the above equation with respect to $s$ gives

By putting $s = 1$ into (3.3), we have

By substituting (3.4) into (3.3), we get

Now, we approximate that $d_{1}(t) = S_{0}^{T}\mathrm{H}(t), d_{2}(t) = S_{1}^{T}\mathrm{H}(t)$ and $s = \mathrm{H}^{T}(s)S$ and putting in (3.5), we get

The above relation can be written as

Approximating $1 = \widehat{S}^{T}\mathrm{H}(s) = \mathrm{H}^{T}(s)\widehat{S}$, the above relation is rewritten as

Setting $\rho_{1} = \widehat{S}S_{0}^{T}+S S_{1}^{T}-SS_{0}^{T}-S \mathrm{H}^{T}(1)\left(R^{2}\right)^{T} \tilde{ h}+\left(R^{2}\right)^{T} \tilde{ h}$, we have

According to (1.1), we need to obtain $D_{s}\ g(s, t)$. Putting the approximations $d_{1}(t), d_{2}(t)$ and the relation (3.4) into (3.2) yields

The above relation can be written as

Putting $1 = H^{T}(s)\widehat{S}$ into the above relation, we get

Setting $\rho_{2} = \widehat{S} S_{1}^{T}-\widehat{S} S_{0}^{T}- \widehat{S}\mathrm{H}^{T}(1)\left(R^{2}\right)^{T} \tilde{ h}+R^{T} \tilde{ h}$, we have

Applying ${^{AB}\! I_{t}^{\gamma}}$ to (1.1), putting $g(s, t)\simeq \mathrm{H}^{T}(s)\rho_{1}\mathrm{H}(t), D_{s}g(s, t)\simeq \mathrm{H}^{T}(s)\rho_{2}\mathrm{H}(t)$, $D_{ss}\ g(s, t)\simeq\mathrm{H}^{T}(s) \tilde{ h}\mathrm{H}(t)$ and approximating $\omega\left(s, t\right)\simeq\mathrm{H}^{T}(s)\rho_{3}\mathrm{H}(t)$ in (1.1) yields

Now approximating $d_{0}(s)\simeq \mathrm{H}^{T}(s)S_{2}, \vartheta_{1}(s) \simeq S_{3}^{T}\mathrm{H}(s), \vartheta_{2}(s)\simeq S_{4}^{T}\mathrm{H}(s)$ and using Theorem 2.4, the above relation can be rewritten as

By Theorem 2.2, the above relation can be written as

Now approximating $1 = \widehat{S}^{T}\mathrm{H}(t)$, we have

We can write the above relation as

Therefore we have

By solving the obtained system, we find $h_{ij}$, $i, j = 1, 2, \cdots, N+1$. Consequently, $g(s, t)$ can be calculated by using (3.7).

4.   Convergence analysis

Set $I = (a, b)^{n}, \, n = 2, 3$ in $\mathbb{R}^{n}$. The Sobolev norm is given as

where $D_{l}^{(k)}u$ and $H^\epsilon (I)$ are the $k$-th derivative of $g$ and Sobolev space, respectively. The notation $\vert g\vert_{H^{\epsilon; N}}$ is given as ^[35]

Theorem 4.1 (See ^[36]). Let $g(s, t)\in H^{\epsilon}(I)$ with $\epsilon\geq 1$. Considering $P_{N}g(s, t) = \sum\limits_{r = 1}^{N+1} \sum\limits_{n = 1}^{N+1}a_{r, n}P_{r}(s) P_{n}(t)$ as the best approximation of $g(s, t)$, we have

and if $1\leq \iota\leq \epsilon$, then

with

Lemma 4.1. The AB derivative can be written by using the fractional order RL integral as follows:

Proof. According to the definitions of the AB derivative and the RL integral, the proof is complete.

Theorem 4.2. Suppose that $0 < \gamma\leq 1, \, \left|\vartheta_{1}(s)\right|\leq \tau_{1}, \, \left|\vartheta_{2}(s)\right|\leq \tau_{2}$ and $g(s, t)\in H^{\epsilon}(I)$ with $\epsilon\geq 1$. If $E(s, t)$ is the residual error by approximating $g(s, t)$, then $E(s, t)$ can be evaluated as

where $1\leq \iota\leq \epsilon$ and $\varrho_{1}$ is a constant number.

Proof. According to (1.1),

and

Substituting Eqs (4.1) and (4.2) in $E(s, t)$ yields

and then

Now, we must find a bound for $\|{^{AB}\! D}_{t}^{\gamma}(g(s, t)-g_{N}(s, t))\|_{L^{2}(I)}$. In view of ^[26], and by using Lemma 4.1, in a similar way, we write

Therefore,

where $\frac{\Phi(\gamma)}{1-\gamma}E_{\gamma, 2}(\varpi)\leq \delta_{1}$. Thus, from (4.4), we can write

where $|g|^*_{H^{\epsilon; N}(I)} = C N^{\vartheta(\iota)-\epsilon}|g|_{H^{\epsilon; N}(I)}$. By Theorem 4.1,

and

where $|D_{s}g|^*_{H^{\epsilon; N}(I)} = C N^{\vartheta(\iota)-\epsilon}|D_{s}g|_{H^{\epsilon; N}(I)}$. Taking $\varrho_{1} = \max\{\delta_{1}+ \tau_{1}, \tau_{2}\}$ and substituting (4.5)–(4.7) into (4.3); then, the desired result is obtained.

5.   Test problems

In this section, the proposed technique which is described in Section 3 is shown to be tested using some numerical examples. The codes are written in Mathematica software.

Example 5.1. Consider (1.1) with $\vartheta_{1}(s) = -1, \vartheta_{2}(s) = 0.1$ and $\omega\left(s, t\right) = 0$. The initial and boundary conditions can be extracted from the analytical solution $g(s, t) = \tau_{0}e^{\tau_{1}t-\tau_{2}s}$ when $\gamma = 1$. Setting $\tau_{0} = 1, \tau_{1} = 0.2, \tau_{2} = \dfrac{\vartheta_{1}(s)+\sqrt{\vartheta_{1}^{2}(s)+4\vartheta_{2}(s)\tau_1}}{2\vartheta_{2}(s)}$, considering $N = 3$ and using the proposed technique, the numerical results of the TF-ADE are reported in Tables 1 and 2, and in Figures 1–3.

Example 5.2. Consider (1.1) with $\vartheta_{1}(s) = s, \vartheta_{2}(s) = \dfrac{s^{2}}{2}$ and $\omega\left(s, t\right) = 0$. The initial and boundary conditions can be extracted from the analytical solution $g(s, t) = s E_{\alpha}(t^{\alpha})$. By setting $N = 5$ and using the proposed technique, the numerical results of the TF–KE are as reported in Figures 4–6.

6.   Conclusions

Time fractional Kolmogorov equations and time fractional advection-diffusion equations have been used to model many problems in mathematical physics and many scientific applications. Developing efficient methods for solving such equations plays an important role. In this paper, a proposed technique is used to solve TF-ADEs and TF-KEs. This technique reduces the problems under study to a set of algebraic equations. Then, solving the system of equations will give the numerical solution. An error estimate is provided. This method was tested on a few examples of TF-ADEs and TF-KEs to check the accuracy and applicability. This method might be applied for system of fractional order integro-differential equations and partial differential equations as well.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

The authors would like to thank for the support from Scientific Research Fund Project of Yunnan Provincial Department of Education, No. 2022J0949. The authors also would like to thank the anonymous reviewers for their valuable and constructive comments to improve our paper.

Conflict of interest

The authors declare there is no conflicts of interest.