A Collocation Technique via Pell-Lucas Polynomials to Solve Fractional Differential Equation Model for HIV/AIDS with Treatment Compartment

1  Introduction

Fractional differential equations are quite popular among scientists in order to model various stable physical phenomena. Recently, fractional order derivatives have attracted great attention due to their numerous applications in nonlinear complex systems arising in fluid mechanics, damping laws, electrical networks, signal processing, diffusion-reaction processes, relaxation processes, electrochemistry, mathematical biology, physics and various important phenomena in other branches of science and engineering. Fractional derivatives provide more precise models of real-life problems than integer order derivatives. Because some systems exhibit memory, history, or non-local effects that may be difficult to model with the help of integer order derivatives, many natural phenomenon problems are modeled utilizing fractional calculus. In recent years, different types of powerful techniques such as the generalized differential transformation method [1], the Adomian decomposition method [2], the homotopy analysis method [3], the homotopy analysis transformation method [4], the modified Laplace transform method [5] and the homotopy perturbation transform method [6,7] have been introduced to find the approximate solution of the fractional model of this type of differential equations.

Fractional differential equations are frequently used in mathematical modeling. Mathematical modeling can provide an understanding of the mechanisms underlying the transmission and spread of disease, help identify key factors in the disease transmission process, recommend effective control and preventive measures, and provide an estimate of the severity and potential scale of the epidemic. Mathematical biology is a multidisciplinary field of study. Mathematical biology is also used as an important tool in medicine, as well as in differential equations, basic sciences, and many fields of engineering such as genetic engineering, biomedical engineering, and clinical engineering. Lately, mathematical biology has been successfully applied to several important fields in medicine including epidemiology, genetics, drug design and discovery, biofluids, cardiovascular diseases, immunology, microbiology, neuroscience, oncology, virology, and more. Thus, medical phenomena are better understood, and practical action paths are found. The improvements in healthcare and quality of life are achieved because there are important contributions such as early diagnosis, effective medicines, and control of epidemics. From past to present, many models such as the logistic equation and the Lotka-Volterra equations based on prey and predator populations [8], HIV (human immunodeficiency virus) infection model [9–12], SIR model (Susceptible-Infected-Removed) and COVID-19 model [13], etc., have been developed in the field of medicine using mathematical biology.

Recently, the studies on mathematical modeling of human immunodeficiency virus (HIV) have become very popular. HIV, which causes the acquired immunodeficiency syndrome (AIDS), devastates the human body’s skill to fight infections. This disease is very hazardous and can be mortal if left untreated. The first AIDS case was detected in 1981 [14]. In 2017, the US Center for Disease Control and Prevention (CDC) [15] declared that if HIV/AIDS is not treated with antiretroviral drugs, HIV infection progresses in several phases. According to some official reports, it is observed that between 2000 and 2016, HIV-related deaths in Africa decreased by one-third thanks to ART (Antiretroviral Therapy) [16]. There are many methods such as the exponential collocation method (ECM) [17], the Bessel collocation method (BCM) [18], the differential transform method (DTM) [19], the optimization method [20], the Galerkin-like method [21], Laplace Adomian decomposition method (LADM) [22], homotopy perturbation method (HPM) [23], variational iteration method (VIM) [24] for solving the HIV infection model. In addition, there are many methods such as homotopy perturbation method, variational iteration method and the Adomian decomposition method (ADM) [25], the fourth kind Chebyshev wavelet method [26], shifted Legendre collocation method [27], the homotopy analysis method [28], LADM [29,30], the Legendre wavelet approach [31], the septic B-spline scheme [32], the fractional approximation method [33] for solving fractional order HIV models. On the other hand, there are some studies conducted for the HIV/AIDS epidemic model with antiretroviral therapy (ART). Luo et al. [34] gave the global stability of disease-free equilibrium and the endemic equilibrium, investigated the long-time stochastic dynamic of the model, and gave some numerical simulations. Wang et al. [35] presented the optimal control problem and conducted numerical simulations. Chen et al. [36] developed a type-2 fuzzy logic controller for antiretroviral therapy of HIV infection and performed simulations for two strategies. Ali et al. [37] constructed the existence and uniqueness conditions of the HIV/AIDS model by utilizing Schaefer and Banachtype fixed point theorems. The model’s qualitative analysis investigated, determined the existence/uniqueness of the solution to the HIV/AIDS model [38], established the stability results for the system by incorporating the Ulam-Hyers method, and applied Newton’s polynomial and the Toufik-Atangana numerical method.

