# A Higher-Order Galerkin Time Discretization and Numerical Comparisons for Two Models of HIV Infection

Authors: Attaullah; Şuayip Yüzbaşı; Sultan Alyobi; Mansour F. Yassen; Wajaree Weera
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

Human immunodeficiency virus (HIV) infection affects the immune system, particularly white blood cells known as CD4+ T-cells. HIV destroys CD4+ T-cells and significantly reduces a human’s resistance to viral infectious diseases as well as severe bacterial infections, which can lead to certain illnesses. The HIV framework is defined as a system of nonlinear first-order ordinary differential equations, and the innovative Galerkin technique is used to approximate the solutions of the model. To validate the findings, solve the model employing the Runge-Kutta (RK) technique of order four. The findings of the suggested techniques are compared with the results obtained from conventional schemes such as MuHPM, MVIM, and HPM that exist in the literature. Furthermore, the simulations are performed with different time step sizes, and the accuracy is measured at various time intervals. The numerical computations clearly demonstrate that the Galerkin scheme, in contrast to the Runge-Kutta scheme, provides incredibly precise solutions at relatively large time step sizes. A comparison of the solutions reveals that the obtained results through the Galerkin scheme are in fairly good agreement with the RK4 scheme in a given time interval as compared to other conventional schemes. Moreover, having performed various numerical tests for assessing the efficiency and computational cost (in terms of time) of the suggested schemes, it is observed that the Galerkin scheme is noticeably slower than the Runge-Kutta scheme. On the other hand, this work is also concerned with the path tracking and damped oscillatory behaviour of the model with a variable supply rate for the generation of new CD4+ T-cells (based on viral load concentration) and the HIV infection incidence rate. Additionally, we investigate the influence of various physical characteristics by varying their values and analysing them using graphs. The investigations indicate that the lateral system ensured more accurate predictions than the previous model.


---

## Body

## 1. Introduction

Despite recent scientific advancements and substantial health-care initiatives, HIV/AIDS infection remains one of the most devastating diseases in human history. So still, several communities are still seriously impacted by this disease. At the moment, the global spread of HIV infection has an influence on the rising prevalence of other infectious diseases. Infectious diseases are illnesses that can be diagnosed clinically and have a large impact on the human community. They are caused by bacterial pathogens such as microorganisms, retroviruses, fungal, and parasitic agents or by harmful proteins known as prions. The most common are bacterial infectious diseases and tuberculosis, as well as virus-borne HIV and influenza [1, 2]. The hepatitis C virus (HCV) is highly widespread, with almost 71 million cases globally, and coinfections with HIV are also extremely frequent, with approximately 2.3 million coinfections [3]. Infection with HIV destroys the body’s immune system; affects bodily organs such as the brain, kidney, and heart; and leads to death. Unfortunately, there is presently no treatment for this infectious illness. However, there are effective retroviral medications for improving patients’ health complications. HIV is spread through sexual contact, sharing needles, and direct interaction with virus-infected blood or other bodily fluids, as well as mother-to-child transmission during delivery [2]. HIV is the virus that causes a disease known as immunodeficiency syndrome (AIDS), which weakens the body’s capability to fight infection and leaves it vulnerable to attacks from contagious diseases. To improve health-care therapy, a detailed overview of the inflammatory and coagulation processes in HIV infection is essential [4]. In humans, CD4+ T-cells are HIV’s targeted cells, and these cells play an essential role in the immune system’s response. Their depletion can have a wide range of implications, potentially disrupting the immune system’s functionality. Therefore, the reduction in the number of these cells is utilised as a symptomatic and illness stage indicator. In illnesses such as allergies and autoimmune disorders, regulatory T-cells (Treg) play significant roles in maintaining self-tolerance and immunological regulation. These cells also have a role in suppressing efficient immune responses to microbial organisms and cancerous cells [5]. Since the investigation is about the density of CD4+ T-cells, we will refer to T-cells and I-cells as healthy/infected CD4+ T-cells throughout the rest of the manuscript.Numerous investigations have been conducted in order to identify the transmission dynamics of biological processes and infections. Cao et al. [6] presented a comprehensive approach to identifying successful treatment approaches and druggable targets by investigating regulatory interactions of cell phenotypic switching using precisely calculated probability landscapes of response networks. Mathematical formulation of infectious illnesses, such as HIV/AIDS coinfection and acute infection, is identified as a vital tool for understanding the dynamics of disease elimination and/or control [7]. The significance of immunological responses to HIV in modelling was highlighted by [8]. Omondi et al. [9] developed and evaluated a model containing a sex population that examined the HIV infection dynamics in both males and females. Vaidya et al. [10] examined how morphine affects resistivity to target cells, viral dynamics, and viral replication in a steady state, T-cell loss, and the efficiency of antiretroviral treatment (ART). They also looked into how morphine can boost basic reproduction rates, which poses significant challenges to HIV prevention. Vaidya et al. [11] established a viral infection model and investigated the role of antibody responses in reducing viral infection rates during infection. Calshaw and Raun [12] developed a model which consists of populations of T-cells, I-cells, and virus particles. Yusuf et al. [13] introduced a model and analysed the transmission behaviours of HIV infection. Singh [14] suggested a fractional order HIV model, with plasma densities of T-cells, I-cells, and free virus. They examined the dynamics of fractional HIV infection in T-cells as well as the impact of treatment. Thirumalai et al. [15] investigated the fractional order differential equations to analyse the combination of drug treatments for HIV infection. AlShamrani [16] suggested and evaluated an adaptive immune response model for HIV infection. Zaka Ullah and Baleanu [17] presented the SICA model in a fractional structure for the transmission dynamics of HIV/AIDS. Thirumalai et al. [18] studied HIV infection for T-cells using a model based on a system of differential equations. They evaluated and assessed an approximate solution of the mathematical model and presented a count of T-cells, I-cells, and viruses present at any given instant. Khan and Odinsyah [19] investigated HIV dynamics using a fractional operator of the Caputo type. They utilised the data of HIV infection cases in Indonesia from 2006 to 2018 to approximate the HIV model and determine the basic reproduction number. They proceeded by providing a mathematical description of the HIV/AIDS population. Chen and Zheng [20] explored the formation of spatial patterns in a predator-prey model with predator-taxis under homogenous zero-flux boundary conditions. They used the predator-taxis coefficient as the Turing bifurcation’s possible critical value to play its part in constructing the spatial pattern. Chen and Wu [21] explored the diffusive nutrient-microorganism dynamics in a spatially heterogeneous environment. Chen et al. [22] studied the pattern formation and secondary instabilities of the Gierer-Meinhardt model, an activator-inhibitor system, including the Eckhaus instability and zigzag instability. By performing a linear stability analysis at the unique positive equilibrium, conditions on the Hopf bifurcation and Turing instability are discovered. Chen et al. [23] suggested a predator-prey model with a delayed ratio. They explored the system’s Turing-Hopf bifurcation and established amplitude equations using the multiple time scale methodology. Theoretical and numerical simulations are used to find spatiotemporal solutions. Abdel et al. [24] examined the accurate travelling wave solutions of the fractional model of HIV-1 for T-cells and the effects of medication therapy. Naik et al. [25] suggested a Caputo-like fractional order model of HIV-1 that incorporates interactions between cancer cells and shows the chaotic behaviour. Raza et al. [26] developed a delayed model to study HIV dynamics in the community as well as the stability of a susceptible-infectious-immune system with a delay period.In the literature, a variety of traditional analytical and numerical schemes are used to find the approximate solution of various models to describe the HIV dynamics. The homotopy analysis method (HAM) was applied by Ismail and Alomari [27], Ongun [28] performed the Laplace Adomian decomposition technique (LADM), Doğan [29] implemented the multistep LADM, and Yüzbaşı and Karaçayır [30] developed a novel computational approach based on exponential polynomials related to the Galerkin scheme for the HIV infection model. Recently, Sohaib [31] have studied the Legendre wavelet collocation method and the continuous Galerkin-Petrov scheme to obtain the solutions to the HIV model. They made comparison of the outcomes and highlighted the precision and efficacy of the proposed schemes in contrast to other approaches employed to the model. Yüzbaşı et al. [32] presented an exponential approach for numerical solutions to the HIV model based on T-cells. The technique involves exponential polynomials and collocation points. Haq et al. [33] investigated the dynamics of a fractional order smoking cessation model and presented an approximate solution based on the Laplace transformation. Attaullah et al. [34] have applied the Galerkin time discretization scheme to the following HIV models and compared the results with the Runge-Kutta method. Yüzbaşı [35] developed the Bessel collocation approach and implemented it on the HIV T-cell model. The suggested technique comprises of resolving the problem to the system of nonlinear algebraic equations by expanding the approximate solutions with unknown coefficients employing Bessel polynomials. For the numerical investigation, Ahmed et al. [36] proposed an extended temporal model of HIV with medication treatment impact using the backward Euler and Crank-Nicolson techniques. The HIV model of T-cells was explored by Yüzbaşı [37], and an operational technique was suggested to numerically address the model problem. Merdan [38] utilised the homotopy perturbation method (HPM) to examine the numerical outcomes of HIV model. Attaullah [39] approximates the solution of the HIV infection model with full logistic proliferation and variable source terms employing the Galerkin scheme. The modified multivariational iteration method (MVIM) is used by Merdan et al. [40] in order to approximate the suggested model. Multistage homotopy perturbation method (MuHPM) is used by Leal et al. [41] and demonstrated the path tracking behaviour of HIV dynamics.


