Novel analysis of nonlinear dynamics of a fractional model for tuberculosis disease via the generalized Caputo fractional derivative operator (case study of Nigeria)

1.   Introduction

In mathematical biology, predicting the time evolution of a biomass or a population over a spatial domain is a very important problem. Often, Lotka-Voterra-type equations are used to describe population dynamics, whether studying competition or predator-prey models. Such systems have been studied in many settings. Conditions under which two species can coexist have been treated theoretically. Numerical models have been developed that hope to mimic behavior of the system represented in the models. In the case that dispersion is occurring to each species, biomass or population studies over a domain $\Omega_T = (0, T) \times \Omega$ for some $\Omega \subset \mathbb{R}^d$, generally with $d = 1, 2$, or 3 may take the general form

where $L$ is a linear operator on $\Omega$, distinct from the Laplacian, that can take various forms. If $a_1, a_2 > 0$ in system (1.1), then some degree of local diffusion of each species is modeled by the equations. The $b_1 = b_2 = 0$ case has been considered by multiple authors, such as in ^[1,2,3,4] and in their accompanying references.

In real world settings, populations in competition may demonstrate diffusion, dispersion, or some degree of both local and nonlocal behaviors. Competition for common resources may not only follow in their immediate neighborhood, but also in the entire spatial domain. In addition, this competition is not necessarily occurring only between individuals at the same location, but also between individuals at different locations; see ^[5] and references therein for an excellent motivation for and summary of models of nonlocal dispersion operators. Following similar reasoning, when $b_1 = b_2 = 1$, the operator $L$ given by

has been motivated as an accurate reflection of dispersion between species for suitable $J$, and will be used here. The authors in ^[6] show the derivation of the nonlinear dispersion operator as $J*u-(J*1)u$, but then point out that it is unrealistic to use a convolution term to model biological species in bounded domains. Further, the restriction $L$ in (1.2) of $J*u-(J*1)u$ arises naturally in cases of hostile surroundings, a periodic environment, or under reflected boundary conditions. In ^[7], these assumptions are shown to be unnecessary, where the dispersion operator $L$ in (1.2) is used to model nonlocal interaction. There, the authors show that for $L$, total internal energy is conserved and free energy decreases along trajectories so that $L$ a suitable choice to reflect dispersion as a replacement for local dispersion modeled by the Laplacian for a symmetric interaction kernel $J$. Thus, operator $L$ is also used both in nonlocal Allen-Cahn-type equations (c.f. ^[8]) or when describing population changes with nonlocal dispersion in biological systems (as in ^[9] and references therein).

Indeed, when restricted to $\Omega_T\subset (0, T)\times \mathbb{R}^d$, kernel function $J$ in (1.2) becomes a measure of the probability that population members at all positions x affect those at ${\bf y}\in R^d$, and vice-versa. Hence,

may be used to describe the rate at which members of a species are leaving ${\bf x}\in \Omega$ at time $t$ to travel to all other sites ${\bf y}\in \mathbb{R}^d$. Thus, for $b_1 = b_2 = 1$, various forms of system (1.1) using $L$ in (1.2) have also been studied in depth, including traveling wave solutions, spreading speed, stability of traveling waves, and of entire solutions as in ^[10,11,12], and references therein.

There has been a good deal of work on Lotka-Volterra models that include local diffusion, and more recently, nonlocal dispersion has been treated in such systems (see ^{[1,2,3,4,13,14]}, along with references therein). In ^[15], an implicit approach to the numerical analysis of the system is introduced that mimics the dynamical properties of the true solution. In addition, it is proven that the scheme introduced there is uniquely solvable and unconditionally stable. The asymptotic behavior of the difference scheme is studied by constructing upper and lower solutions for the difference scheme. The convergence rate of the numerical solution to the true solution of the system is also given.

Following notation in (1.1) and (1.2) is

where $\Omega$, $\Omega_T$, and $L$ are defined previously. Here, $u(t, {\bf x})\ge 0$ and $v(t, {\bf x})\ge 0$ denote the population densities of the competing species for time $t\ge 0$ and ${\bf x}\in\Omega\subset \mathbb{R}^d$ for $d = 1, 2$. In addition, we assume that the system whose solutions describe the density of each species that began with (1.1) has been nondimensionalized, so that in (1.3), $a, b, K_1, K_2 > 0$ depend on the initial choice of constants in (1.1). To our knowledge, this system has not been analyzed with an unconditionally stable, nonstandard numerical method whose convergence rate can be given. This is the goal of this contribution.