Spectral methods are one of the common methods used for solving mathematical models because they are extremely sensitive. Inasmuch as it uses linear combinations of the orthogonal polynomials as basis functions, these methods lead to accurate approximate solutions [39,40]. Basically, three types of these methods can be identified. These are Tau [41], collocation [42,43] and Galerkin [44] and these are based on orthogonal systems such as Bernstein polynomials, Bessel polynomials, Chebyshev polynomials, Hermite polynomials, Legendre polynomials, Laguerre polynomials, Pell-Lucas polynomials, etc. Compared to other approximation methods, the presented method requires fewer calculations. With the presented method, an approximate solution can be achieved with small N values selected for problems that do not have an exact solution. However, in some methods in the literature, more iterations are required to obtain approximate solutions for systems of the nonlinear equation or even nonlinear equations [45].

On the other hand, there are many methods based on the Pell-Lucas polynomials (PLPs) on the solutions of Fredholm-type delay integro-differential equations [46], functional differential equations [47], population models [8], Fredholm-Volterra integro-differential equations [48], parabolic-type partial integro-differential equations [49], nonlinear Lane-Emden pantograph differential equations [50] and an SIR model on the spread of the novel coronavirus (nCoV-2019) pandemic [13] and very effective results have been obtained from these methods. However, there is no method based on the Pell-Lucas polynomial solutions of the fractional order HIV/AIDS epidemic model with a treatment compartment (FHEMTC) in the literature. This reveals the importance and novelty of this study. For this reason, this study is also important as it is a new application of the method by developing it for fractional problems. In this paper, we consider the following FHEMTC [14,15]:

 CDt(p)S(t)=Λ−βS(t)I(t)−μ1S(t)−dS(t) CDt(p)I(t)=βS(t)I(t)+α1T(t)−dI(t)−k1I(t)−k2I(t) CDt(p)A(t)=k1I(t)−(δ1+d)A(t)+α2T(t) CDt(p)T(t)=k2I(t)−α1T(t)−(d+δ2+α2)T(t) CDt(p)R(t)=μ1S(t)−dR(t)S(0)=S0,I(0)=I0,A(0)=A0,T(0)=T0,R(0)=R0.(1)

Representations of parameters and variables in (1) are given in Table 1.

In this study, we examine the approximate solutions based on PLPs of (1):

SN(t)=∑j=0Nc1,jQj(t),IN(t)=∑j=0Nc2,jQj(t),AN(t)=∑j=0Nc3,jQj(t),TN(t)=∑j=0Nc4,jQj(t),RN(t)=∑j=0Nc5,jQj(t).(2)

Here, ci,j (i=1,2,3,4,5) are the unknown coefficients, N>0 and Qj(t) are PLPs and these polynomials are defined by [51,52]

Qj(t)=∑p=0⟦i/2⟧2j−2pjj−p(j−pp)tj−2p.(3)

Please see [51,52] to learn more about PLPs.

The flow of this paper is created as follows: In Section 2, the required matrix relations of the method are presented. In Section 3, the numerical method for FHEMTC is presented. In Section 4, the error analysis is given. In Section 5, the parameters and initial conditions in the HIV/AIDS epidemic model are determined and Pell-Lucas collocation method (PLCM) is applied to the obtained model. In Section 6, a brief conclusion of the article is presented.

2  Preliminaries and Matrix Relations