### 1.1. Fundamental Objectives

The fundamental aim of the present paper is to examine the well-known HIV infection model of T-cells characterized by Leal et al. [41] based on three coupled nonlinear ordinary differential equations. The mentioned model is solved by applying the Galerkin discretization scheme known as the continuous Galerkin-Petrov method, briefly cGP(2) method. This method has superiority to some extent over the other traditional techniques. Each time interval in the cGP(2) method has two unknowns that can be computed by solving a 2×2 block system. The suggested approach is accurate to three orders over the whole time period and even demonstrates super convergence to four orders in discrete time points. To illustrate the reliability of the Galerkin scheme, the model is solved using the Runge-Kutta technique of order four, shortened as the RK4 method. Our second aim is to compare the solutions of the Galerkin method with the findings of other techniques, e.g., MuHPM [42], MVIM [40], and HPM [38]. All the findings are visualised through different graphs. The comparative study was performed to validate the accuracy and computational cost (in terms of time) of the Galerkin and RK4 methods for the same and different time step sizes. Therefore, to having an idea about the accuracy and reliability of the approximate solutions obtained from the proposed technique, the well-known classical RK4 method is also implemented. Estimate the absolute errors between the results of both schemes to better understand and authenticate our proposed method’s applicability in practical problems. The assessment of the absolute values of errors between the Galerkin scheme and RK4 scheme solutions shows that the solutions provided are more accurate in comparison to the solutions obtained through conventional schemes. The graphical and tabular outcomes demonstrated that the proposed scheme is very accurate and achieves highly accurate results at a larger step size in comparison to the existing techniques employed for the mentioned model. The Galerkin scheme is an adaptive scheme that achieves the same accuracy at a larger step size at a lower cost. The detailed analysis of the aforementioned model demonstrates that the Galerkin scheme is more authentic and accurate than the previous approaches employed for the model. Our third aim is to develop it by the insertion of a new dynamic supply term dependent on the viral density for the generation of new T-cells and a mass action term in the HIV model. Leal et al. [41] suggested the HIV infection model with a stable parameter for the production of new T-cells from the thymus. However, HIV has the potential to invade T-cells in the thymus, which might affect the generation of new cells from the thymus. Therefore, the model discussed in [41] is extended by introducing a variable source term (“sV=0.5Is+sI/I+V” discussed in [43, 44]) based on viral density. In addition, the mass action term kVT is introduce to explain how infection occurs when the virus interacts with T-cells, causing free virus to be eradicated at a rate of −kVT, where k indicates the infection rate. Moreover, the variation of different clinical parameters involved in the model is observed by changing their values. A computer program written in mathematical software is used to perform these computations.


### 1.2. The Outline of the Paper

The rest of the article is structured as follows: Section4 introduces the fundamental formulation of the HIV model. Section 3.1 provides the Galerkin method implemented for the model. Numerical results, behaviours of different parameters, and comparison of the solutions of the Galerkin method and RK4 method with other classical techniques applied to the model are discussed in Section 4.1. Section 5 gives the conclusion and summary of the article.


## 1.1. Fundamental Objectives

The fundamental aim of the present paper is to examine the well-known HIV infection model of T-cells characterized by Leal et al. [41] based on three coupled nonlinear ordinary differential equations. The mentioned model is solved by applying the Galerkin discretization scheme known as the continuous Galerkin-Petrov method, briefly cGP(2) method. This method has superiority to some extent over the other traditional techniques. Each time interval in the cGP(2) method has two unknowns that can be computed by solving a 2×2 block system. The suggested approach is accurate to three orders over the whole time period and even demonstrates super convergence to four orders in discrete time points. To illustrate the reliability of the Galerkin scheme, the model is solved using the Runge-Kutta technique of order four, shortened as the RK4 method. Our second aim is to compare the solutions of the Galerkin method with the findings of other techniques, e.g., MuHPM [42], MVIM [40], and HPM [38]. All the findings are visualised through different graphs. The comparative study was performed to validate the accuracy and computational cost (in terms of time) of the Galerkin and RK4 methods for the same and different time step sizes. Therefore, to having an idea about the accuracy and reliability of the approximate solutions obtained from the proposed technique, the well-known classical RK4 method is also implemented. Estimate the absolute errors between the results of both schemes to better understand and authenticate our proposed method’s applicability in practical problems. The assessment of the absolute values of errors between the Galerkin scheme and RK4 scheme solutions shows that the solutions provided are more accurate in comparison to the solutions obtained through conventional schemes. The graphical and tabular outcomes demonstrated that the proposed scheme is very accurate and achieves highly accurate results at a larger step size in comparison to the existing techniques employed for the mentioned model. The Galerkin scheme is an adaptive scheme that achieves the same accuracy at a larger step size at a lower cost. The detailed analysis of the aforementioned model demonstrates that the Galerkin scheme is more authentic and accurate than the previous approaches employed for the model. Our third aim is to develop it by the insertion of a new dynamic supply term dependent on the viral density for the generation of new T-cells and a mass action term in the HIV model. Leal et al. [41] suggested the HIV infection model with a stable parameter for the production of new T-cells from the thymus. However, HIV has the potential to invade T-cells in the thymus, which might affect the generation of new cells from the thymus. Therefore, the model discussed in [41] is extended by introducing a variable source term (“sV=0.5Is+sI/I+V” discussed in [43, 44]) based on viral density. In addition, the mass action term kVT is introduce to explain how infection occurs when the virus interacts with T-cells, causing free virus to be eradicated at a rate of −kVT, where k indicates the infection rate. Moreover, the variation of different clinical parameters involved in the model is observed by changing their values. A computer program written in mathematical software is used to perform these computations.


## 1.2. The Outline of the Paper

The rest of the article is structured as follows: Section4 introduces the fundamental formulation of the HIV model. Section 3.1 provides the Galerkin method implemented for the model. Numerical results, behaviours of different parameters, and comparison of the solutions of the Galerkin method and RK4 method with other classical techniques applied to the model are discussed in Section 4.1. Section 5 gives the conclusion and summary of the article.


## 2. Basic Mathematical Structure of the HIV Infection Model

This section concerns the model suggested by Leal et al. [41]. Our first aim in this study is to apply the continuous Galerkin-Petrov time discretization scheme [31, 34] to the following HIV model and to compare outcomes with the methods existing in literature to validate the accuracy and efficiency of the scheme and the MATLAB code. The model is as follows:
(1)dTdt=s−μTT+rT1−T+ITmax,dIdt=kVT−μII,dVdt=NμII−μVV+kVT,where the state variables Tt,It, and Vt show the density of T-cells, I-cells, and the HIV virus in the blood, respectively. The diagrammatic illustration of the HIV model and the structure of the HIV virus are given in Figure 1. The initial conditions of state variables, different parameters, and constants with their explanations are provided in Table 1.Figure 1
The diagrammatic description of the model and HIV virus.
(a)(b)Table 1
Description of state variables, parameters, and their values [41] (units: day-1 mm-3).
DescriptionsValuesDependent variablesT0Population of T-cells0.10I0Population of I-cells0.0V0Population of HIV particles0.10Parameters and constantssThe supply rate of T-cells0.1rThe growth rate of the T-cell population3μTThe killing rate of T-cells0.02μIThe killing rate of I-cells0.3μVThe destruction rate of free virus2.4kThe infection rate of T-cells by free virus2.4×10−5NThe density of virus released by I-cellsVariesTmaxThe maximum density of T-cells1500


## 3. The Numerical Schemes

This section discusses the numerical schemes used to solve the above-mentioned model. Nowadays, the Galerkin approach has been used effectively to handle a wide range of complicated problems in engineering and science, seen for example [31, 45–49]. We implemented the technique on the HIV model discussed in [41]. In addition, we employ the conventional Runge-Kutta approach to assess the precision of the Galerkin scheme findings.


### 3.1. The Continuous Galerkin-Petrov (cGP) Technique

