A Numerical Investigation Based on Exponential Collocation Method for Nonlinear SITR Model of COVID-19

1  Introduction

The world has recently contracted the dangerous and deadly disease coronavirus 2019 (COVID-19), which is an almost uncontrollable respiratory infection. The new coronavirus was reported in Wuhan, China in December 2019. The World Health Organization (WHO) announced the global pandemic of the infection in March 2020. Its disease pattern is called COVID-19 [1]. This disease that is highly contagious and through droplets can be transmitted from person to person. Many people are infected with COVID-19 and the number of victims of this virus is increasing every day. The COVID-19 spread ratio is very problematic and a concern for the whole world. The main cause of the spread of this virus is the contact of an infected person with healthy people, because studies have shown that this infection is usually caused by the transmission of blood cells through coughing or sneezing. These blood cells can stay in the air for a long time and cause infections for others. However, it is very challenging for scientists to study the preventive measures that can be taken to control the spread of the virus and to develop a vaccine to fight the virus. How to control the disease is one of the most important challenges if it is transmitted to animals or birds. Researchers are developing a variety of tactics to study the growing behavior of the coronavirus, and are studying ways to end COVID-19. Much research has been done to determine the conditions under which this deadly virus can be controlled. Scientists have found that COVID-19 is one of the most important outbreaks that attacks the respiratory system. One of the main reasons for the prevalence of COVID-19 is due to the transmission of germs through the respiratory cells in humans, and this virus is considered as a vector of transmission. The World Health Organization (WHO) has warned that the outbreak of the coronavirus could spread more rapidly if control measures were not implemented in a timely manner. WHO has advised people to stay away from infectious people or animals, including fever or any respiratory problems, and to recommend the use of surgical masks and the continued use of hand sanitizers in public places to protect oneself from infection. Social distance or less crowded places can reduce the risk of spreading the corona virus because it is more likely to spread in crowded places. The use of epidemiological computational simulation models plays a main role in estimating the transmission parameters and guessing the effective behavior of the infection. These models have been shown to be useful in depicting the growth rate or rate of decay of viruses over time. And they are very useful in providing control measures that can be adapted to reduce the spread of the disease. Many results and applications of COVID-19 simulation models are accepted and published [2–9]. Providing a mathematical model is another matter, while predicting the consequences of the disease is also very disappointing. One method of analyzing the behavior of diseases is divisional modeling, which can be used for mathematical models related to affective diseases. An appropriate vaccine must be discovered to prevent further damage of the virus. Otherwise, the world must be prepared to face many new challenges such as human casualties, food shortages, poverty, unemployment and so on.

One way to study the effective growth of this deadly outbreak is to use computer simulation block models. In different researches, different numerical and analytical approaches have been considered to perform the outputs of SITR models [10–13].

Differential equations have a remarkable role in several scientific and engineering phenomena that have always been considered at physical and technical applications and they are appeared in various areas as mathematics, physics and engineering sciences [14–17].

In recent years, Yüzbaşi et al. have applied the collocation method based on exponential approximation to solve some problems like pantograph equation, the linear neutral delay differential, Fredholm integro-differential difference equations and so on [18–21].

The aim of the present work is to find a numerical solution for a nonlinear SITR model using the new dynamic parameters of COVID-19, along with numerical analysis to better understand the expansion using numerical approaches through exponential basis.

The organization of this article is structured as follows: In Section 2, mathematical modelling of the SITR system is presented. In Section 3, we express briefly required mathematical elementary and matrix relations for exponential functions of the method. In Section 4, we present the numerical implementation of the matrix operation of the method. Section 5 involves the error analysis. For this purpose, assuming the residual of solution of Eq. (1). Section 6 contains numerical examples, where approximate solutions corresponding to various N values are obtained using the proposed method. Numerical experiments are examined to illustrate the efficiency and accuracy of the method, and results are reported. Finally, the last section consists of a brief conclusion.

2  Mathematical Modeling

Mathematical modelling of the dynamical system is an interesting field of study that has attracted most researchers’ attention. Dynamical systems have a wide range of applications, including population growth models, biomedical research, biological systems, and engineering.