In this section, we define the fundamental facts on fractional derivatives and then we express the terms of FHEMTC (1) in matrix form.

Lemma 2.1. PLPs (3) are written the following matrix form [50]:

QN(t)=FN(t)KN,(4)

where, QN(t)=[Q0(t)Q1(t)⋯QN(t)], FN(t)=[1tt2⋯tN]. When N is even and odd, we can, respectively, write the matrix KNT as follows:

KNT=[200⋯002111(10)0⋯02021(11)02222(20)⋯0⋮⋮⋮⋱⋮20NN2(N2N2)022NN+22(N+22N−22)⋯2NNN(N0)],

and

KNT=[200⋯002111(10)0⋯02021(11)02222(20)⋯0⋮⋮⋮⋱⋮021NN+12(N+12N−12)0⋯2NNN(N0)].

Proof. If the vector FN(t) is multiplied by KN from the right side, then (4) is obtained. ■

Lemma 2.2. The assumed solution forms (2) can also be written in matrix form as

SN(t)=FN(t)KNC1,N,IN(t)=FN(t)KNC2,N,AN(t)=FN(t)KNC3,N,TN(t)=FN(t)KNC4,N,RN(t)=FN(t)KNC5,N,(5)

where, Ck,N=[ck,0ck,1⋯ck,N]T(k=1,2,3,4,5) and other matrices are as in Lemma 2.1.

Proof. If QN(t)=FN(t)KN is multiplied, separately, by C1,N, C2,N, C3,N, C4,N and C5,N from the right, we have (5). ■

Definition 2.1. Using [53], we describe the fractional derivative of h(t) in the Caputo sense as

 CDt(p)h(t)=Jm−pDmh(t)=1Γ(m−p)∫0t(t−x)m−p−1h(m)(x)dx,