As described previously, the system is used to model the two species competing with each other for the same prey, where both species are continuously distributed in time $t$ throughout a region $\Omega$, with each exhibiting free movement in the form of nonlocal dispersion. The method introduced to discretize (1.3) and the study of its properties will follow a similar development to the one used in ^[16], where properties of a single integro-differential equation that models an Ising spin system, with a convolution term that involves the Kac potential, are discussed in detail (see ^{[17,18,19,20]}).

In Section 2, we introduce the difference scheme used for the approximation of (1.3) over $\Omega_T\subset(0, T)\times\mathbb{R}$. We prove existence of the numerical solution to the scheme and that this solution is stable, independent of the choice of $\Delta t$ and $\Delta {\bf x}$. We give the convergence rate of the numerical scheme to the true solution. In Section 3, we present some results of numerical experiments that confirm the stability and convergence of the proposed difference scheme in one- and two-dimensional spatial domains $\Omega$. In Section 4, we provide a summary of the results.

2.   A nonstandard numerical scheme

Analysis will be carried out over a domain $\Omega$ in one-dimensional space. All results carry over naturally to higher-dimensional space. For $t > 0$ we introduce time step $t_k = k\Delta t$ for $k = 0, 1, 2, \ldots,$ where $\Delta t$ is of fixed size to be determined later. On the interval $\Omega = (-L, L)\subset \mathbb{R}$ we define the partition

Using $u^k_i$ and $v_i^k$ to represent the numerical approximation to the true solutions $u$ and $v$ to (1.3) at $(t_k, x_i)$, our choice of difference scheme for $n = 1$ in (1.3) is nonstandard to invoke desirable properties that will be established later, namely

for $k = 0, 1, 2, \ldots$ and for $0\leq i\leq N$. Throughout, for convenience, the discretization of $Lu$ as given in (1.2) for (2.1) will be denoted by

a similar expression is used for $Lv$. We also introduce the initial conditions in (1.3) as

for $i = 0, 1, 2, \ldots, N$, where for all $i$, $u_0(x_i), v_0(x_i)\ge 0$.

Solving (2.1) for $u_i^{k+1}$ and $v_i^{k+1}$ gives the iteration scheme for $n = 1$ as

for $k = 0, 1, 2, \ldots$ and $i = 0, 1, 2, \ldots, N$.

Although we present and prove theorems for $n = 1$, it is useful for programming numerical models to show the method of approximating (1.3) for $n = 2$ as well. In this case, we choose $\Omega = (-L, L)\times(-W, W) \subset \mathbb{R}^2$, and partitions

and

where $\triangle x = 2L/M$ and $\triangle y = 2W/N.$

The difference scheme for (1.3) includes

together with

all for $i = 0, 1, \ldots, M$ and $j = 0, 1, \ldots, N$, where in (2.4),

From (2.3) and (2.4), we arrive at the explicit finite difference scheme

for $i = 0, 1, \ldots, M$, $j = 0, 1, \ldots, N$, and $k = 0, 1, 2, \ldots.$

Theorem 2.1. For $n = 1$ and for the initial conditions in (1.3), let $m_1 = \max u_0(x)$, $m_2 = \max v_0(x)$, $M_1 = \max\{K_1, m_1\}$, and $M_2 = \max\{K_2, m_2\}$. Then for all $k = 0, 1, 2, \ldots$ and for $i = 0, 1, 2, \ldots, N$,

Hence, the numerical scheme (2.2) is unconditionally nonnegative and unconditionally stable.

Proof. We proceed by induction. For $k = 0$,

for $i = 0, 1, 2, \ldots, N$, so that (2.6) holds.

Assume now that (2.6) holds for some $k\in\mathbb{N}$. Then for $k+1$,

Similarly, $v_i^{k+1}\le M_2,$ so the result holds for $k+1$ if it holds for $k$. Therefore, by mathematical induction, (2.6) is valid for all $k = 0, 1, 2\ldots$ and for $i = 0, 1, 2, \ldots, N$. □

We now turn to the question of convergence of the difference equations (2.2) to the true solution of (1.3).

Theorem 2.2. If $u, v\in C^{1, 2}([0, T]\times \overline{\Omega})$ are solutions to (1.3), then the solution of (2.2) converges to $u$ and $v$ as $\Delta t, \Delta x\rightarrow 0$, uniformly on $[0, T]$, with rate $\mathcal{O}(\Delta t+\Delta x^2)$.