The system of ODE’s for the considered model can be written follows:Findu~:0,tmax⟶V=ℝd like as follows:
(2)dtu~t=Ft,u~t∀t∈I,u~0=u0,where dt shows the time derivative of u~t, I=0,T is the total interval, and u~t=u~1t,u~2t,u~3t∈V⇒u~0=u~10,u~20,u~30∈V are the initial values of u~t at t=0. We also assume that u~1t,u~2t,u~3t=Tt,It,Vt which implies that u~10,u~20,u~30=T0,I0,V0. The function F=f1,f2,f3 is nonlinear and is described as F:I×K⟶K.Theweak formulation (see [31, 34, 45–49] for explanations) of a the problem (2) is follows: find u~∈X such that u~0=u~0 and
(3)∫Idtu~t,vtdt=∫IFt,u~t,νtdtforallν∈Y,where X and Y represent the solution and test space, respectively, to explain the time discretization of a variational type problem (2).To characterize the functiont⟶u~t, we describe the space CI,K=C0I,K that is the space of continuous functions u~: I⟶K equipped with the norm as follows:
(4)u~CI,K=supt∈Iu~K.We will use the spaceL2I,K as the space of discontinuous functions which is given by
(5)L2I,K=u~:I⟶K:u~L2I,K=∫Iu~K2dt<∞1/2.In time discretization, we split the intervalsI into N subintervals Iτ=tτ−1,tτ, where τ=1,⋯,N and 0=t0<t1<t2<⋯<tN−1<tN=T. The parameter j indicates the time discretization parameter, as well as the maximum time step size j=max1≤τ≤Njτ, where jτ=tτ−tτ−1 is the length of nth time interval Iτ. The following set of time intervals Mj=I1,⋯,IN will be called the time mesh. We find out the solution u~: I⟶K on each time interval Iτ by a function u~j: I⟶K which is a piecewise polynomial of some order l w.rt time. The time-discrete solution space for u~j is Xjl⊂X and is defined by
(6)Xjl=u~∈CI,K:u~Iτ∈ℙlIτ,KforallIτ∈Mj,where
(7)ℙlIτ,K=u~:Iτ⟶K:u~t=∑s=0lUsts,forallt∈Iτ,Us∈K,∀s.The discrete test space foru~ȷ is Yȷl⊂Y and is defined by
(8)Yjk=v∈L2I,K:vIτ∈ℙk−1I,K∀Iτ∈Mj,which is composed ofl−1 piecewise polynomials (see [31, 34, 45–49] for details) and is discontinuous at the time step end nodes. We multiply equation (2) by the test function νȷ∈Yȷk and integrate over the interval I. We get the discrete-time problem: find u~ȷ∈Xȷk such that u~ȷ0=0 and
(9)∫Iu~ȷ′t,νȷtdt=∫IFt,u~ȷt,νȷtdt∀νȷ∈Yȷl,where 〈·and·〉 represent the usual inner product in L2I,K. This discretization is known as the exact continuous Galerkin-Petrov method or simply the “exact cGP(l) scheme” of order k. The Galerkin-Petrov name is due to the fact that the solution space Xȷl is different from the test space Yȷl. The term “exact” denotes that the time integral on the right-hand side of equation (9) is determined exactly. Because the discrete test space Yȷl is discontinuous, equation (9) could be computed by a time marching technique in which local problems on the interval are handled successively. Therefore, we select the test function νȷt=νφt with arbitrary time independent ν∈K and a scalar function φ:I⟶ℝ which is zero on IIτ and a polynomial of order less than or equal to l−1 on the time interval Iτ=tτ−1,tτ. Then, we get from (9) the Iτ problem of the exact cGP scheme of order l: find u~ȷIτ∈ℙlIτ,K such that
(10)∫Iτdtu~ȷt,νφtdt=∫IτFt,u~ȷt,νφtdt∀ν∈K∀φ∈ℙl−1Iτ,with the initial conditionu~ȷIτtτ−1=u~ȷIτ−1tτ−1 for τ≥2 and u~ȷIτtτ−1=u~0 for τ=1.In case of a nonlinear functionF·,·, we need to calculate the integrals numerically on the right-hand side of equation (10). The l+1-point Gauβ-Lobatto formula is exact if the function to be integrated has a polynomial of degree less than or equal to 2l−1. As a result, this formula is applied to the integral on the left-hand side of (10) which will give the exact value. Then, the “Iτ problem of the numerically integrated cGPl method” is as follows: find u~ȷIτ∈ℙlIτ,K such that u~ȷtτ−1=u~τ−1,
(11)∑s=0lw^sdtu~ȷtτ,sφtτ,s=∑s=0kw^sFtτ,s,u~ȷtτ,sφtτ,s∀φ∈ℙl−1Iτ,where ws are the weights and t^∈−1,1,s=0,1,2,3,⋯,l represent the nodes on the reference interval. To find u~ȷ on each time interval Iτ, we use a polynomial ansatz to illustrate it as follows:
(12)u~ȷt=∑s=0kUτsϕτ,st∀t∈Iτ,where the coefficients Uτs are the components of K and the functions ϕτ,s∈ℙlIτ are the Lagrange basis functions (see [31, 34, 45–49] for more details) with respect to l+1 suitable nodal points tτ,s∈Iτ satisfying conditions mentioned below:
(13)ϕτ,stτ,r=δr,s,r,s=0,1,2,⋯,l,where δr,s is the Kronecker delta that is defined as
(14)δr,s=1ifr=s,0ifr≠s.Like in [50], the tτ,s have been defined as the quadrature points of l+1-point Gauβ-Lobatto formula (for detail information, see [31, 34, 45–49]) on the interval Iτ. For the selection of initial conditions, we can set tτ,0=tτ−1 which denotes the initial condition for equation (10)
(15)Uτ0=u~ȷIτ−1ifτ≥2,forτ=1⇒Uτ0=u~0.We define the basis functionsϕτ,s∈ℙlIτ via the affine reference transformation (see [31, 34, 45–49] for detail explanations) T¯:I^⟶Iτ, where I^=−1,1 and
(16)t=T¯t^=tτ+tτ−12+ȷτ2t^∈Iτ∀t^∈Iτ,τ=1,2,3⋯,N.Letϕs^∈ℙlI^, s=0,1,⋯,l, demonstrate the basis functions that meet the requirements
(17)ϕ^st^r=δr,s,r,s=0,1,⋯,l,where t^0=−1 and t^r, r=1,2,⋯,l, are the quadrature points for the reference interval I^. Then, we define the basis functions on the given interval Iτ by the mapping
(18)ϕτ,st=ϕs^t^witht^=T¯τ−1t=2ȷτt+tτ−1−tτ2∈I^.Likewise, the test basis functionsφτ,r are described by the suitable reference basis functions φ^∈ℙl−1I^, i.e.,
(19)φτ,rt=φr^T¯τ−1t∀t∈Iτ,r=1,2,⋯,l.From the representation (12), we get for dtu~ȷ(20)dtu~ȷt=∑s=0kUτsϕτ,s′t∀t∈Iτ.By putting (20) in (10), we get
(21)∫Iτdtu~ȷt,νφtdt=∫Iτ∑s=0lUτs,νϕs′tφt.The integral is now transformed into the reference intervalhatI and computed using the l+1-point Gauβ-Lobatto quadrature formula which leads, for each test basis function φ∈ℙl−1 and for all ν∈K,
(22)∫I^τ∑s=0lUτs,νϕ^s′t^φ^t^dt^=∫I^τFωτt^,∑s=0lUτst^,νφ^t^dt^∀ν∈K,⇒∑μ=0lw^μ∑s=0lUτs,νϕ^s′t^μφ^t^μ=∑μ=0lw^μFωτt^μ,∑s=0lUτsϕ^st^μ,νφ^t^μ,where w^μ are the weights and t^μ∈−1,1 are the points of integration with t^0=−1 and t^l=1. If we choose the test functions φτ,i∈ℙl−1Iτ such that
(23)φ^tμ^=w∧−1δr,μr,μ=1,2,⋯,l.Now find the coefficients that are unknownUτs∈K for s=1,⋯,l,(24)∑s=0lαr,sUτs=ȷτ2Ftτ,r,Uτs+βrFtτ,0,Uτs∀i=1,2,⋯,l,where Uτs=Uτ−1s for τ>1 and U10=u~0 for τ=1 and
(25)αr,s=ϕ^s′t^r+βrϕ^s′t0,βr=w^0φ^rt^0.We will discuss thecGP(k) method for the cases l=1 and l=2 in the following subsections.


#### 3.1.1. The cGP(1) Method

We used the two point Gauβ-Lobatto formula with tτ,0=tτ−1,tτ,1=tτ and weights w^0=w^1=1 which gives the well-known trapezoidal rule. We obtain α1,0=−1,α1,1=1, and β1=1. For the single coefficient Uτ1=u~ȷtτ∈K, the problem leads to the following block equation:
(26)α1,1Uτ1−ατ,0Uτ0=ȷτ2Ftτ,Uτ1+Ftτ−1,Uτ0.


#### 3.1.2. The cGP(2) Method