In this approach, for the SITR model is designed for the new COVID-19 dynamics, the population is allocated to various compartments with particular labels, susceptible (S), infected (I), treatment (T), and recovered (R). The S is divided into S1(t) and S2(t) subclasses. Subclass S1(t) identifies people who have not yet developed COVID-19 disease, and subclass S2(t) identifies people who are not yet infected but have some serious illnesses or older. No proper vaccination for the disease has been discovered yet, so the entire individuals are at risk of getting contaminated by this infection. That is why the whole population is often put into the susceptible class. From the susceptible class, the entities of that class can get infected by making contact with the infectious person and then they can join the infected class as well as exposed class.

The SITR model of the novel COVID-19 dynamics is given by the following nonlinear 1st order differential equations and the description of each compartment is given in Table 1.

{S1′(t)=B−βI(t)S1(t)−δβT(t)−αS1(t),S1(0)=φ1,S2′(t)=B−βI(t)S2(t)−δβT(t)−αS2(t),S2(0)=φ2,I′(t)=−μI(t)+βI(t)(S1(t)+S2(t))−αI(t)+βδT(t)+σI(t),I(0)=φ3,T′(t)=μI(t)−ρT(t)−αT(t)+ψT(t)+εT(t),T(0)=φ4,R′(t)=−αR(t)+ρT(t),R(0)=φ5.(1)

Each equation is describing the transmission behavior of individuals in the respective compartments. By this transmission number of individuals can vary in each of the five compartments [22–30]. The description of the transition rates in each cell is given in Table 2. The parameter of interest with appropriate settings for particular descriptions of the dynamics are provided in Tables 1 and 2.

We try to find the approximate solution of system of differential Eq. (1) with initial conditions as series of exponential functions. Exponential functions or exponential polynomials are based on the linearly independent exponential basis set

ℬ={1,e−t,e−2t,…}.(2)

To begin with, we assume that the unique solutions of system (1) can be expressed as a exponential series of the form

S1(t)=∑n=0∞ane−nt,S2(t)=∑n=0∞bne−nt,I(t)=∑n=0∞cne−nt,

T(t)=∑n=0∞dne−nt,R(t)=∑n=0∞ene−nt,(3)

then by truncation these power series after the (N+1)st term, their linear combinations defined by the expansions

S1,N(t)=∑n=0Nane−nt,S2,N(t)=∑n=0Nbne−nt,IN(t)=∑n=0Ncne−nt,

TN(t)=∑n=0Ndne−nt,RN(t)=∑n=0Nene−nt,(4)

which, coefficients an,bn,cn,dn,en are unknown and N is an arbitrary positive integer.

Remark. We solve the system (Sys) in time interval [0, 1] and for larger intervals [0, R], convert it to interval [0, 1] by using an iterative method. This iterative process is such that to solve the system in t∈[0,R] by initial conditions (IC) at t=0, we first solve the system in the interval t∈[0,1]. Then we use the results for t=1 with a transfer as the initial conditions and solve the system at [1, 2]. And we continue this iterative procedure until t=R by step size h=1.

{Sys(t)=0,IC(0)=k0,Sys(t−1)=0,IC(1)=k1,⋮Sys(t−(R−1))=0,IC(R−1)=kR−1.(5)

After solving these systems, we merge the functions as multivariate function.

3  Preliminaries and Matrix Relations

In this section, we outline operational matrices of the exponential method we will use in order to solve system (1).

In the first step, we create the differentiation matrices which are the basic tools of the current approach. Differentiation matrices make this method more suitable for managing high-order differential equations. By constructing an operational matrix, it is easy to derive high-order derivatives of the unknown in terms of values at collocation points.

Firstly, we inscribe the approximated solution uN(t) defined by linear combination (4) of Eq. (1) in the matrix form as

u(t)=E(t)A,(6)

where

E(t)=[1e−te−2t…e−Nt],

and

A=[a1a2…aN]T.

Taking advantage of the linearity of expansion (4), we can compute the derivative of u by differentiating the basic functions. The derivatives of u are obtained as follows:

u′(t)=∑n=1N(−nan)e−nt,(7)

and for higher order derivatives of u we present a matrix form. Next, we explain how to create a differentiation matrix through the method, and we extract and create a matrix D so that the equations are in the collocation points.

The derivative of the approximate solution can also be expressed as a product of matrices. Namely, E(t) has a relation with its first derivative E′(t) that is demonstrated by

E′(t)=E(t)D,

where the operational (differentiation) matrix D corresponding to above relation is represented

D=[000⋯00−10⋯000−2⋯0⋮⋮⋮⋱⋮000…−N],

and that, after repeating the procedure k

Ek(t)=E(t)Dk,k=0,1,2,…,(8)

holds for any nonnegative integer k, that D0 is the identity matrix which its dimension is: (N+1)×(N+1).

Note that

Dk=[000⋯00(−1)k0⋯000(−2)k⋯0⋮⋮⋮⋱⋮000…(−N)k].

By using of the matrix relations (6) and (8), we can write matrix representation as

u(k)(t)=E(t)DkA,k=0,1,…,m.(9)

After replacing the collocation points {xi}i=0N in (9), we turn into the following system of matrix equations as

u(k)(ti)=E(ti)DkA,i=0,1,…,N,

that, in the matrix form, we have

U(k)=E.Dk.A,

where

E=[E(t0)E(t1)⋮E(tN)]=[1e−t0e−2t0…e−Nt01e−t1e−2t1…e−Nt1⋮⋮⋮⋱⋮1e−tNe−2tN…e−NtN],U(k)=[u(k)(t0)u(k)(t1)⋮u(k)(tN)].

4  Implementation of Matrix Operation

In this section, we explain how to use the exponential collocation method for problem (1). For this purpose, we use the following procedure. The basis of this method is based on the calculation of unknown coefficients using the collocation points.

To acquire an exponential series solution of Eq. (1) under the conditions, the operational matrix method is applied as follows:

S1,N(t)=E(t)Aa,S2,N(t)=E(t)Ab,IN(t)=E(t)Ac,

TN(t)=E(t)Ad,RN(t)=E(t)Ae,(10)

where

E(t)=[1e−te−2t…e−Nt],

and

Aa=[a1a2⋯aN]T,Ab=[b1b2⋯bN]T,Ac=[c1c2⋯cN]T,

Ad=[d1d2⋯dN]T,Ae=[e1e2⋯eN]T.

The relationship between E(t) and its first derivative can also be written

S1,N′(t)=E′(t)Aa,S2,N′(t)=E′(t)Ab,IN′(t)=E′(t)Ac,

TN′(t)=E′(t)Ad,RN′(t)=E′(t)Ae,(11)

that by applying operational matrices, we obtain the following matrix forms:

S1,N′(t)=E(t)DAa,S2,N′(t)=E(t)DAb,IN′(t)=E(t)DAc,

TN′(t)=E(t)DAd,RN′(t)=E(t)DAe.(12)

Now let’s create the U(t) and U′(t) matrices as follows:

U(t)=E¯(t)A¯,andU′(t)=E¯(t)D¯A¯,(13)

where

U(t)=[S1(t)S2(t)I(t)T(t)R(t)],U′(t)=[S1′(t)S2′(t)I′(t)T′(t)R′(t)],E¯(t)=[E(t)00000E(t)00000E(t)00000E(t)00000E(t)],

and

D¯=[D00000D00000D00000D00000D],A¯=[AaAbAcAdAe].

The system (1) can be expressed as

{S1′(t)+αS1(t)+βI(t)S1(t)+δβT(t)=B,S2′(t)+αS2(t)+βI(t)S2(t)+δβT(t)=B,I′(t)+(μ+α−σ)I(t)−βI(t)(S1(t)+S2(t))−βδT(t)=0,T′(t)+(ρ+α−ψ−ε)T(t)−μI(t)=0,R′(t)+αR(t)−ρT(t)=0,(14)

that it can be rewritten as following matrix form:

U′(t)+PU(t)+MV1(t)+NV2(t)=q,(15)

where

U(t)=[S1(t)S2(t)I(t)T(t)R(t)],U′(t)=[S1′(t)S2′(t)I′(t)T′(t)R′(t)],P=[α00δβ00α0δβ000μ+α−σ−βδ000−μρ+α−ψ−ε0000−ρα],

and

V1(t)=[I(t)S1(t)],V2(t)=[I(t)S2(t)],M=[β0−β00],N=[0β−β00],q=[BB000].

To determine the unknown coefficients, we use the collocation points

ti=RNi,i=0,1,…,N.(0≤ti≤R),(16)

and by replacing the collocation points into Eq. (15), we obtain a system of matrix equations as follows:

U′(ti)+PU(ti)+MV1(ti)+NV2(ti)=q,i=0,1,…,N,(17)

and this system corresponds to the following matrix form:

U′¯+P¯U¯+M¯V¯1+N¯V¯2=Q,(18)

where

U′¯=[U′(t0)U′(t1)⋮U′(tN)],P¯=[P0…00P…0⋮⋮⋱⋮00…P],U¯=[U(t0)U(t1)⋮U(tN)],

and

M¯=[M0…00M…0⋮⋮⋱⋮00…M],N¯=[N0…00N…0⋮⋮⋱⋮00…N],V¯1=[V1(t0)V1(t1)⋮V1(tN)]=I¯S¯1,

and

V¯2=[V2(t0)V2(t1)⋮V2(tN)]=I¯S¯2,Q=[qq⋮q].

Now, let us find the relations between the matrix A and other matrices in system (18) as follows:

U¯=E¯A¯,U′¯=E¯D¯A¯,V1¯=E¯AcE¯Aa,V2¯=E¯AcE¯Ab.(19)

After the substitution of the above relations, we obtain to the following matrix equation:

E¯D¯A¯+P¯E¯A¯+M¯E¯AcE¯Aa+N¯E¯AcE¯Ab=Q,

by previous definition A¯=[AaAbAcAdAe]T we obtain to

(E¯D¯+P¯E¯+M¯E¯AcE¯+N¯E¯AcE¯)A¯=Q.(20)

Briefly, (20) can also be presented as follows:

WA¯=Q,

where

W=E¯D¯+P¯E¯+M¯E¯AcE¯+N¯E¯AcE¯.(21)

Here, (20) conforms to a nonlinear system of the (N+1)×(N+1) algebraic equations with the unknown coefficients A¯.

By placing the colocation points at the algebraic equations and considering initial conditions, the system be solved and the coefficients be determined.

5  Error Analysis

Since most systems do not have exact solutions and most similar programs usually do not have the ability to generate numerical solutions, we must verify the accuracy of the numerical results using a method. We can confirm the accuracy and effectiveness of the numerical results by considering the residual error functions. The residual method is a general class of methods developed to obtain the approximate solution of differential equations. In the residual method, an approximate solution based on the overall behavior of the dependent variable is considered. The assumed solution is often chosen to satisfy the boundary conditions. This assumed solution is then replaced by the differential equation. Since the assumed solution is only approximate, it generally does not satisfy the differential equation and therefore leads to an error or what we call the residual. The residual then disappears in the mean sense throughout the solution domain to produce a system of algebraic equations.

If S1,N(t),S2,N(t),IN(t),TN(t) and RN(t) be numerical solutions of unknown functions, the residual error functions are defined as

{Res1(t)=|S1,N′(t)+αS1,N(t)+βIN(t)S1,N(t)+δβTN(t)−B|,Res2(t)=|S2,N′(t)+αS2,N(t)+βIN(t)S2,N(t)+δβTN(t)−B|,Res3(t)=|IN′(t)+(μ+α−σ)IN(t)−βIN(t)(S1,N(t)+S2,N(t))−βδTN(t)|,Res4(t)=|TN′(t)+(ρ+α−ψ−ε)TN(t)−μIN(t)|,Res5(t)=|RN′(t)+αRN(t)−ρTN(t)|.(22)

6  Numerical Applications

The rate of increase in the number of infections depends on the product of the number of infected and susceptible persons. The system (1) explains the dramatic increase in infection rates worldwide. The travel of infected individuals around the world has led to an increase in the number of infected people, which in turn has led to a further increase in the susceptible individuals. This creates a positive feedback loop that leads to a rapid increase in the number of actively infected persons. Therefore, during the period of increase, the number of susceptible people increases and consequently the number of infected people also increases.

In this section, we present some applications of the method that are presented in different values for the parameters β,μ and α as Contact Rate, Recovery Rate and Death Rate, respectively.

Case (I): β=0.25,μ=0.08,α=0.20.

Case (II): β=0.30,μ=0.10,α=0.25.

These results show the effectiveness of the present method in achieving good accuracy with fewer collocation points and less computational time. All the calculations have been performed using MAPLE. By applying the suggested method for N=4,6 in 0≤t≤1, we obtain the approximate solutions of S1(t),S2(t),I(t),T(t) and R(t). Different appropriate values of epidemic parameters are given in Table 3.

The residual of solutions has been used to show that this method is efficient and reasonably accurate.

The coefficients of the approximate numerical outcomes of the method for N=4 of Case (I) and Case (II) are given in Tables 4 and 6, respectively. The residual errors for Case (I) and Case (II) are reported in Tables 5 and 7, respectively. Fig. 1 shows the graph of residual error functions for Case (I).

Figure 1: Graph of residual errors for Case (I) for N=4 and N=6

In Fig. 2, graph of functions S1,N(t),S2,N(t),IN(t),TN(t) and RN(t) for Case (I) for N=4 and N=6 are presented that show the behaviour of main functions in SITR model. Fig. 3 shows the graph of residual error functions for Case (II). In Fig. 4, graph of functions S1,N(t),S2,N(t),IN(t),TN(t) and RN(t) for Case (I) for N=4 and N=6 are presented that show the behaviour of main functions in SITR model.

Figure 2: Graph of functions S1,N(t),S2,N(t),IN(t),TN(t) and RN(t) for Case (I) for N=4 and N=6

Figure 3: Graph of residual errors for Case (II) for N=4 and N=6

Figure 4: Graph of functions S1,N(t),S2,N(t),IN(t),TN(t) and RN(t) for Case (II) for N=4 and N=6

These graphs show that contact rates are primarily a source of large increase in the number of susceptible individuals, but begin to decrease over time. This is because higher contact rates cause more people to become infected and move to an infected class; therefore, the number of individuals in the susceptible class decreases. About S2(t), graphs clearly show that the number of people with the disease increases with increasing contact rate. Decreasing the contact rate leads to a small increase in the number of infected people, while increasing the contact rate suddenly increases the number of infected people. About I(t), graphs show the growing behavior of susceptible individuals with higher recovery rates. And about T(t), they show that as the recovery rate increases, so does the number of people in the infectious class. It can be noted that when the recovery rate is low, fewer people recover from the virus due to the death of infected people.

The figures also show that with a high death rate, the number of recovered individuals suddenly decreases. Because when the death rate is high, so many people in the affected and recovered classes lose their lives, resulting in a decrease in people in all the classes. Since the death rate is so high, the infected and recovered people are almost gone.

Figs. 5 and 6 show the graph of residual error functions for Case (I) and Case (II) on interval t∈[0,10], respectively. For solving the system on [0, 10], we apply iterative process (5). When the 10-day future change of COVID-19 is examined, it is observed from Figs. 2 and 4 that the number of infected individuals IN(t) from COVID-19 increased by decreasing until the 5th day (up to t=5) and then remained constant until the 10th day (t=10).

Figure 5: Graph of residuals for Case (I), N=6 and t∈[0,10]

Figure 6: Graph of residual errors for Case (II)N=6 and t∈[0,10]

7  Conclusion

In this paper, the exponential approximation is used to solve the numerical investigation of a nonlinear SITR model that represents the dynamics of the new disease of COVID-19. The method is based on exponential functions and the collocation method as an operational matrix. As seen, there is no concern about approximating higher-order derivatives of the unknowns. Also, to show the accuracy and efficiency of the method, two cases with different values have been examined. Through the examples provided, we realize that the obtained numerical results have very good residual errors. Moreover, it is realized that errors decrease when N values increase. As shown in the results obtained from computations, we conclude that implementation of this method will be very easy with less computational costs for similar problems. The results indicate the effectiveness of the present method in achieving good accuracy with fewer collocation points and less computational time.

Acknowledgement: The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.

Funding Statement: The authors received no specific funding for this study.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.