Proof. Let $(u(t, x), v(t, x))$ represent the solution pair to (1.3), where $u, v\in C^{1, 2}([0, T]\times \overline{\Omega})$. Set

We will prove the convergence claim based on $u$, then the same will follow for $v$ by the symmetry of equations in (1.3). Let $\Delta t = T/K$, so that $t_k = k\Delta t$ for $k = 0, 1, 2, \ldots, K$. From (1.3) and (2.1) we have

where $R_u$ is a function with $R_u(\Delta t, \Delta x) = \mathcal{O}(\Delta t+\Delta x^2)$. Let

for $k = 0, 1, 2, \ldots, K$ and $i = 0, 1, 2, \ldots, N.$ Then $X^0_i = 0, \ Y_i^0 = 0$, for $i = 0, 1, 2, \ldots, N.$ Using (2.1) in conjunction with (2.7),

so that from (2.8),

We turn to upper bounds on each of the terms in (2.9). To accomplish this, for each $k$, $k = 0, 1, 2, \ldots, K$, we will use $W_u^k = \max_{i} \left| U_i^k-u_i^k\right|$ and $W_v^k = \max_{i} \left|V_i^k-v_i^k\right|$. Setting $C_1 = \max_{i}(J*1)_i$,

and

Now, since $\left|u_i^k\right|$ and $\left|v_i^k\right|$ are uniformly bounded by $M_1$ and $M_2$, from (2.1),

where $C(M_1, M_2, C_1)$ is a constant that depends only on $M_1$, $M_2$, and $C_1$. Hence there exists $C_2$ such that

and

where $K_1$ is that constant related to carrying capacity in (1.3), so that for some constants $C_3$, $C_4$, and $C_5$,

and thus there exist constants $C_6$, $C_7$, and $C_8$ with

Substituting (2.10) and (2.11) into (2.9), we obtain

for $i = 0, 1, \ldots, N$, where $C_8$ and $C_9$ are constants independent of $i$ and $k$, and where $\triangle t^2$-terms are absorbed into $\Delta t\, R_u(\Delta t, \Delta x)$. Therefore, for each $k$, $k = 0, 1, 2, \ldots, K$,

Similarly, there exist $C_{10}$ and $C_{11}$ with

where, as with $R_u$, $R_v$ is a function with $R_v(\Delta t, \Delta x) = \mathcal{O}(\Delta t+\Delta x^2)$. Now setting $Z^k = W_u^k+W_v^{k}$, from (2.13) and (2.14), there exists a constant $C_0$ with

for $k = 0, 1, 2, \ldots, K$. Set

Then since $Z^0 = 0$, using (2.15) and (2.16) and iterating,

for $k = 0, 1, 2, \ldots, K-1$. Since $e^x\geq 1+x$, it follows that $e^{K x}\geq (1+x)^K$, so that for all $k$, $k = 0, 1, 2, \ldots, K-1$, and again using $D$ from (2.16),

Thus, for $k = 0, 1, 2, \ldots, K-1$,

so that $Z^{k}\rightarrow 0$ for $k = 0, 1, 2, \ldots, K$ as $\Delta t\rightarrow 0$, $\Delta x\rightarrow 0$. This completes the proof. □

Remark. Similar results hold for $n = 2$ in Theorems 2.1 and 2.2, and their proofs.

3.   Numerical results

In this section we finish by presenting some results of computational experiments that verify the stability and convergence of the proposed difference scheme, confirming that the numerical solutions preserve the properties of the theoretical solution as well as those guaranteed by Theorem 2.2. Since there is no exact solution to compare with the approximation generated by the difference scheme, we use fix $\triangle x$ and compute for various $\triangle t$ values, then vice-versa. We compare the results in tables. We also present graphical results for dimensions $n = 1$ in (2.2) and $n = 2$ in (2.5).

${\rm{C}}{\small{\rm{ASE}}}\; {\rm{I}}$. For $n = 1$, we test method (2.1) for $\Omega = (-1, 1)$, $\epsilon = 0.1$, $a = 0.4$, $b = 0.6$, $K_1 = K_2 = 1$, $u_0(x) = 0.2\cos(2\pi x)+1$, and $v_0(x) = 0.3\sin(2\pi x)+1$, where $J(x) = (\epsilon \sqrt{\pi})^{-1}\exp\left(-x^2/\epsilon^2\right)$. We call these approximations $(u(t, x), v(t, x))$. Their convergence to steady state solutions is demonstrated in Figures 1 and 2 for $\Delta t = \Delta x = 0.05$, as an example, since convergence is independent of time and space steps and graphs look much the same for any reasonable choices of small $\Delta t$ and $\Delta x$.

