On the comparative performance of fourth order Runge-Kutta and the Galerkin-Petrov time discretization methods for solving nonlinear ordinary differential equations with application to some mathematical models in epidemiology

1.   Introduction

Epidemiology is the foundation of public health and is defined as the study of the distribution and determinants of diseases or disorders within groups of people, and the development of knowledge on how to prevent and control them. Epidemiological research helps us understand not only who has a disorder or disease but why and how it was brought to this individual or region. Nowadays, epidemiologists use the insights gathered in their research to determine how illness within a population affects our society and systems on a larger scale, and, in turn, they provide recommendations for interventions. The hepatitis B virus (HBV) became widespread, and epidemiologists around the world worked to control the spread. The epidemiology of HBV infection is diverse, with population prevalence, age and mode of acquisition and likelihood of progression to chronic infection mutually interdependent. Hepatitis B is a life-threatening liver disease caused by the HBV ^[1,2], which has a circular genome made of partly double-stranded DNA and is difficult to eradicate after infection due to the development of cccDNA ^[1,3]. HBV may induce both acute sickness and chronic liver infection, with sustained blood levels of HBV surface antigen (HBsAg), IgG anti-core antigen (anti-HBc), and HBV DNA ^[3]. Chronic infection may progress to cirrhosis or liver cancer ^[4]. According to the World Health Organization, around 2 billion individuals worldwide have been infected with the virus, with approximately 350 million having HBV. Each year, HBV causes approximately 600, 000 fatalities. Hepatitis B is prevalent in China, with an estimated 93 million people afflicted ^[5]. Due to the high risk of infection and the high number of fatalities associated with HBV, it is critical to improve our knowledge of the dynamics of HBV illness. HBV is spread by contact with the blood or other bodily fluids of the infected individuals ^[2]. Generally, the average age at which people get affected is influenced by the population's infection incidence ^[5]. Transmission from mother to baby occurs throughout pregnancy, delivery and the postpartum period. Vaccination, including passive immunoprophylaxis and postnatal proactive vaccination, is effective in protecting infants against contracting HBV during labor and delivery and in the postnatal period, but it cannot halt HBV prenatal infection ^{[6,7,8,9,10,11]}. Mathematical models may assist us in gaining insight into disease transmission, evaluating the success of different preventative methods, and ultimately bringing the illness under control. Mathematical modeling of a dynamical system is an interesting study area that has grabbed the interest of a majority of researchers. Dynamical systems have a wide range of applications, including population growth models, biomedical research, biological systems, and engineering. Aniji et al. ^[12] established a mathematical model of hepatitis B virus infection for antiviral treatment that depicts the probable interaction between HBV and target liver cells. Zada et al. ^[13] evaluated and analyzed the hepatitis B epidemic model with optimal control. Zhao et al. ^[14] designed a mathematical model of hepatitis B viral transmission and addressed its potential use in China's immunization plan. Means et al. ^[15] explored the Hepatitis B Virus's geographical impacts and potential function. Khatun et al. ^[16] explored viral dynamics, which aids in the understanding of the dynamics of HIV, hepatitis B virus (HCB), hepatitis C virus (HCV), and many other viruses. Additionally, they explored the hepatitis B virus model with Cytotoxic T Lymphocyte (CTL) immunological responses. Zhang et al. ^[17] developed a mathematical model to explain how hepatitis B spreads. The stability of equilibria and illness persistence are investigated. The study demonstrated that the model's dynamics are entirely determined by the basic reproductive number. Khan et al. ^[18] reported a detailed analysis of the Hepatitis B model using the Caputo-Fabrizio fractional derivative. The required findings are obtained by establishing the criteria for the existence and uniqueness of the solution to the considered model using fixed point theory. The continuous Galerkin-Petrov ($\mathrm{cGP}$) discretization scheme is a very popular and powerful approach to handling complex dynamical systems. Numerous researchers addressed so many interesting properties and some feasible practical applications of the proposed scheme's implementation. Attaullah and Sohaib ^[19] employed a continuous Galerkin-Petrov ($\mathrm{cGP}$) discretization scheme with a mathematical model of Human Immunodeficiency Virus (HIV) infection, which contained a system of nonlinear differential equations. They validated the results of the proposed scheme in comparison to other traditional schemes. Hussain et al. ^[20] utilized the Galerkin method with the Chen model and performed a comparison between the findings of the Galerkin technique and the fourth-order Runge-Kutta (RK4) method for the presented problem. They concluded that the Galerkin scheme achieved more accurate results at larger time step sizes than the RK4 method. Hussain ^[21] implemented a novel class of higher-order Galerkin temporal discretization techniques for non-stationary incompressible flow problems. Attaullah et al. ^[22] implemented the Galerkin time discretization scheme for the transmission dynamics of HIV infection with nonlinear supply rate and discussed the impact of different physical parameters on the population dynamics of healthy/infected cells and virus particles. Schieweck ^[23] described A-stable discontinuous Galerkin-Petrov time discretization of higher-order. For the heat equation, higher-order Galerkin temporal discretizations and fast multigrid solvers were used by Hussain et al. ^[24]. For the non-stationary Stokes equations, Hussain et al. ^[25] demonstrated an accurate and efficient higher-order Galerkin time-stepping approach. For nonlinear systems of ordinary differential equations, higher-order variational time discretizations were utilized by Matthies and Schieweck ^[25]. Attaullah et al. ^[26] numerically simulated the HIV model using a Galerkin approach. To validate the numerical results, a detailed evaluation of the outcomes achieved by the Galerkin scheme is contrasted with the findings of the RK4 technique. Various simulations were carried out to evaluate the findings of the Galerkin and RK4 approach, demonstrating the credibility of the $\mathrm{cGP}$(2) technique. The comparison analysis and graphical plots illustrate that the Galerkin scheme is more appropriate than the RK4 method for achieving high accuracy at larger step sizes. For the approximate solution of the HIV model with a variable source term and full logistic proliferation, Attaullah et al. ^[27] implemented the Galerkin-technique and compared the results with the findings of the RK4 method. Khan et al. ^[28] described a mathematical model for the transmission dynamics of acute and chronic hepatitis B and developed an optimal control method for hepatitis B spread in a community. Numerous models of HBV transmission dynamics have been developed to examine the effect of preventative and control methods such as vaccination, antiviral therapy and the incidence rate ^{[29,30,31,32]}. These models have provided valuable insight into the effectiveness of different control strategies. However, the majority of these models have been studied primarily analytically. The mathematical model proposed in the paper is the extension of the model presented by Khan et al. ^[28] by including a non-linear saturated incidence rate ^[32] and investigated numerically. We carry out a comprehensive investigation of the various influence of several clinical parameters. The main contributions are as follows:

$1.$Implemented Galerkin scheme and the RK4 scheme correctly to the HBV model and presented the solutions numerically and graphically.

$2.$The significant effects of the various physical parameters are visualized accurately.

$3.$To authenticate the validity and reliability of the suggested scheme, employed the proposed schemes to other mathematical models (Chen System and the HIV infection model) existing in the literature.

$4.$The comparative study is performed to validate the accuracy (adopting various step sizes) and computational cost (in terms of time) of Galerkin and RK4-scheme.

$5.$Finally, invoked the methods to the problem having an exact solution and confirm the accuracy and time consuming with different step sizes.

2.   The mathematical model of Hepatitis B Virus

The mathematical model addressed here is as follows:

In the this model the population $\tilde{T}(t)$ is divided into four compartments: susceptible individuals $\tilde{S}(t),$ infected individuals $\tilde{I}_{1}(t),$ infected individuals with chronic hepatitis B (CHB), $\tilde{I}_{2}(t)$ and recovered individuals $\tilde{R}(t).$ The initial conditions and the detailed information of all parameters included in the model are presented in Table 1.

3.   The numerical schemes

3.1.   The continuous Galerkin-Petrov method

The Galerkin technique has been widely employed in recent years to handle a wide range of nonlinear problems in science and engineering, such as ^{[23,24,25,33,34,35,37,43,44]}. In this paper, we are using this approach for the models discussed in ^{[19,20,27,28,38]}. The system of ODEs for the considered model can be written as follows:

Find $\mathbf{\tilde{\mathbf u}}:[0, t_{max}]\to \mathbf{V} = {\mathbb R}^d$ such that

where $d_{t}$ shows the derivative of ${\mathbf {\tilde{u}}}(t)$ w.r.t. time, $I$ = [0, T] is the time interval, ${\mathbf {\tilde{u}}}(t) = (\tilde{u}_1(t), \; \tilde{u}_2(t), \; \tilde{u}_3(t))\; \in\; V \; \Rightarrow {\mathbf {\tilde{u}}}(0) = (\tilde{u}_1(0), \tilde{u}_2(0), \tilde{u}_3(0))\in V$ are the initial values of ${\mathbf {\tilde{u}}}(t)$ at time $t = 0$. We also assume that $(\tilde{u}_1(t), \; \tilde{u}_2(t), \; \tilde{u}_3(t)) = (T(t), \; I(t), \; V(t)),$ which implies that $(\tilde{u}_1(0), \; \tilde{u}_2(0), \; \tilde{u}_3(0)) = (T(0), I(0), V(0))$. The function $F = (f_1, f_2, f_3)$ is nonlinear and is defined as F : $I$ $\times$ $K$ $\rightarrow$ $K.$

The weak formulation of Problem (3.1) is: Find $\tilde{u}$ $\in$ X such that ${\mathbf {\tilde{u}}}(0) = {\mathbf {\tilde{u}}}_0$ and

where $\langle\cdot, \cdot\rangle$ represents the usual inner product in $L^2(I, K),$ and $X,$ $Y$ represent the solution space and test space, respectively. To explain the time discretization of variational type Problem(3.1).

Characterize the function $t$ $\rightarrow {\mathbf {\tilde{u}}}(t),$ and we define the space $C(I, K) = C^0(I, K)$ as the space of continuous functions $\tilde{u}$ : $I$ $\rightarrow$ $K$ equipped with the norm

We will use the space $L^2(I, K)$ as the space of discontinuous functions, which is given by

In time discretization, we split the intervals $I$ into $N$ sub-intervals $I_n = [t_{n-1}, \; t_n]$, where $n = 1, \cdots, N$ and $0 = t_0 < t_1 < t_2 < \cdots < t_{N-1} < t_N = T$. The symbol $\tau$ will indicate the time discretization parameter, as well as the maximum time step size $\tau = \max \limits_{1\leq n\leq N}\tau_n$, where $\; \tau_n = t_n-t_{n-1},$ the length of nth time interval $I_n$. The following set of time intervals ${ M_{\tau}} = \{\bar {I}_1, \cdots, \bar{I}_N\}$ will be called the time-mesh. We find out the solution ${\mathbf {\tilde{u}}}$ : $I$ $\rightarrow$ $K$ on each time interval $I_n$ by a function ${\mathbf {\tilde{u}}}_\tau$ : $I$ $\rightarrow$ $K$ which is a piecewise polynomial of some order $l$ w.r.t time. The time-discrete solution space for ${\mathbf {\tilde{u}}}_\tau$ is $X_\tau^l \subset X$ and is defined by

where

The discrete test space for ${\mathbf {\tilde{u}}}_\tau$ is $Y_\tau^l\subset Y$ and is defined by