Three-point Gauβ-Lobatto formula (Simpson’s rule) is used to define the quadratic basis functions with weights w^0=w^2=1/3, w^1=4/3, t^0=−1, t^1=0, and t^2=1. Then, we get
(27)αr,s=−541142−42,βr=12−1,r=1,2,s=0,1,2.Thus, the system to be solved forUτ1,Uτ2∈K from the known Uτ0=Uτ−12 becomes
(28)α1,1Uτ2=−α1,0Uτ0+ȷ2τFtτ,1,Uτ1+β1Ftτ,0,Uτ0,α2,1Uτ1+α2,2Uτ2=−α2,0Uτ0+ȷτ2Ftτ,2,Uτ2+β2Ftτ,0,Uτ0,where Uτ0 indicate the initial condition at the current time interval.


### 3.2. Runge-Kutta Method

This is a well-known technique of order four, established by Kutta [51] and extensively utilised for initial value problems (see [52, 53] for more details).


### 3.3. Comparison and Computation Analysis

This section demonstrates the significance of the Galerkin and RK4 methods to the system of ODE’s 4.1–4.3 in the aforementioned model. The initial conditions and parameter values are given in Table1. In order to determine the effectiveness of the suggested technique, we contrasted the numerical and graphical solutions achieved by the Galerkin scheme with those obtained by the RK4 method and other conventional methods, i.e., MuHPM [41], MVIM [40], and HPM [38], provided in Tables 2–4 for Tt, It, and Vt. After the evaluation, it was revealed that the proposed approach produced more precise findings than the solutions obtained through conventional approaches. Furthermore, numerical computations are carried out using the Galerkin scheme and the RK4 scheme with the same and different time step sizes. The absolute errors between the outcomes were assessed using both approaches, as demonstrated in Tables 5 and 6. From comparison, it is clearly visible that the Galerkin scheme achieves more precise findings at a larger step size as compared to the RK4 technique. Figure 2 shows the results of both methods for t∈0,1 and t∈0,70 with ȷ=0.1 steps. The figures showed that the results overlapped over each other throughout the time period. In addition, the analysis of the numerical costs for different step sizes in terms of time shows that the Galerkin method is more time-consuming as in comparison to the RK4 method is presented in Table 7 and Figure 3(f). The polar and bar graphs illustrate the absolute errors between the solutions of the Galerkin and RK4 approaches depicted in Figure 4. The mesh graphs for the solutions of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] are visualised in 3. Overall, we may conclude that the numerical approach reported in this article can be relied on to produce reasonably flexible and accurate results when applied to comparable problems.Table 2
The findings of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] for Tt.
tiMuHPM [41]MVIM [40]HPM [38]RK4 methodPresent solutions0.00.100000000.100000000.1000000000.10000000000000.10000000000000.10.146358950.146359100.1463589000.14635633342150.14635849282220.20.208807720.208808100.2087991000.20880067887670.20880649648410.40.406239230.406240800.4056066000.40621367496460.40623478431240.60.764420340.764428700.7564485000.76435081453700.76440824443910.81.414038301.414094001.3642150001.41387024891631.41400906112661.02.591575702.591921002.3786790002.59119519036622.59150945893731.24.723924184.725783004.0065120004.72309837340384.72378033501621.48.578327308.587223006.5213030008.57658959800018.57802431662681.615.522792515.561670010.2735700015.51922813062515.5221702037551.827.961354328.119610015.7007900027.95421977630627.9601062983772.050.008070750.614770023.3374000049.99418511292750.0056374736263.0605.328400758.304100121.4166000605.19715873002605.309626004254.01387.07811-38782.630428.18300001387.03656989361387.07346477905.01484.41661-1.2822E+71179.6330001484.41297411281484.416218519410578.768500-7.535E+1831121.36000578.01860976666578.803221775333042.9620000-7.954E+656710438.00042.98826240709642.9794781819207051.1260000-8.86E+1594.49161E+0850.73161364454350.729213455592Table 3
The findings of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] for It.
tiMuHPM [41]MVIM [40]HPM [38]RK4 methodPresent solutions0.00.000000000.000000000.000000000.00000000000000.00000000000000.10.000002860.000002860.000002860.00000286453320.00000286484710.20.000006030.000006030.000006030.00000603187890.00000603254640.40.000013160.000013150.000013150.00001315648600.00001315797850.60.000021220.000021220.000021220.00002122067760.00002122315970.80.000030180.000030170.000029940.00003017281020.00003017646611.00.000040040.000040020.000039180.00004003141580.00004003645811.20.000050880.000050840.000048440.00005087001740.00005087669631.40.000062830.000062700.000056970.00006281620700.00006282482191.60.000076080.000075740.000063540.00007606842280.00007607934211.80.000090960.000090110.000066490.00009094185170.00009095554992.00.000107990.000105970.000063460.00010796867950.00010798581413.00.000305680.00021559-0.00018160.00030559579170.00030566741484.00.002109200.00040533-0.00148370.00210841864440.00210907783795.00.020119310.00074045-0.00550380.02011076952760.020117982587410.0828.5948000.01387534-0.2285212828.38422383941828.5668668426630.0775.2937001625.28280-64.142730775.29918583941775.2857747503970.0594.7040002.2292E+13-4621.4330594.00643606878594.00643962536Table 4
The findings of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] for Vt.
tiMuHPM [41]MVIM [40]HPM [38]RK4 methodPresent solutions0.00.100000000.100000000.100000000.10000000000000.10000000000000.10.078663150.078663180.078663150.07866381455020.07866326355950.20.061879800.061879910.061878250.06188084740160.06187998051660.40.038294840.038295960.038199470.03829613043710.03829505758290.60.023704500.023710290.022681580.02370570312840.02370470746760.80.014680330.014700420.009255110.01468131432290.01468049325491.00.009100820.00915723-0.01048470.00910157907680.00910094475311.20.005653260.00579375-0.04982590.00565382162980.00565335183931.40.003525410.00385124-0.12948780.00352581552480.00352547818281.60.002214750.00293882-0.28016880.00221503826520.00221480224191.80.001410460.00297863-0.54509450.00141065416330.00141049332042.00.000920390.00426490-0.98256560.00092052056450.00092041462753.00.000329410.13649810-9.17492000.00032936511840.00032939832824.00.001410905.40689000-43.5530800.00141038094990.00141082558615.00.01292825214.291400-143.651800.01292275834970.012927397503710.0668.7895902.095E+10-5456.7680668.57158158359668.7579298116730.0983.9820001.916E+42-150636600983.92573498469983.9243818601770.0769.4900001.602E+106-1.0827E+8768.96288598031768.95708309906Table 5
The absolute errors between the findings of the Galerkin and RK4 schemes forTt, It, and Vt with same step size and t∈0,70.
cGP2j=0.01−RK4j=0.01tiTtItVt0.00.00000000000000E-100.00000000000000E-140.00000000000000E-110.12.84775258929670E-103.46123001238724E-144.45784936742299E-110.27.67201135953854E-107.50830062186403E-147.01317059981221E-110.31.55009055591293E-091.21598593011158E-138.27457546925814E-110.42.78372330742016E-091.74456149260862E-138.67752675270950E-110.54.68632355143939E-092.34046319786763E-138.53053669369608E-110.67.57295037789874E-093.00841803183537E-138.04949266963728E-110.71.18961085426861E-083.75390284178614E-137.38314721859812E-110.81.83025996580710E-084.58311269117091E-136.63192608468810E-110.92.77127207848338E-085.50296584017097E-135.86177772987639E-111.28.98785703640215E-088.88763672032833E-133.77432252118837E-111.41.89002486550294E-071.17644883093106E-122.69146562675848E-111.63.87287119707480E-071.52607589651126E-121.85752377021564E-111.87.73954734967219E-071.95361702347036E-121.23150862337690E-112.01.50285070077416E-062.48521210934657E-127.65154150976499E-123.01.44064168807745E-051.06163661303844E-116.34902492992889E-124.04.45008413407777E-069.61197484453114E-116.36393082328879E-115.03.44141881214455E-071.02496256620954E-096.58341856724087E-10102.90578410613307E-052.45547904569321E-052.57475288663045E-05301.43229695481750E-062.60882575275900E-065.73396846448304E-07704.48719525536490E-071.54975691657455E-061.31251681523281E-06Table 6
The absolute errors forTt, It, and Vt between the Galerkin and RK4 techniques for t∈0,70 with different step sizes.
cGP2j=0.01−RK4j=0.0001tiTtItVt0.00.00000000000000E-110.00000000000000E-150.00000000000000E-120.15.86432846727547E-117.12909557572966E-158.69847249784783E-120.21.57959229030169E-101.54760669814828E-141.36827840724330E-110.33.19075987853523E-102.50826784472732E-141.61410537713458E-110.45.72855207714440E-103.60121817349207E-141.69235597757833E-110.59.64069601927520E-104.83457641053866E-141.66324766459613E-110.61.55727908346393E-096.21799635977353E-141.56892346336868E-110.72.44506304092340E-097.76247081478616E-141.43842576738606E-110.83.75948450148655E-099.48020382628059E-141.29135296206284E-110.95.68795943678424E-091.13844887293310E-131.14057773642484E-111.08.49506598399330E-091.34896515615809E-139.94224494399454E-121.21.83776105444622E-081.83649035106863E-137.31735158576718E-121.43.84712617318428E-082.42403387401731E-135.19390173231438E-121.67.82565194867857E-083.12689110715247E-133.55596108114753E-121.81.54546054176308E-073.96344591803606E-132.32396979264748E-122.02.94468271988535E-074.95918884842010E-131.40495156741166E-123.02.23191591430805E-061.76651054895061E-121.14225539834048E-124.05.47712943443912E-071.65905887118678E-111.07716444271155E-115.04.52982931165025E-081.77432183529058E-101.13838096935104E-10104.84750069063011E-064.06695653509814E-064.55263102594472E-06302.30107175980265E-074.45262458015350E-071.09476786747109E-07707.22985902257278E-082.61271452473011E-072.25629946726258E-07Figure 2
The geometrical comparison between the findings of the Galerkin and RK4 schemes forTt, It, and Vt for t∈0,70 and t∈0,1 with j=0.1.
(a)(b)(c)(d)(e)(f)Table 7
The numerical cost in terms of the CPU time (in seconds) for the cGP(2) scheme and RK4 scheme fort∈0,70.
CPU time (in seconds)1/jcGP(2) methodRK4 method700.8513710.0584631000.9589540.0540762001.6077080.0677535002.6876630.08498110004.4738710.12477515005.2455370.19959320006.0210620.234840500010.2441010.5153861000017.5105890.8591492000031.6600671.945279Figure 3
The mesh graphs of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] for the model and the computation cost in terms of CPU time (f).
(a)(b)(c)(d)(e)(f)Figure 4
The polar and bar graphs of errors between the findings of the Galerkin and RK4 schemes and bar graphs for the Galerkin and RK4 schemes forTt, It, and Vt.
(a)(b)(c)(d)(e)(f)