${\rm{C}}{\small{\rm{ASE}}}\; {\rm{I}}{\small{\rm{A}}}$: Hold $\triangle x = 0.05$. Let $(u(t, x), v(t, x))$ denote the numerical solution under the parameters as chosen above corresponding to $\triangle t$, while $(u_1(t, x), v_1(t, x))$ corresponds to $\triangle t_1$. Table 1 shows the maximum absolute errors, $\max|u(t, x)-u_1(t, x)|$ and $\max|v(t, x)-v_1(t, x)|$, at $t = 5$ across $\Omega$.

We note that the reduction of $\Delta t$ by a factor of $0.5$ reduces the error by $\mathcal{O}(\Delta t)$, as predicted by Theorem 2.2.

${\rm{C}}{\small{\rm{ASE}}}\; {\rm{I}}{\small{\rm{B}}}$: We fix $\triangle t = 0.1$ and vary $\triangle x$. As before, we let $(u(t, x), v(t, x))$ represent the numerical solutions corresponding to $\triangle x$, and let $(u_1(t, x), v_1(t, x))$ represent the numerical solutions corresponding to $\triangle x_1$. We compare the error differences $\max|u(t, x)-u_1(t, x)|$ and $\max|v(t, x)-v_1(t, x)|$ at $t = 5$ across $\Omega$.

We note that the reduction of $\Delta x$ by a factor of $0.5$ reduces the error by a factor of $\mathcal{O}(\Delta x^2)$, or about 0.25, as stated in Theorem 2.2.

${\rm{C}}{\small{\rm{ASE}}}\; {\rm{II}}$. For $n = 2$, let $\Omega = (-1, 1)\times (-1, 1)$, and let $\epsilon = 0.1$, $a = 0.4$, $b = 0.6$, $u_0(x, y) = 0.4+0.2\cos(2\pi x)\cos(2\pi y)$, and $v_0(x, y) = 0.5+0.3\sin(2\pi x)\sin(2\pi y)$, where $J_{\epsilon}(x, y) = \frac{1}{\epsilon^2\pi}\exp(-\frac{x^2+y^2}{\epsilon^2})$.

We first show graphs of some approximate solutions generated by the two-dimensional method (2.5) for $\Delta t = 0.25$, and $\Delta x = \Delta y = 0.2$ in Figures 3 and 4. As in the one-dimensional case, since convergence is independent of time and space steps and graphs look much the same for any reasonable choices of small $\Delta t$ and $\Delta x$, we have chosen these values as a representative of any such reasonable choice.

${\rm{C}}{\small{\rm{ASE}}} \; {\rm{II}}{\small{\rm{A}}}$: We hold $\triangle x = \triangle y = 0.1$ and compare accuracy for various $\triangle t$-values in (2.5). Denote the numerical solution $(u(t, x, y), v(t, x, y))$ as the one generated by (2.5) corresponding to $\triangle t$ and $(u_1(t, x, y), v_1(t, x, y))$ corresponding to $\triangle t_1$. We compare the differences $\max|u(t, x, y)-u_1(t, x, y)|$ and $\max|v(t, x, y)-v_1(t, x, y)|$ at $t = 5$ across $\Omega$ in Table 3.

${\rm{C}}{\small{\rm{ASE}}} \; {\rm{II}}{\small{\rm{B}}}$: Finally, we carry out the same accuracy test for fixed $\Delta t = 0.25$ and various $\Delta x = \Delta y$ and $\Delta x_1 = \Delta y_1$ values at time $t = 5$. The results are displayed in Table 4.

All approximations in the tables show convergence at the rates predicted by Theorem 2.2.

4.   Conclusions

The foregoing results have motivated the use of a Lotka-Volterra-type equation with operator $L$ that reflects intra-species dispersion, or nonlocal interaction, with competition between species whose populations are given by $u$ and $v$. A nonstandard numerical scheme was introduced that is stable, independent of the choice of time step, and that yields biologically sensible (nonnegative) numerical approximations to populations $u$ and $v$ of this system. Moreover, this nonstandard scheme was shown to be convergent to the solution of the proposed system and the order of convergence given. Because its convergence was established, it is possible to state with confidence that accurate solutions to the system are shown in the numerical experiments that were offered to confirm the results.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank the reviewers for their help and insightful suggestions for improvement of this paper.

Conflict of interest

The authors declare there are no conflicts of interest.