which is composed of $l-1$ piecewise polynomials and is discontinuous at the time step end nodes. We multiply Eq (3.1) by the test function $v_\tau \in Y_\tau^k$ and integrate it over the whole interval $I$. We get the discrete-time problem: Find ${\mathbf {\tilde{u}}}_\tau \in X_\tau^k$ such that ${\mathbf {\tilde{u}}}_\tau(0) = 0$ and

where $\langle\cdot, \cdot\rangle$ represents the usual inner product in $L^2(I, K).$ This discretization is known as the exact continuous Galerkin-Petrov method or briefly the "exact $\mathrm{cGP}$(l)-method'' of order k. The Galerkin-Petrov name is due to the fact that the solution space $X_\tau^l$ is different from the test space $Y_\tau^l$. The term "exact" denotes that the time integral on the right-hand side of the Eq (3.5) is determined exactly. Because the discrete test space $Y_\tau^l$ is discontinuous, Eq (3.5) can be solved in a time marching procedure in which local problems on the time interval are handled successively. Therefore, we choose the test function $v_\tau(t) = v\psi(t)$ with arbitrary time-independent $v\in K$ and a scalar function $\psi : I \rightarrow \mathbb{R}$ which is zero on $I|_{I_n}$ and a polynomial of the order less than or equal to $l-1$ on the time interval $I_n = [t_{n-1}, \; t_n]$. Then, we get from (3.5) the $I_n$-problem of the exact continuous Galerkin-Petrov method of order l: Find ${\mathbf {\tilde{u}}}_\tau |_{I_n}\in \mathbb{P}_{l}(I_n, K)$ such that

with the initial condition ${\mathbf {\tilde{u}}}_\tau |_{I_n}(t_{n-1}) = {\mathbf {\tilde{u}}}_\tau |_{I_{n-1}}(t_{n-1})$ for $n\geq 2$ and ${\mathbf {\tilde{u}}}_\tau |_{I_n}(t_{n-1}) = \tilde{u}_0$ for $n = 1$.

In the usual case of a nonlinear function $F\langle\cdot, \cdot\rangle$, (where $\langle\cdot, \cdot\rangle$ denotes the usual inner product in $L^2$ norm) we need to calculate the integrals numerically on the right hand side of the Eq (3.6). The $(l+1)$-point Gauss-Lobatto formula is exact if the function is to be integrated as 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 (3.6) will give the exact value. Then, the "$I_n$-problem of the numerically integrated $\mathrm{cGP}(l)$ method'' is: Find ${\mathbf {\tilde{u}}}_\tau |_{I_n}\in\mathbb{P}_{l}(I_n, K)$ such that ${\mathbf {\tilde{u}}}_\tau(t_{n-1}) = {\mathbf {\tilde{u}}}_{n-1}$,

where $\hat{w}_s$ are the weights and $\hat{t}\in[-1, 1], s = 0, 1, 2, 3, \ldots, l$ represent the nodes on the reference interval. To find ${\mathbf {\tilde{u}}}_\tau$ on each time interval $I_n$, we represent it by a polynomial ansatz

where the coefficients $U^{s}_{n}$ are the components of $K$ and the real-valued functions $\phi_{n, s} \in \mathbb{P}_{l}(I_n)$ are the Lagrange basis functions with respect to $l+1$ suitable nodal points $t_{n, s}\in I_n$ satisfying conditions mentioned below:

where $\delta_{r, s}$ is the Kronecker delta that is defined as:

Like in ^[36], the $t_{n, s}$ have been defined as the quadrature points of $(l+1)$-point Gauss-Lobatto formula on the time interval $I_{n}$. For the selection of initial conditions we can set $t_{n, 0} = t_{n-1}$ which implies the initial condition for (3.6)

We define the basis functions $\phi_{n, s}\in\mathbb{P}_{l}(I_n)$ via the affine reference transformation $\bar{T} : \hat{I}\rightarrow I_n$ where $\hat{I}$ = [-1, 1] and