## 3.1. The Continuous Galerkin-Petrov (cGP) Technique

The system of ODE’s for the considered model can be written follows:Findu~:0,tmax⟶V=ℝd like as follows:
(2)dtu~t=Ft,u~t∀t∈I,u~0=u0,where dt shows the time derivative of u~t, I=0,T is the total interval, and u~t=u~1t,u~2t,u~3t∈V⇒u~0=u~10,u~20,u~30∈V are the initial values of u~t at t=0. We also assume that u~1t,u~2t,u~3t=Tt,It,Vt which implies that u~10,u~20,u~30=T0,I0,V0. The function F=f1,f2,f3 is nonlinear and is described as F:I×K⟶K.Theweak formulation (see [31, 34, 45–49] for explanations) of a the problem (2) is follows: find u~∈X such that u~0=u~0 and
(3)∫Idtu~t,vtdt=∫IFt,u~t,νtdtforallν∈Y,where X and Y represent the solution and test space, respectively, to explain the time discretization of a variational type problem (2).To characterize the functiont⟶u~t, we describe the space CI,K=C0I,K that is the space of continuous functions u~: I⟶K equipped with the norm as follows:
(4)u~CI,K=supt∈Iu~K.We will use the spaceL2I,K as the space of discontinuous functions which is given by
(5)L2I,K=u~:I⟶K:u~L2I,K=∫Iu~K2dt<∞1/2.In time discretization, we split the intervalsI into N subintervals Iτ=tτ−1,tτ, where τ=1,⋯,N and 0=t0<t1<t2<⋯<tN−1<tN=T. The parameter j indicates the time discretization parameter, as well as the maximum time step size j=max1≤τ≤Njτ, where jτ=tτ−tτ−1 is the length of nth time interval Iτ. The following set of time intervals Mj=I1,⋯,IN will be called the time mesh. We find out the solution u~: I⟶K on each time interval Iτ by a function u~j: I⟶K which is a piecewise polynomial of some order l w.rt time. The time-discrete solution space for u~j is Xjl⊂X and is defined by
(6)Xjl=u~∈CI,K:u~Iτ∈ℙlIτ,KforallIτ∈Mj,where
(7)ℙlIτ,K=u~:Iτ⟶K:u~t=∑s=0lUsts,forallt∈Iτ,Us∈K,∀s.The discrete test space foru~ȷ is Yȷl⊂Y and is defined by
(8)Yjk=v∈L2I,K:vIτ∈ℙk−1I,K∀Iτ∈Mj,which is composed ofl−1 piecewise polynomials (see [31, 34, 45–49] for details) and is discontinuous at the time step end nodes. We multiply equation (2) by the test function νȷ∈Yȷk and integrate over the interval I. We get the discrete-time problem: find u~ȷ∈Xȷk such that u~ȷ0=0 and
(9)∫Iu~ȷ′t,νȷtdt=∫IFt,u~ȷt,νȷtdt∀νȷ∈Yȷl,where 〈·and·〉 represent the usual inner product in L2I,K. This discretization is known as the exact continuous Galerkin-Petrov method or simply the “exact cGP(l) scheme” of order k. The Galerkin-Petrov name is due to the fact that the solution space Xȷl is different from the test space Yȷl. The term “exact” denotes that the time integral on the right-hand side of equation (9) is determined exactly. Because the discrete test space Yȷl is discontinuous, equation (9) could be computed by a time marching technique in which local problems on the interval are handled successively. Therefore, we select the test function νȷt=νφt with arbitrary time independent ν∈K and a scalar function φ:I⟶ℝ which is zero on IIτ and a polynomial of order less than or equal to l−1 on the time interval Iτ=tτ−1,tτ. Then, we get from (9) the Iτ problem of the exact cGP scheme of order l: find u~ȷIτ∈ℙlIτ,K such that
(10)∫Iτdtu~ȷt,νφtdt=∫IτFt,u~ȷt,νφtdt∀ν∈K∀φ∈ℙl−1Iτ,with the initial conditionu~ȷIτtτ−1=u~ȷIτ−1tτ−1 for τ≥2 and u~ȷIτtτ−1=u~0 for τ=1.In case of a nonlinear functionF·,·, we need to calculate the integrals numerically on the right-hand side of equation (10). The l+1-point Gauβ-Lobatto formula is exact if the function to be integrated has a polynomial of degree less than or equal to 2l−1. As a result, this formula is applied to the integral on the left-hand side of (10) which will give the exact value. Then, the “Iτ problem of the numerically integrated cGPl method” is as follows: find u~ȷIτ∈ℙlIτ,K such that u~ȷtτ−1=u~τ−1,
(11)∑s=0lw^sdtu~ȷtτ,sφtτ,s=∑s=0kw^sFtτ,s,u~ȷtτ,sφtτ,s∀φ∈ℙl−1Iτ,where ws are the weights and t^∈−1,1,s=0,1,2,3,⋯,l represent the nodes on the reference interval. To find u~ȷ on each time interval Iτ, we use a polynomial ansatz to illustrate it as follows:
(12)u~ȷt=∑s=0kUτsϕτ,st∀t∈Iτ,where the coefficients Uτs are the components of K and the functions ϕτ,s∈ℙlIτ are the Lagrange basis functions (see [31, 34, 45–49] for more details) with respect to l+1 suitable nodal points tτ,s∈Iτ satisfying conditions mentioned below:
(13)ϕτ,stτ,r=δr,s,r,s=0,1,2,⋯,l,where δr,s is the Kronecker delta that is defined as
(14)δr,s=1ifr=s,0ifr≠s.Like in [50], the tτ,s have been defined as the quadrature points of l+1-point Gauβ-Lobatto formula (for detail information, see [31, 34, 45–49]) on the interval Iτ. For the selection of initial conditions, we can set tτ,0=tτ−1 which denotes the initial condition for equation (10)
(15)Uτ0=u~ȷIτ−1ifτ≥2,forτ=1⇒Uτ0=u~0.We define the basis functionsϕτ,s∈ℙlIτ via the affine reference transformation (see [31, 34, 45–49] for detail explanations) T¯:I^⟶Iτ, where I^=−1,1 and
(16)t=T¯t^=tτ+tτ−12+ȷτ2t^∈Iτ∀t^∈Iτ,τ=1,2,3⋯,N.Letϕs^∈ℙlI^, s=0,1,⋯,l, demonstrate the basis functions that meet the requirements
(17)ϕ^st^r=δr,s,r,s=0,1,⋯,l,where t^0=−1 and t^r, r=1,2,⋯,l, are the quadrature points for the reference interval I^. Then, we define the basis functions on the given interval Iτ by the mapping
(18)ϕτ,st=ϕs^t^witht^=T¯τ−1t=2ȷτt+tτ−1−tτ2∈I^.Likewise, the test basis functionsφτ,r are described by the suitable reference basis functions φ^∈ℙl−1I^, i.e.,
(19)φτ,rt=φr^T¯τ−1t∀t∈Iτ,r=1,2,⋯,l.From the representation (12), we get for dtu~ȷ(20)dtu~ȷt=∑s=0kUτsϕτ,s′t∀t∈Iτ.By putting (20) in (10), we get
(21)∫Iτdtu~ȷt,νφtdt=∫Iτ∑s=0lUτs,νϕs′tφt.The integral is now transformed into the reference intervalhatI and computed using the l+1-point Gauβ-Lobatto quadrature formula which leads, for each test basis function φ∈ℙl−1 and for all ν∈K,
(22)∫I^τ∑s=0lUτs,νϕ^s′t^φ^t^dt^=∫I^τFωτt^,∑s=0lUτst^,νφ^t^dt^∀ν∈K,⇒∑μ=0lw^μ∑s=0lUτs,νϕ^s′t^μφ^t^μ=∑μ=0lw^μFωτt^μ,∑s=0lUτsϕ^st^μ,νφ^t^μ,where w^μ are the weights and t^μ∈−1,1 are the points of integration with t^0=−1 and t^l=1. If we choose the test functions φτ,i∈ℙl−1Iτ such that
(23)φ^tμ^=w∧−1δr,μr,μ=1,2,⋯,l.Now find the coefficients that are unknownUτs∈K for s=1,⋯,l,(24)∑s=0lαr,sUτs=ȷτ2Ftτ,r,Uτs+βrFtτ,0,Uτs∀i=1,2,⋯,l,where Uτs=Uτ−1s for τ>1 and U10=u~0 for τ=1 and
(25)αr,s=ϕ^s′t^r+βrϕ^s′t0,βr=w^0φ^rt^0.We will discuss thecGP(k) method for the cases l=1 and l=2 in the following subsections.


