Dynamics of a general model of nonlinear difference equations and its applications to LPA model

1.   Introduction

Fractional differential equations have been applied in modelling of many fields of physics and processes such as earthquakes, optics, finance, hydrology, traffic flow, fluid mechanics, fractional kinetics, mathematical biology, measurement of visco-elastic material properties, electrical network, electro-chemistry electro-magnetic, signal processing, control theory, acoustics, material sciences ^{[2,5,8,9,17,22,39,41,42,43,44,45,46]}. The nonlocal property is the most important advantage of these equations showing state of a complex system does not depend only on its current state but also depends on its all-previous states. Due to this reason, the fractional calculus has becoming more and more popular.

The fractional diffusion equations (FDEs) are generalized forms of classical diffusion equations. These equations have been solved by many researchers. For example, authors of ^{[6,11,24,26,32]} have been applied finite element method, compact FDM, Crank-Nicolson FDM, B-spline based method, implicit difference approximation, respectively, to solve FDEs. Khader ^[14] considered an efficient method based upon Chebyshev approximations. Murio ^[20] developed an implicit unconditionally stable method. Sun et al. ^[25] applied a semi-analytical FEM to solve a class of these equations. Tadjeran et al. ^[27] applied aforesaid method with spatial extrapolation to solve a class of variable coefficients FDEs. Lin and Xu ^[15] constructed a method based on FDM and Legendre spectral method while Ghanbari and Atangana ^[7] presented a method based on the product-integration rule for solving aforesaid equations. Çelik and Duman ^[3] examined a Crank-Nicolson method with the Riesz fractional derivative. Zhai and Feng ^[35] introduced a block-centered FDM for solving these equations.

Additionally, Zhai et al. ^[34,36] constructed unconditionally stable FDMs to solve the time-space fractional diffusion and three-dimensional time-fractional subdiffusion equations, respectively. Wu and Zhai ^[33] proposed a high order FDM to solve the 2D time-fractional convection-dominated diffusion equation while Zhai et al. ^[37,38] proposed a high-order compact FDM and ADI method to solve the 3D fractional convection-diffusion equation. Furthermore, Hanert ^[10] proposed a flexible numerical method to discretize the space-time FDEs, where pseudo-spectral expansion has been used for discretization of the time derivative and either high order pseudo-spectral or low order FE and FD have been used for space derivative. Murillo and Yuste ^[19] considered explicit difference scheme, where three-point centered formula has been used to approximate the spatial derivative. Verma et al. ^[31] studied nonlinear diffusion equations analytically and numerically using the classical Lie symmetry method.

In this paper, we consider a Caputo time-fractional diffusion equation

subject to the initial condition

the Dirichlet boundary conditions

where ${}^{c} D_{t}^{\alpha} u\left(x, t_n \right)$ denotes the Caputo time-fractional derivative ^[40] as

2.   Discretization of the problem

First, the Caputo time-fractional derivative is discretized. For this, we partition the time domain $[0, T]$ uniformly as $0 = t_0 < t_1, ..., < t_N = T$, where interval $\Delta t = t_n-t_{n-1} = \frac{T}{N}$ for $n = 1, 2, \dots, N$. Here N is the number of time mesh. By the definition of Caputo derivative at time $t = t_n$, the time-fractional derivative is approximated as

where $\bar{b}_l = {\left(l+1\right)}^{1-\alpha }-l^{1-\alpha }\ \forall \ 0\le l \le N$ and ${T_r}^n_{\Delta t}$ gives the truncation error $\le \bar{c}_u\Delta t^{2-\alpha }$, where the constant $\bar{c}_u$ is only related to $u(x, t)$.

Lemma 3.1 The below relations should be hold for the coefficients $\bar{b}_l$

$(1)$$\bar{b}_0 = 1$,

$(2)$$\bar{b}_l > 0, \quad \forall \ \ 0\le l\le N,$

$(3)$$\bar{b}_{l-1} > \bar{b}_l, \ \forall \ 1\le l\le N$.

Proof. For the proof of Lemma 3.1, see ^[16].