Let $\hat{\phi_s}\in\mathbb{P}_{l}(\hat{I})$, s = 0, 1, $\ldots, l,$ demonstrate the basis functions that meet the requirements

where $\hat{t_0} = -1$ and $\hat{t_r}$, $r = 1, 2, \ldots, l$, are the quadrature points for the reference interval $\hat{I}$. Then, we define the basis functions on the original time interval $I_n$ by the mapping

Likewise, the test basis functions $\psi_{n, r}$ are defined by the suitable reference basis functions $\hat{\psi}\in\mathbb{P}_{l-1}(\hat{I})$, i.e.,

From the representation (3.8), we get for $d_{t}\tilde{u}_{\tau}$

By putting (3.13) in (3.6), we get

The integral is now transformed into the reference interval $hatI$ and computed using the $(l+1)$-point Gauss-Lobatto quadrature formula, which leads for each test basis function $\psi\in\mathbb{P}_{l-1}$ and for all $v\in K$,

where $\hat{w}_\mu$ are the weights and $\hat{t}_\mu\in [-1, 1]$ are the points of integration with $\hat{t}_0 = -1$ and $\hat{t}_l = 1$. If we choose the test functions $\psi_{n, i}\in\mathbb{P}_{l-1}(I_n)$ such that

Now find the coefficients that are unknown $U_n^s\in K$ for $s = 1, \ldots, l,$

where $U^s_n = U^s_{n-1}$ for $n > 1,$ $U^0_1 = {\mathbf {\tilde{u}}}_0$ for $n = 1$ and

We will discuss the $\mathrm{cGP}$(k) method for the cases $l = 1\; and\; l = 2$ in the following.

3.1.1.   The cGP(1) method

We used the two point Gauss-Lobatto formula with $t_{n, 0} = t_{n-1}, \; t_{n, 1} = t_n$ and weights $\hat{w}_0 = \hat{w}_1 = 1,$ which gives the well-known Trapezoidal rule. We obtain $\alpha_{1, 0} = -1, \alpha_{1, 1} = 1$ and $\beta_1 = 1$. For the single coefficient $U^1_{n} = {\mathbf {\tilde{u}}_\tau}(t_n)\in K$, the problem leads to the following block equation:

3.1.2.   The cGP(2) method

The three-point Gauss-Lobatto formula (Simpson rule) is used to define the quadratic basis functions with weights $\hat{w}_0 = \hat{w}_2 = 1/3, \hat{w}_1 = 4/3$ and $\hat{t}_0 = -1$, $\hat{t}_1 = 0$, $\hat{t}_2 = 1$. Then, we get

Thus, the system to be solved for $\mathbf{U}_n^1, \mathbf{U}_n^2 \in \mathbf{K}$ from the known $\mathbf{U}_n^0 = \mathbf{U}_{n-1}^2$ becomes:

where $\mathbf{U}_n^0$ indicate the initial conditions at the current time interval. For the cGP(k)-method, the discrete solution space consists of continuous piecewise polynomial functions in time of degree $k > 1$ and the discrete test space of discontinuous Galerkin (dG) having polynomial functions of degree $k - 1$ (see ^{[33,34,35,36]} for more details). In the dG($k$)-method, both the solution and test space are constructed by means of discontinuous polynomial functions of degree $k.$ With respect to the computational costs, which mainly depend on the size of the resulting block system that has to be solved for each time interval, the cGP(k)-method is comparable to the dG($k-1$)-method. However, concerning the discretization error in time, the accuracy of the cGP(k)-method is one order higher than that of the dG($k -1$)-method. Furthermore, it is known that all cGP($k$)-methods are A-stable, while all dG(k)-methods are even strongly A-stable (or L-stable), i.e., the dG-methods have better damping properties with respect to high frequency error components (see ^{[33,34,35,36]} for more details).

3.2.   The Runge-Kutta method