for m−1<p<m, m∈N, t>0, h∈C−1m where D=ddx. When h(t) equals tη, the Caputo fractional derivative of h(t) can be written as [54]

 CDt(p)tη={0,forη∈N0andη<⌈p⌉Γ(η+1)Γ(η+1−p)tη−p,forη∈N0andη≥⌈p⌉orη∉N0andη>⌊p⌋(6)

Here, if c is constant, then  CDt(p)c=0.

Lemma 2.3. The matrix forms of the p-th order fractional derivatives of the assumed solution forms (2) in (1) become

SN(p)(t)=FN(p)(t)KNC1,N,IN(p)(t)=FN(p)(t)KNC2,N,AN(p)(t)=FN(p)(t)KNC3,N,TN(p)(t)=FN(p)(t)KNC4,N,RN(p)(t)=FN(p)(t)KNC5,N.(7)

Here,

FN(p)(t)=[0Γ(2)Γ(2−p)t1−pΓ(3)Γ(3−p)t2−p⋯Γ(N+1)Γ(N+1−p)tN−p](8)

and other matrices are given in Lemma 2.2.

Proof. The p-th order fractional derivative of FN(t) can be written as in (8) with the help of (6). Then, the p-th order fractional derivative of the assumed solution form (2) is taken. The proof is completed by using (8) in the derived terms. ■

Lemma 2.4. The nonlinear term in (1) can be expressed the following matrix form:

SN(t)IN(t)=(FN(t)KNC1,N)(FN(t)KNC2,N).(9)

Here, all matrices are as in Lemma 2.2.

Proof. Through Lemma 2.2, we write SN(t)=FN(t)KNC1,N and IN(t)=FN(t)KNC2,N. Consequently, the matrix multiplication of these two terms gives us the proof of Lemma. ■

Lemma 2.5. The initial conditions in (1) can be written with matrix relations

LNC1,N=S0,LNC2,N=I0,LNC3,N=A0,LNC4,N=T0,LNC5,N=R0.(10)

Here, LN=FN(0)KN. Also, other matrices are given in Lemma 2.2.

Proof. By substituting t→0 in (5), the proof is completed. Let’s note that we denote the matrix product FN(0)KN by LN. ■

3  Numerical Method

In this part, PLCM is given for the approximate solution of FHEMTC (1) by utilizing collocation points.

Theorem 3.1. Supposed that the approximate solutions of the HIV/AIDS epidemic model (1) are in form (2), then, the following matrix relations are obtained:

FN(p)(t)KNC1,N=Λ−β(FN(t)KNC1,NFN(t)KNC2,N)−μ1FN(t)KNC1,N−dFN(t)KNC1,NFN(p)(t)KNC2,N=β(FN(t)KNC1,NFN(t)KNC2,N)+α1FN(t)KNC4,N−dFN(t)KNC2,N−k1FN(t)KNC2,N −k2FN(t)KNC2,NFN(p)(t)KNC3,N=k1FN(t)KNC2,N−(δ1+d)FN(t)KNC3,N+α2FN(t)KNC4,NFN(p)(t)KNC4,N=k2FN(t)KNC2,N−α1FN(t)KNC4,N−(d+δ2+α2)FN(t)KNC4,NFN(p)(t)KNC5,N=μ1FN(t)KNC1,N−dFN(t)KNC5,N.(11)

Here, all matrices are given in Lemmas 2.2–2.4.

Proof. Firstly, Eq. (7) are written instead of fractional-order derivative terms in model (1), Eq. (9) is written instead of nonlinear term in model (1) and Eq. (5) are written instead of solution forms in the model (1). Hence, the system (11) is obtained, so the proof is completed. ■

Definition 3.1. The collocation points in interval [0,h] are defined by

tk=hNk,k=0,1,…,N.(12)

Theorem 3.2. Assume that the approximate solutions of (1) are given in (5). In that case, FHEMTC (1) becomes as follows:

G0C1,N+D1,0C2,N+(μ1+d)M0C1,N=Λ,G0C2,N+D2,0C2,N−α1M0C4,N=0,G0C3,N−k1M0C2,N+(δ1+d)M0C3,N−α2M0C4,N=0,G0C4,N−k2M0C2,N+(α1+d+δ2+α2)M0C4,N=0,G0C5,N−μ1M0C1,N+dM0C5,N=0,G1C1,N+D1,1C2,N+(μ1+d)M1C1,N=Λ,G1C2,N+D2,1C2,N−α1M1C4,N=0,G1C3,N−k1M1C2,N+(δ1+d)M1C3,N−α2M1C4,N=0,G1C4,N−k2M1C2,N+(α1+d+δ2+α2)M1C4,N=0,G1C5,N−μ1M1C1,N+dM1C5,N=0,⋮GNC1,N+D1,NC2,N+(μ1+d)MNC1,N=Λ,GNC2,N+D2,NC2,N−α1MNC4,N=0,GNC3,N−k1MNC2,N+(δ1+d)MNC3,N−α2MNC4,N=0,GNC4,N−k2MNC2,N+(α1+d+δ2+α2)MNC4,N=0,GNC5,N−μ1MNC1,N+dMNC5,N=0,LNC1,N=S0,LNC2,N=I0,LNC3,N=A0,LNC4,N=T0,LNC5,N=R0(13)

where,

Gi=FN(p)(ti)KN,D1,i=βFN(ti)KNC1,NFN(ti)KN,D2,i=−βFN(ti)KNC1,NFN(ti)KN+(d+k1+k2)FN(ti)KN,Mi=FN(ti)KN.(14)

Also, other matrices are given in Theorem 3.1. and Lemma 2.5.

Proof. The collocation points (12) are substituted into (11) and thus the following system of algebraic equations is obtained:

FN(p)(t0)KNC1,N=Λ−β(FN(t0)KNC1,NFN(t0)KNC2,N)−μ1FN(t0)KNC1,N−dFN(t0)KNC1,NFN(p)(t0)KNC2,N=β(FN(t0)KNC1,NFN(t0)KNC2,N)+α1FN(t0)KNC4,N−dFN(t0)KNC2,N−k1FN(t0)KNC2,N−k2FN(t0)KNC2,NFN(p)(t0)KNC3,N=k1FN(t0)KNC2,N−(δ1+d)FN(t0)KNC3,N+α2FN(t0)KNC4,NFN(p)(t0)KNC4,N=k2FN(t0)KNC2,N−α1FN(t0)KNC4,N−(d+δ2+α2)FN(t0)KNC4,NFN(p)(t0)KNC5,N=μ1FN(t0)KNC1,N−dFN(t0)KNC5,NFN(p)(t1)KNC1,N=Λ−β(FN(t1)KNC1,NFN(t1)KNC2,N)−μ1FN(t1)KNC1,N−dFN(t1)KNC1,NFN(p)(t1)KNC2,N=β(FN(t1)KNC1,NFN(t1)KNC2,N)+α1FN(t1)KNC4,N−dFN(t1)KNC2,N−k1FN(t1)KNC2,N−k2FN(t1)KNC2,NFN(p)(t1)KNC3,N=k1FN(t1)KNC2,N−(δ1+d)FN(t1)KNC3,N+α2FN(t1)KNC4,NFN(p)(t1)KNC4,N=k2FN(t1)KNC2,N−α1FN(t1)KNC4,N−(d+δ2+α2)FN(t1)KNC4,NFN(p)(t1)KNC5,N=μ1FN(t1)KNC1,N−dFN(t1)KNC5,N⋮FN(p)(tN)KNC1,N=Λ−β(FN(tN)KNC1,NFN(tN)KNC2,N)−μ1FN(tN)KNC1,N−dFN(tN)KNC1,NFN(p)(tN)KNC2,N=β(FN(tN)KNC1,NFN(tN)KNC2,N)+α1FN(tN)KNC4,N−dFN(tN)KNC2,N−k1FN(tN)KNC2,N−k2FN(tN)KNC2,NFN(p)(tN)KNC3,N=k1FN(tN)KNC2,N−(δ1+d)FN(tN)KNC3,N+α2FN(tN)KNC4,NFN(p)(tN)KNC4,N=k2FN(tN)KNC2,N−α1FN(tN)KNC4,N−(d+δ2+α2)FN(tN)KNC4,NFN(p)(tN)KNC5,N=μ1FN(tN)KNC1,N−dFN(tN)KNC5,N.(15)

Now, the relations in (14) are used in the system (15). Then, by writing this system and conditions (10) as a single system, we obtain (13), which completes proof. ■

Corollary 3.1. When nonlinear algebraic system (13) is solved through MATLAB, the unknown coefficients Ci,N (i=1,2,3,4,5) in (5) are calculated. Thus, the approximate solutions of FHEMTC (1) are obtained.

4  Error Analysis

This section has two purposes. The first purpose is to get an error bound for the presented method. The second purpose is to create an error estimation method. Let S(t), I(t), A(t), T(t), R(t) be the exact solutions and SN(t), IN(t), AN(t), TN(t), RN(t) be the approximate solutions based on PLPs of FHEMTC (1).

Theorem 4.1. (Upper Boundary of Errors) Assume that the generalized Maclaurin series are represented by SNMac(t)=FN(t)C~1,N, INMac(t)=FN(t)C~2,N, ANMac(t)=FN(t)C~3,N, TNMac(t)=FN(t)C~4,N, RNMac(t)=FN(t)C~5,N. Then, the errors of the approximate solutions based on PLPs of (1) are bounded by

‖S(t)−SN(t)‖∞≤λN(‖C~1,N‖∞+‖KN‖∞‖C1,N‖∞)+hN+1(N+1)!‖S(N+1)(ct)‖∞,‖I(t)−IN(t)‖∞≤λN(‖C~2,N‖∞+‖KN‖∞‖C2,N‖∞)+hN+1(N+1)!‖I(N+1)(ct)‖∞,‖A(t)−AN(t)‖∞≤λN(‖C~3,N‖∞+‖KN‖∞‖C3,N‖∞)+hN+1(N+1)!‖A(N+1)(ct)‖∞,‖T(t)−TN(t)‖∞≤λN(‖C~4,N‖∞+‖KN‖∞‖C4,N‖∞)+hN+1(N+1)!‖T(N+1)(ct)‖∞,‖R(t)−RN(t)‖∞≤λN(‖C~5,N‖∞+‖KN‖∞‖C5,N‖∞)+hN+1(N+1)!‖R(N+1)(ct)‖∞.(16)

Here, ‖FN(t)‖∞≤max{hN,1}:=λN, △Ci,N=‖Ci,N+1‖∞−‖Ci,N‖∞ (i = 1, 2, 3, 4, 5). Also, the coefficient matrices of SNMac(t), INMac(t), ANMac(t), TNMac(t), RNMac(t) are indicated by C~i,N (i = 1, 2, 3, 4, 5), respectively.

Proof. Firstly, using the triangle inequality, we can write

‖S(t)−SN(t)‖∞≤‖S(t)−SNMac(t)‖∞+‖SNMac(t)−SN(t)‖∞,‖I(t)−IN(t)‖∞≤‖I(t)−INMac(t)‖∞+‖INMac(t)−IN(t)‖∞,‖A(t)−AN(t)‖∞≤‖A(t)−ANMac(t)‖∞+‖ANMac(t)−AN(t)‖∞,‖T(t)−TN(t)‖∞≤‖T(t)−TNMac(t)‖∞+‖TNMac(t)−TN(t)‖∞,‖R(t)−RN(t)‖∞≤‖R(t)−RNMac(t)‖∞+‖RNMac(t)−RN(t)‖∞.(17)

Secondly, using Lemma 2.2, we know the matrix representations of the Pell-Lucas polynomial solutions from Lemma 2.2. Therefore, we have

‖SNMac(t)−SN(t)‖∞=‖FN(t)(C~1,N−KNC1,N)‖∞≤‖FN(t)‖∞(‖C~1,N‖∞+‖KN‖∞‖C1,N‖∞),‖INMac(t)−IN(t)‖∞=‖FN(t)(C~2,N−KNC2,N)‖∞≤‖FN(t)‖∞(‖C~2,N‖∞+‖KN‖∞‖C2,N‖∞),‖ANMac(t)−AN(t)‖∞=‖FN(t)(C~3,N−KNC3,N)‖∞≤‖FN(t)‖∞(‖C~3,N‖∞+‖KN‖∞‖C3,N‖∞),‖TNMac(t)−TN(t)‖∞=‖FN(t)(C~4,N−KNC4,N)‖∞≤‖FN(t)‖∞(‖C~4,N‖∞+‖KN‖∞‖C4,N‖∞),‖RNMac(t)−RN(t)‖∞=‖FN(t)(C~5,N−KNC5,N)‖∞≤‖FN(t)‖∞(‖C~5,N‖∞+‖KN‖∞‖C5,N‖∞).(18)

Owing to 0≤t≤h, the term ‖FN(t)‖∞ can be written as follows:

‖FN(t)‖∞≤max{hN,1}:=λN.(19)

Substituting (19) into (18), the following inequalities are obtained:

‖SNMac(t)−SN(t)‖∞≤λN(‖C~1,N‖∞+‖KN‖∞‖C1,N‖∞),‖INMac(t)−IN(t)‖∞≤λN(‖C~2,N‖∞+‖KN‖∞‖C2,N‖∞),‖ANMac(t)−AN(t)‖∞≤λN(‖C~3,N‖∞+‖KN‖∞‖C3,N‖∞),‖TNMac(t)−TN(t)‖∞≤λN(‖C~4,N‖∞+‖KN‖∞‖C4,N‖∞),‖RNMac(t)−RN(t)‖∞≤λN(‖C~5,N‖∞+‖KN‖∞‖C5,N‖∞).(20)