Next, we use CExpB-spline collocation technique to discretize the space derivatives. First, we partition the space domain $[a, b]$ uniformly as $a = x_0 < x_1, ..., x_M = b$ with space size $\Delta x = h = x_{m+1}-x_m = \frac{b-a}{M}\,$ where $m = 0, 1, \dots, M$. For discrete form, we denote $u_m^n = u(x_m, t_n)$ for $m = 0, 1, ..., M$ and $n = 0, 1, ..., N$. The CExpB-spline functions ${Ep}_m\left(x\right)$ for $m = -1, \ 0, \ 1, \dots, M+1$ are defined as ^[28,29,30]:

where

In Eq (2.2), the free parameter $\hat \kappa$ have been used for obtaining the different forms of CExpB-spline functions. The set of ${Ep}_m\left(x\right)\ \forall \ m = -1, \ 0, \ 1, \dots, M+1$ forms a basis over the problem domain. We assume that the approximation $u^M$ to the exact $u(x, t)$ at the point $(x_m, t_n)$ is expressed in terms of linear combinations of the CExpB-spline functions and unknown time-dependent quantities as follows:

where $C_i\left(t\right)$ are the unknown quantities which we have to evaluate for the approximated solution $u^M\left(x, t\right)$ at $\left(x_i, t_j\right)$. Since each CExpB-spline covers four elements, so each element is covered by four CExpB-splines. So the variation of the $u^M(x, t)$, over the element, can be written as:

Using Eq (2.4), the $u(x, t)$ and its first two derivatives at the knots in terms of $C_m^n$ are given as:

Now, we modify the CExpB-spline basis functions which generate a new set of CExpB-spline basis functions. The procedure, for modifying the basis functions, is given as ^[1,12,13]:

where $\{\hat{Ep}_0(x), \hat{Ep}_1(x)\dots \dots \hat{Ep}_M(x)\}$ forms a basis over the problem domain. Next, the modified form of the approximated solutions, as a linear combination of modified CExpB-spline functions, is given by

Next, the initial vector $\left(C^0_m, \ m = 0, \ 1, \dots, \ M\right)$ can be obtained by initial condition and boundary values of its derivatives as:

Eq (2.10) yields a $(M+1) \times (M+1)$ system as:

Now, discretizing the Eq (1.1) by using (2.1) and Crank-Nicolson method, we have

Using Eqs (2.4) and (2.6) in Eq (2.12) and simplifying the terms, we have

where

For $m = 0$, from Eq (2.13), we get

For m = M, from Eq (2.13), we get

At time step $t_{n}, \ n = 1, 2, ..., N$, the Eqs (2.13)–(2.15) can be formulated into a $(M+1)\times(M+1)$ linear system as follows:

where

3.   Stability analysis

The von-Neumann stability analysis ^[4,18,21,23] have been done in this section. Taking $f\left(x, t\right) = 0$ in Eq (1.1) and discretizing, we get

Simplifying and rearranging the terms, we get

which can be written as

where

and

Now we consider one Fourier mode out of the full solution $C_m^n = \delta ^n e^{ki \phi}$ as trial solutions at a given point $x_i$, where $\phi = \theta h$. The $\theta$ and $h$ are the mode number and element size respectively, and $k = \sqrt{-1}$. Substituting the trial solution in above equation and simplifying the terms, we get

where $\bar F = \sum\limits_{l = 1}^{n - 1} {\left({\bar b_{l-1} - \bar b_l } \right)} \delta ^{n - l} + \bar b_{n - 1} \delta ^0.$

Now we define

Using (3.4), we get

Now, using above equation and the values of $A, \ B, \ D$ and $E$ in Eq (3.3), we have

From (3.6), we get $|\delta| \le 1$ , and hence the method is unconditionally stable for the discretized system of the Caputo time-fractional diffusion equation.

4.   Computational results and discussions

In this section, we consider numerical examples of the Caputo time-fractional diffusion equation in order to check the accuracy and efficiency of the method. We use the following error norms :

The ROC is calculated by