### 3.1.1. The cGP(1) Method

We used the two point Gauβ-Lobatto formula with tτ,0=tτ−1,tτ,1=tτ and weights w^0=w^1=1 which gives the well-known trapezoidal rule. We obtain α1,0=−1,α1,1=1, and β1=1. For the single coefficient Uτ1=u~ȷtτ∈K, the problem leads to the following block equation:
(26)α1,1Uτ1−ατ,0Uτ0=ȷτ2Ftτ,Uτ1+Ftτ−1,Uτ0.


### 3.1.2. The cGP(2) Method

Three-point Gauβ-Lobatto formula (Simpson’s rule) is used to define the quadratic basis functions with weights w^0=w^2=1/3, w^1=4/3, t^0=−1, t^1=0, and t^2=1. Then, we get
(27)αr,s=−541142−42,βr=12−1,r=1,2,s=0,1,2.Thus, the system to be solved forUτ1,Uτ2∈K from the known Uτ0=Uτ−12 becomes
(28)α1,1Uτ2=−α1,0Uτ0+ȷ2τFtτ,1,Uτ1+β1Ftτ,0,Uτ0,α2,1Uτ1+α2,2Uτ2=−α2,0Uτ0+ȷτ2Ftτ,2,Uτ2+β2Ftτ,0,Uτ0,where Uτ0 indicate the initial condition at the current time interval.


## 3.1.1. The cGP(1) Method

We used the two point Gauβ-Lobatto formula with tτ,0=tτ−1,tτ,1=tτ and weights w^0=w^1=1 which gives the well-known trapezoidal rule. We obtain α1,0=−1,α1,1=1, and β1=1. For the single coefficient Uτ1=u~ȷtτ∈K, the problem leads to the following block equation:
(26)α1,1Uτ1−ατ,0Uτ0=ȷτ2Ftτ,Uτ1+Ftτ−1,Uτ0.


## 3.1.2. The cGP(2) Method

Three-point Gauβ-Lobatto formula (Simpson’s rule) is used to define the quadratic basis functions with weights w^0=w^2=1/3, w^1=4/3, t^0=−1, t^1=0, and t^2=1. Then, we get
(27)αr,s=−541142−42,βr=12−1,r=1,2,s=0,1,2.Thus, the system to be solved forUτ1,Uτ2∈K from the known Uτ0=Uτ−12 becomes
(28)α1,1Uτ2=−α1,0Uτ0+ȷ2τFtτ,1,Uτ1+β1Ftτ,0,Uτ0,α2,1Uτ1+α2,2Uτ2=−α2,0Uτ0+ȷτ2Ftτ,2,Uτ2+β2Ftτ,0,Uτ0,where Uτ0 indicate the initial condition at the current time interval.


## 3.2. Runge-Kutta Method

This is a well-known technique of order four, established by Kutta [51] and extensively utilised for initial value problems (see [52, 53] for more details).


## 3.3. Comparison and Computation Analysis