This method is very well-known of order four established by Kutta ^[39] and extensively utilized for initial value problems (see ^[40] for more details). The Runge-Kutta method of order four is used to solve numerically the first order initial value problems. Let

be the initial value problem with the initial condition $y(a) = \alpha$, and let $N > 0$ be an integer and set $h = \frac{b-a}{N}$ is the step size. Partition the whole interval into $N$ sub-intervals with mesh points ${t_i = a+ih}$, for $i = 0, 1, 2, \cdots, N-1$. Then, the Runge-Kutta method of order four is described as

where

The Runge Kutta method of order four (RK-4) agrees with the Taylor series method up to terms of $\mathcal{O}(h^4)$. This method can be extended to solve a system of $n$ first-order differential equations. The generalization of the method is as follows.

Let

be the nth-order system of first-order initial value problems with the initial conditions

Use the notation $y_i^j,$ for each $i = 0, 1, 2, \cdots, N$ and $j = 1, 2, \cdots, n,$ to denote an approximation to $y^j(t_i).$ That is, $y_i^j$ approximates the $jth$ solution $y^(t)$ of (3.24) at the $ith$ mesh points $t_i.$ For the initial condition, set

Suppose that the values $y_i^1, \; y_i^2, \cdots, y_i^n$ have been computed. We obtain $y_{i+1}^1, \; y_{i+1}^2, \cdots, y_{i+1}^n$ by first calculating

for each $j = 1, 2, \cdots, n;$ and then

for each $j = 1, 2, \cdots, n.$ The values $k_1^1, k_1^2, \cdots, k_1^n$, must be computed before any of the terms of the form $k_2^j$ can be determined.

3.3.   Numerical simulations and discussions

This section relates to the implementation of the $\mathrm{cGP}$(2)-scheme and RK4 scheme for the model presented and figure out the correct results. All the parameters, values are presented in Table 1. The findings and absolute errors achieved utilizing cGP(2) scheme (at $\tau$ = 0.1) and RK4 method (at $\tau$ = 0.1) were compared and presented in Tables 2–5 and in Figure 2. The findings achieved utilizing the presented approach are considerably comparable to those obtained using the RK4-method with differing computational time. From comparison one could see that the results are not precise but similar up to nine digits approximately. The results demonstrated that the proposed methodologies yielded reliable and precise solutions and that both results overlapped significantly during the stipulated time period. We also solved the model using cGP(4) scheme (at $\tau$ = 0.2) and RK4 (at $\tau$ = 0.1) and found the absolute errors between the findings. Comparing the absolute values of the errors presented in Tables 6–9, we see that the higher order method cGP(4) achieves with a quite large time step size higher accuracy than the lower order method cGP(2) with smaller time step displayed in Tables 2–5. Moreover, the analysis of the numerical costs for different step sizes in terms of time show that the $\mathrm{cGP}$(2)-method is more time consuming as in comparison to RK4-method as presented in Table 10 and Figure 2e. Furthermore, Figure 1 depicts the influence of varying the parameters $\xi_1$ and $\psi_2$ on the dynamical behavior of the state variables. By boosting the value of $\xi_1$, we observed that the population of susceptible individuals grows, and the population of infected individuals with acute HBV, infected individuals with chronic HBV and recovered individuals declines. By increasing $\psi_2$, the concentration of susceptible and recovered people increases, whereas the population of acutely infected and chronically infected individuals declines. Also, graphical visualizations of error solutions and errors phase portraits using $\mathrm{cGP}$(2)-scheme (at $\tau$ = 0.0005 where $\tau$ is the step size) and RK4-scheme (at $\tau$ = 0.0005) are shown in Figure 3 for $t \in [0, 10].$ The findings demonstrate that the cGP(2) scheme outperforms the RK4 scheme in terms of time step accuracy while incurring higher computational costs.

3.4.   The Basic information of the computer used for simulations

● Windows: 8.1 single language