where, $E^{h_1}$ and $E^{h_2}$ represent the errors at space mesh sizes $h_1$ and $h_2$, respectively and $E^{k_1}$ and $E^{k_2}$ represent the errors at time mesh sizes $k_1$ and $k_2$, respectively.

4.1.   Example 1

First, we consider the Caputo time-fractional diffusion equation

with the exact solution

where

We solve the Caputo time-fractional diffusion Eq (4.3) for different values of $M$ and $N$. Table 1 shows the comparison of the present method with the finite element method ^[6] and cubic B-spline collocation method ^[24] for different values of $M$, fixed $\Delta t = 0.001$ at the time-fractional order $\alpha = 0.5\ $ in terms of $L_2$ and $L_{max}$ error norms. This table shows that the present method gives better results than the results obtained by those available in ^[6,24]. Moreover, one can notice from Table 1 that the error norms are decreasing as we increase the space as well as time mesh sizes. The analytical and approximate solutions together with the absolute errors are presented in Figure 1 for $\hat \kappa = 11$, $\Delta x = 0.001$, and $\alpha = 0.5$ at $\Delta t = 0.2$, $\Delta t = 0.1$, and $\Delta t = 0.05$, respectively. Figure 2 shows the comparison between analytical and approximate solutions for $M = 64$, $\Delta t = 0.001$, and $\alpha = 0.5$ at different $t$. From these figures, we notice that there is an excellent agreement between analytical and approximate solutions. Also, the absolute errors are decreasing on increasing the time interval.

4.2.   Example 2

Next, we consider the Caputo time-fractional diffusion equation

with the exact solution

where

We solve the Caputo time-fractional diffusion Eq (4.4) for various values of $M$. Table 2 shows the maximum absolute error norms $L_{max}$ for different values of time-fractional orders i. e. $\alpha = 0.2, \ 0.5\ $ and $0.8$ at fixed time mesh size of $N = 50$ for different space mesh sizes. One can see from Tables 2 and 3 that the error norms are decreasing as we increase the mesh sizes. As we can notice from this table that the present method is of order $O\left((\Delta x)^2, (\Delta t)^{1+\alpha}\right)$. The approximate solutions for time levels $t =$ 0.1, 0.4, 0.6 and 0.8 are plotted in Figure 3 for $M = 32, \ N = 50$ and time-fractional order $\alpha = 0.5$. One can notice that the amplitude of the approximate solutions is increasing on incresing the time $t$.

4.3.   Example 3

Now, we consider the Caputo time-fractional diffusion equation

where

subject to the initial and boundary conditions

and with the exact solution

Now, we solve the Caputo time-fractional diffusion Eq (4.5) for different values of $M$ and $N$. Table 4 shows the comparison in term of $L_{max}\ $error norms obtained by the present method and the method based on the product-integration (PI) rule presented in ^[7] for $\Delta x = \frac{1}{101}$, $\beta = 6$, different values of $N$ and $\alpha$ at $t = 1$. As we can see that the present method gives more accurate results than the method presented in ^[7]. Also, we notice that the error norms are decreasing as we increase the time mesh sizes. Figure 4 shows the analytical and approximate solutions together with absolute error norms for $\Delta x = \frac{1}{101}, \ \alpha = 0.85, \ \beta = 5, \ \Delta t = 0.01$ at $T = 0.35$ showing that errors are decreasing on increasing the grid sizes.

5.   Conclusions

A modified CExpB-spline collocation technique has been presented for solving the Caputo time-fractional diffusion equation. The modified CExpB-spline collocation technique is used to discretize the space derivatives. The three examples of the Caputo time-fractional diffusion equation have been considered. The obtained results show that the present method gives more accurate results than the results obtained in ^[6,7,24]. It is observed numerically that the method is second-order accurate in space and $(1+\alpha)$ order in time. The stability analysis shows that the method is unconditionally stable. Moreover, the implementation of the present method is easy and needs low memory storage which is the advantage. The present method can easily be extended for solving higher dimensional fractional PDEs.

Conflict of interest

The authors have no conflict of interest.