This section demonstrates the significance of the Galerkin and RK4 methods to the system of ODE’s 4.1–4.3 in the aforementioned model. The initial conditions and parameter values are given in Table1. In order to determine the effectiveness of the suggested technique, we contrasted the numerical and graphical solutions achieved by the Galerkin scheme with those obtained by the RK4 method and other conventional methods, i.e., MuHPM [41], MVIM [40], and HPM [38], provided in Tables 2–4 for Tt, It, and Vt. After the evaluation, it was revealed that the proposed approach produced more precise findings than the solutions obtained through conventional approaches. Furthermore, numerical computations are carried out using the Galerkin scheme and the RK4 scheme with the same and different time step sizes. The absolute errors between the outcomes were assessed using both approaches, as demonstrated in Tables 5 and 6. From comparison, it is clearly visible that the Galerkin scheme achieves more precise findings at a larger step size as compared to the RK4 technique. Figure 2 shows the results of both methods for t∈0,1 and t∈0,70 with ȷ=0.1 steps. The figures showed that the results overlapped over each other throughout the time period. In addition, the analysis of the numerical costs for different step sizes in terms of time shows that the Galerkin method is more time-consuming as in comparison to the RK4 method is presented in Table 7 and Figure 3(f). The polar and bar graphs illustrate the absolute errors between the solutions of the Galerkin and RK4 approaches depicted in Figure 4. The mesh graphs for the solutions of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] are visualised in 3. Overall, we may conclude that the numerical approach reported in this article can be relied on to produce reasonably flexible and accurate results when applied to comparable problems.Table 2
The findings of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] for Tt.
tiMuHPM [41]MVIM [40]HPM [38]RK4 methodPresent solutions0.00.100000000.100000000.1000000000.10000000000000.10000000000000.10.146358950.146359100.1463589000.14635633342150.14635849282220.20.208807720.208808100.2087991000.20880067887670.20880649648410.40.406239230.406240800.4056066000.40621367496460.40623478431240.60.764420340.764428700.7564485000.76435081453700.76440824443910.81.414038301.414094001.3642150001.41387024891631.41400906112661.02.591575702.591921002.3786790002.59119519036622.59150945893731.24.723924184.725783004.0065120004.72309837340384.72378033501621.48.578327308.587223006.5213030008.57658959800018.57802431662681.615.522792515.561670010.2735700015.51922813062515.5221702037551.827.961354328.119610015.7007900027.95421977630627.9601062983772.050.008070750.614770023.3374000049.99418511292750.0056374736263.0605.328400758.304100121.4166000605.19715873002605.309626004254.01387.07811-38782.630428.18300001387.03656989361387.07346477905.01484.41661-1.2822E+71179.6330001484.41297411281484.416218519410578.768500-7.535E+1831121.36000578.01860976666578.803221775333042.9620000-7.954E+656710438.00042.98826240709642.9794781819207051.1260000-8.86E+1594.49161E+0850.73161364454350.729213455592Table 3
The findings of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] for It.
tiMuHPM [41]MVIM [40]HPM [38]RK4 methodPresent solutions0.00.000000000.000000000.000000000.00000000000000.00000000000000.10.000002860.000002860.000002860.00000286453320.00000286484710.20.000006030.000006030.000006030.00000603187890.00000603254640.40.000013160.000013150.000013150.00001315648600.00001315797850.60.000021220.000021220.000021220.00002122067760.00002122315970.80.000030180.000030170.000029940.00003017281020.00003017646611.00.000040040.000040020.000039180.00004003141580.00004003645811.20.000050880.000050840.000048440.00005087001740.00005087669631.40.000062830.000062700.000056970.00006281620700.00006282482191.60.000076080.000075740.000063540.00007606842280.00007607934211.80.000090960.000090110.000066490.00009094185170.00009095554992.00.000107990.000105970.000063460.00010796867950.00010798581413.00.000305680.00021559-0.00018160.00030559579170.00030566741484.00.002109200.00040533-0.00148370.00210841864440.00210907783795.00.020119310.00074045-0.00550380.02011076952760.020117982587410.0828.5948000.01387534-0.2285212828.38422383941828.5668668426630.0775.2937001625.28280-64.142730775.29918583941775.2857747503970.0594.7040002.2292E+13-4621.4330594.00643606878594.00643962536Table 4
The findings of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] for Vt.
tiMuHPM [41]MVIM [40]HPM [38]RK4 methodPresent solutions0.00.100000000.100000000.100000000.10000000000000.10000000000000.10.078663150.078663180.078663150.07866381455020.07866326355950.20.061879800.061879910.061878250.06188084740160.06187998051660.40.038294840.038295960.038199470.03829613043710.03829505758290.60.023704500.023710290.022681580.02370570312840.02370470746760.80.014680330.014700420.009255110.01468131432290.01468049325491.00.009100820.00915723-0.01048470.00910157907680.00910094475311.20.005653260.00579375-0.04982590.00565382162980.00565335183931.40.003525410.00385124-0.12948780.00352581552480.00352547818281.60.002214750.00293882-0.28016880.00221503826520.00221480224191.80.001410460.00297863-0.54509450.00141065416330.00141049332042.00.000920390.00426490-0.98256560.00092052056450.00092041462753.00.000329410.13649810-9.17492000.00032936511840.00032939832824.00.001410905.40689000-43.5530800.00141038094990.00141082558615.00.01292825214.291400-143.651800.01292275834970.012927397503710.0668.7895902.095E+10-5456.7680668.57158158359668.7579298116730.0983.9820001.916E+42-150636600983.92573498469983.9243818601770.0769.4900001.602E+106-1.0827E+8768.96288598031768.95708309906Table 5
The absolute errors between the findings of the Galerkin and RK4 schemes forTt, It, and Vt with same step size and t∈0,70.
cGP2j=0.01−RK4j=0.01tiTtItVt0.00.00000000000000E-100.00000000000000E-140.00000000000000E-110.12.84775258929670E-103.46123001238724E-144.45784936742299E-110.27.67201135953854E-107.50830062186403E-147.01317059981221E-110.31.55009055591293E-091.21598593011158E-138.27457546925814E-110.42.78372330742016E-091.74456149260862E-138.67752675270950E-110.54.68632355143939E-092.34046319786763E-138.53053669369608E-110.67.57295037789874E-093.00841803183537E-138.04949266963728E-110.71.18961085426861E-083.75390284178614E-137.38314721859812E-110.81.83025996580710E-084.58311269117091E-136.63192608468810E-110.92.77127207848338E-085.50296584017097E-135.86177772987639E-111.28.98785703640215E-088.88763672032833E-133.77432252118837E-111.41.89002486550294E-071.17644883093106E-122.69146562675848E-111.63.87287119707480E-071.52607589651126E-121.85752377021564E-111.87.73954734967219E-071.95361702347036E-121.23150862337690E-112.01.50285070077416E-062.48521210934657E-127.65154150976499E-123.01.44064168807745E-051.06163661303844E-116.34902492992889E-124.04.45008413407777E-069.61197484453114E-116.36393082328879E-115.03.44141881214455E-071.02496256620954E-096.58341856724087E-10102.90578410613307E-052.45547904569321E-052.57475288663045E-05301.43229695481750E-062.60882575275900E-065.73396846448304E-07704.48719525536490E-071.54975691657455E-061.31251681523281E-06Table 6
The absolute errors forTt, It, and Vt between the Galerkin and RK4 techniques for t∈0,70 with different step sizes.
cGP2j=0.01−RK4j=0.0001tiTtItVt0.00.00000000000000E-110.00000000000000E-150.00000000000000E-120.15.86432846727547E-117.12909557572966E-158.69847249784783E-120.21.57959229030169E-101.54760669814828E-141.36827840724330E-110.33.19075987853523E-102.50826784472732E-141.61410537713458E-110.45.72855207714440E-103.60121817349207E-141.69235597757833E-110.59.64069601927520E-104.83457641053866E-141.66324766459613E-110.61.55727908346393E-096.21799635977353E-141.56892346336868E-110.72.44506304092340E-097.76247081478616E-141.43842576738606E-110.83.75948450148655E-099.48020382628059E-141.29135296206284E-110.95.68795943678424E-091.13844887293310E-131.14057773642484E-111.08.49506598399330E-091.34896515615809E-139.94224494399454E-121.21.83776105444622E-081.83649035106863E-137.31735158576718E-121.43.84712617318428E-082.42403387401731E-135.19390173231438E-121.67.82565194867857E-083.12689110715247E-133.55596108114753E-121.81.54546054176308E-073.96344591803606E-132.32396979264748E-122.02.94468271988535E-074.95918884842010E-131.40495156741166E-123.02.23191591430805E-061.76651054895061E-121.14225539834048E-124.05.47712943443912E-071.65905887118678E-111.07716444271155E-115.04.52982931165025E-081.77432183529058E-101.13838096935104E-10104.84750069063011E-064.06695653509814E-064.55263102594472E-06302.30107175980265E-074.45262458015350E-071.09476786747109E-07707.22985902257278E-082.61271452473011E-072.25629946726258E-07Figure 2
The geometrical comparison between the findings of the Galerkin and RK4 schemes forTt, It, and Vt for t∈0,70 and t∈0,1 with j=0.1.
(a)(b)(c)(d)(e)(f)Table 7
The numerical cost in terms of the CPU time (in seconds) for the cGP(2) scheme and RK4 scheme fort∈0,70.
CPU time (in seconds)1/jcGP(2) methodRK4 method700.8513710.0584631000.9589540.0540762001.6077080.0677535002.6876630.08498110004.4738710.12477515005.2455370.19959320006.0210620.234840500010.2441010.5153861000017.5105890.8591492000031.6600671.945279Figure 3
The mesh graphs of the Galerkin method, RK4 method, MuHPM [41], MVIM [40], and HPM [38] for the model and the computation cost in terms of CPU time (f).
(a)(b)(c)(d)(e)(f)Figure 4
The polar and bar graphs of errors between the findings of the Galerkin and RK4 schemes and bar graphs for the Galerkin and RK4 schemes forTt, It, and Vt.
(a)(b)(c)(d)(e)(f)


## 4. Mathematical Formulation of the HIV Infection Model

The study of HIV/AIDS infection requires an understanding of T-cell population dynamics because T-cells are the major target of HIV particles. Since these cells are generated in the bone marrow, premature T-cells are directed to the thymus for further differentiation and maturation before becoming immune-competent T-cells. However, HIV may be able to infect these cells in both bone marrow and the thymus and destroy the supply of new cells [54]. But majority of the HIV models in the literature assumed with the constant rate of generation of new T-cells from the thymus. Herein, we modified the HIV model assumed by Leal et al. [41]. We consider the innovative model with the assumption that the production of new cells from the thymus is a decreasing function based on viral load (“sV=0.5Is+sI/I+V” used in [43, 44, 55]) rather than the constant source term. Further, include mass action term “kVT,” where k is the infection rate. Infection occurs when a virus infects T-cells, reducing them at the rate of “−kVT” and producing virus at the rate of “kVT.” It is observed that the mentioned model gives more accurate and realistic-based findings for T-cells, I-cells, and HIV. The new model is follows:
(29)dTdt=0.5Is+sII+V−μTT+rT1−T+ITmax−kVT,dIdt=kVT−μII,dVdt=NμII−μVV+kVT,where Tt, It, and Vt represent the density of T-cells, I-cells, and virus particles in the blood, respectively. The diagrammatic illustration of the HIV model and the structure of the HIV virus are illustrated in Figure 1. Following confirmation of the aforementioned scheme and validation of the MATLAB code, we applied the aforesaid scheme to a novel established model and visualised the dynamical behaviour of the model with the interaction of major clinical parameters. The initial conditions of state variables, various parameters, and constants with their explanations are provided in Table 1 and I=1mm-3.


### 4.1. Analysis and Discussions