● Processor: Intel(R)Core(TM)i5-5200CPU@2.20GHz 2.20GHz

● Installed memory(RAM): 4.00GB

● System type: 64 bit Operating System, x64-based processor

4.   Example 1: The Chen system

The Chen system was initially established in 1999 by Chen and Ueta ^[41] and described as follows:

where $X(t)$, $Y(t)$ and $Z(t)$ are the dependent variables, and $\alpha,$ $\beta$ and $\gamma$ are positive constant parameters. According to bifurcation analysis assuming the parameters $\alpha = 35$ and $\gamma = 28,$ system of Equations 4.1–4.3 shows non-chaotic and chaotic behaviour for $\beta = 12$ and $\beta = 3,$ respectively ^[20].

For the Chen system solutions, Hussain et al. ^[20] employed and investigated the Galerkin-method for the model of both non-chaotic behavior and chaotic behavior. They numerically compared the Galerkin-method solutions to the classic $\mathrm{RK4}$-method solutions and determined that the Galerkin solutions with larger time steps achieved comparable accuracy to the $\mathrm{RK4}$ solutions with significantly smaller time steps. This opposed the claim that the findings of the Galerkin-method and $\mathrm{RK4}$-method are same.

4.1.   Chen system with non chaotic behavior

Herein, we consider the Chen system with non-chaotic behavior and implement the proposed methods for same and different step sizes to strengthen and confirm the claim that (1) the results of Galerkin and $\mathrm{RK4}$ are not equal, as could be seen from Tables 11–16. The computations are conducted with several time step sizes, and the accuracy is assessed at various time intervals. It could be observed that the results are not similar at $\tau = 0.01$ while identical up to eight digit, when $\tau = 0.0001.$ (2) Tables 17 and Figure 4(g) show the elapsed time of both schemes, and it may be worthy to note that the $\mathrm{RK4}$-scheme required less time as compared to the $\mathrm{cGP}$(2)-scheme. In addition, Figure 4 shows the differences between the outcomes obtained through the Galerkin and $\mathrm{RK4}$ methods for different values of $\tau$ in the non-chaotic case. This confirms that the $RK4$ technique requires less computing cost than the Galerkin scheme.

5.   Example 2: HIV infection model description

In the field of HIV infection of healthy T-cells, mathematical models have become indispensable tools for better understanding the dynamical behaviour, transmission, disease progression, and immunological interactions of HIV. In humans, the HIV virus mainly targets healthy T-cells. Whenever HIV enters the body, it infects a significant number of healthy T-cells, causing their depletion gradually. As a consequence, the body's immune system is destroyed, weakening the host immunological response to viral diseases and eventually leading to acquired immunodeficiency syndrome (AIDS). The reduction in the number of these cells is used as a major indication to monitor the development of HIV infection and the symptoms of AIDS. These cells proliferate at a steady rate from bone marrow and thymus precursors. To characterize the dynamical behavior of HIV infection, several mathematical models have been introduced. The model proposed by Attaullah et al. ^[19] is as follows:

Herein, $\tilde{T}$, $\tilde{I}$ and $\tilde{V}$ describe the densities of healthy/infected T-cells and virus respectively. The explanation of parameters, along their values and initial condition for the state variables are illustrated in Table 18.

The aforementioned model is considered as a model problem, and we implement the suggested methods. Table 19 and Figure 5 show the numerical cost in terms of time, and for all other results the readers are advised to see ^[19] for details. It is worth noting that the $RK4$ approach required very less time than the $\mathrm{cGP}$(2) approach (see Table 19 and Figure 5). This is further supported by the fact that $RK4$ takes less time than $\mathrm{cGP}(2)$.

6.   Example 3: Model problem

In order to demonstrate that the provided temporal discretization technique is accurate and efficient in terms of time, as a test problem with the prescribed exact solution, analyze the implementations of the $\mathrm{cGP}$(2) and $RK4$ approaches. The detail of the considered problem is given on page 332 of ^[42]. The problem consists of two differential equation as follows:

where the initial values are $U(0) = 0$ and $V(0) = 0.$ The exact solution to the aforementioned model is

Performed different numerical tests to assess the accuracy of the suggested techniques utilized an equidistant and unequal time step size $\tau = T/N$ for all tests, with $T = 0.5.$ First, find out the solutions using the $\mathrm{cGP}$(2) and $RK4$ at $\tau = 0.1$ and then using at $\tau = 0.00125$ and $\tau = 0.000625,$ respectively. To determine the numerical solution's accuracy, the numerical solutions were compared with the prescribed analytical exact solution presented in Tables 20 and 21. The computational results clearly demonstrate that the Galerkin-scheme gives highly accurate results at relatively large time step sizes as compared to RK4-scheme. For the same step size ($\tau = 0.1$), the Galerkin-scheme gained the accuracy of $10^{-5}$ while the $RK4$-method achieved the accuracy of $10^{-4}.$ Further, $RK4$ achieved the accuracy of $10^{-5}$ at $\tau = 0.000625,$ while the Galerkin-method can already reach the same precision and accuracy with $\tau = 0.00125.$ This is a clear negation to the claim of ^[28] that Galerkin and $RK4$ have the same results. On the other hand, we find the order of convergence and compared the experimental orders of convergence with the theoretical order of convergence. The cGP(2)-method is of order 3 in the whole time interval and super-convergent of order 4 at the discrete time points. One can see from Table 24 that the experimental orders of convergence for the cGP(2) method are much more visible. From Table 24, we observe that the EOC coincides with the theoretical orders of convergence for corresponding time discretization schemes.

Moreover, we computed the required CPU-time in seconds at different time step sizes for the mentioned above system, as given in Table 25 and plotted in Figure 6 for both methods. It can be seen from Table 25 and Figure 6 that the $RK4$-scheme is much faster than the Galerkin-scheme. Hence, these numerical tests show that the Galerkin approach is less efficient and takes more time as in comparison to the $RK4$-method.

7.   Conclusions and future recommendations

In this work, a higher-order Galerkin and a fourth order classical Runge Kutta scheme were implemented. Various numerical tests were performed to determine the reliability and effectiveness of the described schemes. In addition, the impacts of increasing the parameters $\xi_1$ and $\psi_2$ on the dynamical behavior of the state variables were appropriately discussed. We noticed that increasing the value of $\xi_1$ increases the population of susceptible individuals while decreasing the population of infected individuals with acute HBV, infected individuals with chronic HBV and recovered individuals. By elevating $\psi_2$, the concentration of susceptible and recovered individuals grows, while the number of acutely and chronically infected people decreases. We executed many computational tests and compared the results quantitatively. According to computational simulations, the Galerkin-scheme provides much more accurate solutions at relatively high time step sizes than the $RK4$-method. Afterwards, the accuracy of the cGP(2) technique confirm that it produces more accurate results at relatively longer step sizes than the RK-4 scheme at a comparable numerical cost. The results of numerical studies revealed that the explicit $RK4$-scheme used less CPU time in comparison to the accuracy, but the implicit Galerkin-scheme required more CPU time. We proved that cGP(2) works better on time steps than the RK4 method. The proposed scheme could be applicable for solving complex real-world problems.

In the future, we plan to investigate the numerical and qualitative aspects of the present concept using different fractional order derivatives. The concept of fractional analysis is well known, and it has recently become a prominent study area. The utility of fractional calculus in modeling a wide range of real-world phenomena has been demonstrated. Using fractional order derivatives and integrals, researchers have researched infectious illnesses such as HIV, AIDS, and others.

Acknowledgments

This research received funding support from the NSRF via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation, (grant number B05F650018).

Conflict of interest

The authors declare that they have no conflicts of interest.