In this section of the research, several numerical simulations are carried to illustrate the detailed behaviour of the proposed model. Figures5(a)–5(c) depict that by increasing the quantity of free virus produced by reduction of T-cells, the concentration of T-cells increases at top level initially. But after ten days approximately, its wavelength starts decreasing with time and touches the lowest level. Due to killing of viral particles and cell response, the density of T-cells increases again after initial decline. However, immune cells are directly impacted by viral particles; their number steadily diminishes over time, indicating that T-cells are targeted by viruses. The increasing value of “N” gives rapid rise to the depletion of T-cells at huge level and is more predictable sign for AIDS. The profile of I-cells shows oscillatory behaviour, and the wavelength is increasing by increasing the generation rate of free virus “N.” The viruses attack the immunity cells and convert these T-cells into I-cells. Due to such virus replication and infection of T-cells, the amount of I-cells increases. Meanwhile, the immune system has the capability to identify viruses and destroy them to control the infection. Therefore, the rate of I-cells again decreases after seventy days approximately. The influence of growth rate (“r”) of the T-cell population on the dynamical behaviour of sate variables Tt,ItVt is visualised in Figures 5(d)–5(f). From the graphical view, it could be observed that state variables show considerably the same dynamics by increasing the growth rate and death rate but speeds up the decaying oscillation. The impact of increasing the death rate μT of the T-cell population on the dynamical behaviour of sate variables is depicted in Figures 6(a)–6(c). If the mortality rate of T-cells increases, there is no significant influence on the density of T-cells in the earliest stages. Since the virus attacks T-cells on a continuous basis, a gradual but significant decline occurs after seventy days. Similarly, the time period for the concentration of I-cells and free HIV virus particles is affected slightly. Figures 6(d)–6(f) demonstrate that if the mortality rate of the I-cell population grows, the rate of T-cells reaches a maximum in the early days but begins to drop in concentration with the passage of time, and a continuous pattern of crest and trough with a small wavelength in density of T-cells could be noticed. The population of I-cells diminishes, although time has a vital effect, as seen in the graphs. Since HIV particles are latently present, they affect the number of I-cells, and their patterns display the same behaviour as I-cells. The infection rate oscillates over a seventy-day time-frame, and the death rate of I-cells has a significant effect on the growth of viral particles. By reducing the mortality rate of virus particles, the concentration of T-cells decreases with time, whereas speeding the depletion with longer wavelengths, as illustrated in Figure 7(a). The HIV particles target T-cells; their concentration rises if the mortality rate of HIV particles is lowered. Also, at the end of ten days, the initial and final densities of T-cells are almost identical, whereas the density of I-cells and viral particles grows and the time for growth in I-cell and viral particles decreases, as shown in Figures 7(b) and 7(c).Figure 5
The impact of significant parameters on the HIV infection model for T-cell generation with varying source term (dependent on viral load).
(a)(b)(c)(d)(e)(f)Figure 6
The impact of significant parameters on the HIV infection model for T-cell generation with varying source term (dependent on viral load).
(a)(b)(c)(d)(e)(f)Figure 7
The impact of significant parameters on the HIV infection model for T-cell generation with varying source term (dependent on viral load) and phase diagram.
(a)(b)(c)(d)(e)(f)The phase and chaotic behaviour of the specified model are shown in Figures7(d)–8(b). Many scientific and engineering applications are based on a system’s chaotic nature [34]. It is well acknowledged that there is a substantial tendency toward understanding and representing chaotic system behaviours. The chaotic diagrams demonstrate the feasibility and application of the suggested numerical approach, which may be generalized to novel chaotic systems.Figure 8
The phase diagram of HIV infection model including T-cell generation with varying source term (dependent on viral load).
(a)(b)


## 4.1. Analysis and Discussions

In this section of the research, several numerical simulations are carried to illustrate the detailed behaviour of the proposed model. Figures5(a)–5(c) depict that by increasing the quantity of free virus produced by reduction of T-cells, the concentration of T-cells increases at top level initially. But after ten days approximately, its wavelength starts decreasing with time and touches the lowest level. Due to killing of viral particles and cell response, the density of T-cells increases again after initial decline. However, immune cells are directly impacted by viral particles; their number steadily diminishes over time, indicating that T-cells are targeted by viruses. The increasing value of “N” gives rapid rise to the depletion of T-cells at huge level and is more predictable sign for AIDS. The profile of I-cells shows oscillatory behaviour, and the wavelength is increasing by increasing the generation rate of free virus “N.” The viruses attack the immunity cells and convert these T-cells into I-cells. Due to such virus replication and infection of T-cells, the amount of I-cells increases. Meanwhile, the immune system has the capability to identify viruses and destroy them to control the infection. Therefore, the rate of I-cells again decreases after seventy days approximately. The influence of growth rate (“r”) of the T-cell population on the dynamical behaviour of sate variables Tt,ItVt is visualised in Figures 5(d)–5(f). From the graphical view, it could be observed that state variables show considerably the same dynamics by increasing the growth rate and death rate but speeds up the decaying oscillation. The impact of increasing the death rate μT of the T-cell population on the dynamical behaviour of sate variables is depicted in Figures 6(a)–6(c). If the mortality rate of T-cells increases, there is no significant influence on the density of T-cells in the earliest stages. Since the virus attacks T-cells on a continuous basis, a gradual but significant decline occurs after seventy days. Similarly, the time period for the concentration of I-cells and free HIV virus particles is affected slightly. Figures 6(d)–6(f) demonstrate that if the mortality rate of the I-cell population grows, the rate of T-cells reaches a maximum in the early days but begins to drop in concentration with the passage of time, and a continuous pattern of crest and trough with a small wavelength in density of T-cells could be noticed. The population of I-cells diminishes, although time has a vital effect, as seen in the graphs. Since HIV particles are latently present, they affect the number of I-cells, and their patterns display the same behaviour as I-cells. The infection rate oscillates over a seventy-day time-frame, and the death rate of I-cells has a significant effect on the growth of viral particles. By reducing the mortality rate of virus particles, the concentration of T-cells decreases with time, whereas speeding the depletion with longer wavelengths, as illustrated in Figure 7(a). The HIV particles target T-cells; their concentration rises if the mortality rate of HIV particles is lowered. Also, at the end of ten days, the initial and final densities of T-cells are almost identical, whereas the density of I-cells and viral particles grows and the time for growth in I-cell and viral particles decreases, as shown in Figures 7(b) and 7(c).Figure 5
The impact of significant parameters on the HIV infection model for T-cell generation with varying source term (dependent on viral load).
(a)(b)(c)(d)(e)(f)Figure 6
The impact of significant parameters on the HIV infection model for T-cell generation with varying source term (dependent on viral load).
(a)(b)(c)(d)(e)(f)Figure 7
The impact of significant parameters on the HIV infection model for T-cell generation with varying source term (dependent on viral load) and phase diagram.
(a)(b)(c)(d)(e)(f)The phase and chaotic behaviour of the specified model are shown in Figures7(d)–8(b). Many scientific and engineering applications are based on a system’s chaotic nature [34]. It is well acknowledged that there is a substantial tendency toward understanding and representing chaotic system behaviours. The chaotic diagrams demonstrate the feasibility and application of the suggested numerical approach, which may be generalized to novel chaotic systems.Figure 8
The phase diagram of HIV infection model including T-cell generation with varying source term (dependent on viral load).
(a)(b)


## 5. Concluding Remarks

In this paper, the continuous Galerkin-Petrov time discretization scheme is implemented for the model for HIV infection, which is comprised of three nonlinear ordinary differential equations. Afterwards, implement the RK4 method for solving the model and compare the solutions of both techniques. Moreover, a comprehensive analysis of the outcomes of the Galerkin method and RK4 method with findings of other schemes available in the literature is presented in detail. The assessment of the absolute values of errors between the Galerkin scheme and RK4 scheme solutions shows that the solutions provided are more accurate in comparison to the solutions obtained through conventional schemes. The graphical and tabular outcomes demonstrated that the proposed scheme is very accurate and achieves highly accurate results at a larger step size in comparison to the existing techniques employed for the mentioned model. The Galerkin scheme is an adaptive scheme that achieves the same accuracy at a larger step size at a lower cost. In comparison towards other approaches, our technique yields more accurate results than other methods that have addressed the same problem. The findings demonstrate that the approach employed may be utilised for numerous types of nonlinear systems of differential equations. The performance of the proposed approach demonstrates that these methods are logically effective and reliable for solving nonlinear problems in complex dynamical systems. On the other hand, we examined the dynamical behaviour of the model with a variable supply rate (depending on viral load) of new T-cells. The influence of significant clinical parameters on the dynamics of T-cells, I-cells, and free HIV virus particles is described and observed by varying their values. The main findings are as follows:(1)
Increasing the density virus produced by per I-cellsN increases the concentration of T-cells, I-cells, and virus particles(2)
The decaying oscillatory behaviour is observed in the density of state variables by rising the value ofr(3)
By enhancing the mortality of I-cells, the density of T-cells and virus rises while the population of I-cells diminishes(4)
The increase in population of all state variables is observed by decreasing the mortality rate of HIVIn addition, the phase and chaotic behaviour of the aforementioned model are presented. Many scientific and technical applications are predicated on the chaotic character of a system. It is widely accepted that there is a heavy tendency toward comprehending and portraying chaotic system behaviours. The chaotic diagrams show the practicality and accessibility of the suggested numerical technique, which may be expanded to novel chaotic systems.In the future, we plan to elaborate on the research by incorporating vaccination rates and investigating their impact on the model’s dynamical behaviour.



---

*Source: 3599827-2022-11-09.xml*