{'fraction_non_alphanumeric': 0.03739061256961018, 'fraction_numerical': 0.016131883673649785, 'mean_word_length': 4.581401085347805, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical solution of COVID-19 pandemic model via finite difference and meshless techniques

In the present paper, a reaction–diffusion epidemic mathematical model is proposed for analysis of the transmission mechanism of the novel coronavirus disease 2019 (COVID-19). The mathematical model contains six-time and space-dependent classes, namely; Susceptible, Exposed, Asymptomatically infected, Symptomatic infected, Quarantine, and Recovered or Removed (SEQIaIsR). The threshold number R0 is calculated by utilizing the next-generation matrix approach. In addition to the simple explicit procedure, the mathematical epidemiological model with diffusion is simulated through the operator splitting approach based on finite difference and meshless methods. Stability analysis of the disease free and endemic equilibrium points of the model is investigated. Simulation results of the model with and without diffusion are presented in detail. A comparison of the obtained numerical results of both the models is performed in the absence of an exact solution. The correctness of the solution is verified through mutual comparison and partly, via theoretical analysis as well.


Introduction
In December 2019, a novel virus COVID-19 was discovered. The first outbreak of the virus was reported in Wuhan, Hubei Province, China, with animals as a suspected source for the transmission to humans [1]. Subsequently, due to rapid human-to-human interactions, the transmission of the virus spread to many countries around the world, causing the World Health Organization (WHO) to declare a global health emergency in late January 2020 [2]. Many countries were put on high alert due to the pandemic, which has resulted in insufficient medical contingencies needed to cope with the rising number of the infected cases. In addition to huge economic losses, shortages of medicines and medical gadgets caused enormous human fatalities around the world. According to the WHO report [3], on 22 February 2022, the total number of infected cases reported worldwide was 424,822,073 and deaths were over 4,109,303.
The COVID-19 usually attacks the respiratory system, creating numerous severe health complications. Infected people may have numerous symptoms, for example, continuous cough, difficulty in breathing, and/or high fever [4]. Typical symptoms of people who get COVID-19 include coughing, shortness of breath, sneezing, fever as well as complications such as pneumonia and sore throat [5]. The incubation R. Zarin et al. Using mathematical tools, a real situation is represented in a mathematical formulation to mimic the real situation and predict the possible worst case scenarios. Therefore, this research aims to implement a mathematical model to study the different dynamics of the COVID-19 pandemic. For the numerical solutions of these models, various numerical methods have been reported in the literature. These include, Finite difference method [7], Fourier spectral method (FSM) [8], Finite element method (FEM) [9], and Finite difference method (FDM) [10]. Considerable literature is available on the models, methods, and their applications. Each method has its own advantages and shortcomings in terms of computational cost, accuracy, and convergence. The FDM provides reasonable accuracy on simple geometry but is not recommended for complex geometries. The FEM provides an accurate solution both on simple and complex geometries but has the drawback of construction of an efficient geometry fitted mesh. Global meshless methods based on radial basis functions (RBFs) have the flexibility of providing accurate solutions on complex geometries and are easily extendible to higher dimensions but it has a dense matrix representation. Meshless collocation methods are built upon a radially symmetric function called RBFs. These RBFs can be with or without shape parameter. Finding an optimal value of the shape parameter is still an open problem and solution of the full dense coefficient matrix is challenging. Meshless methods based on Multiquadric (MQ), Inverse multiquadric (IMQ), Gaussian (GS), and Inverse quadratic (IQ) are some of the examples [11].
The meshless approach based on RBFs appeared in the literature frequently for the numerical solution of different state-of-the-art PDEs [12][13][14][15]. In 1990, Kansa established a collocation method based on MQ RBFs to interpolate random data [16] and for approximation of PDEs [17]. Afterward, this approach got attention due to its meshless nature and ability to handle complex geometries in higher dimensions. Subsequently, the effectiveness of the meshless methods enhanced manifolds as a result of the localization of the meshless collocation methods [18][19][20][21]. The key feature of this approach is interpolation on overlapping sub-domains and ease of inversion of a small matrix. Other contributions towards using the local meshless method are shown in [22,23] and the references therein.
Grid is not involved in the meshless formulation, which is the main attribute of this method. Meshless approximation utilizes only one geometric property i.e., pairwise distances between the grid points. In any number of spatial dimensions, distances can be easily computed, which is considered to be one of the distinctive benefits of the meshless collocation method. Due to this added advantage, dealing in the higher dimensions is relatively easier in the case of meshless methods. Not only this but the meshless collocation method can also be used in those situations in which data points are scattered. Unlike the conventional methods, the coefficient matrix is dense and this is one of the main drawbacks of the global meshless method. The coefficient matrix of the global meshless method is not well-conditioned, and if not treated properly with respect to the appropriate selection of the shape parameter, it may lead to inaccurate inversion of the system matrix.
The novelty of the paper is the addition of the diffusion term to convert a system of ordinary differential equations model [24] to a set of partial differential equations model. Three methods: the FDOSM, the MOSM and the MEM, are employed for the numerical solution of the system of PDEs alongwith different sets initial conditions. This is an important addition to the existing stationary model [24] in the context of Covid-19 infection. Comparative performance of the state-of the art meshless method with finite difference method for the COVID-19 model is performed and the idea of Operator Splitting is extended from the finite difference method to the meshless method. Furthermore, both the stationary and dynamic models are considered for variety of initial conditions. In the case of ODE model, the simulation results often do not match with the real expectations. We have tested both the models for each case and have shown that in the case of variable initial conditions, the without diffusion model does not show any progression in the infection of COVID-19. Contrary to the real situation, it is observed during this experimentation that the infection spread is not witnessed by the ODE model in the case of variable initial conditions, even after the passage of 500 days.
Classification of the paper is as follows: Section 2 is devoted to the mathematical model formulation. The equilibrium points and reproductive number are obtained in Section 3. In Section 4, a stability analysis for the system is presented. For the numerical solutions of the system, a meshless explicit method (MEM), a finite difference operator splitting method (FDOSM), and meshless operator splitting method (MOSM) are presented in Section 5. In Section 6, three sets of different initial conditions are presented to test the efficacy of the methods. Graphical simulation results of the proposed numerical methods are presented in Section 7. Concluding remarks are provided in the last section.

Model formulation
In the literature, the existing COVID-19 virus infection models are mostly stationary [24,25]. However, more realistic models must include space dependency while modeling the spread of the virus to control disease spread through the mobility of individuals in space. Ahmed, I. et al. [24] proposed the following lumped parameter COVID-19 model: The initial conditions (ICs) for the above model are given as follows: In the proposed SEQI I R pandemic model the total population ( ) is classified into Susceptible ( ), Exposed ( ), without symptoms (asymptomatic) infected individuals ( ), with symptoms (symptomatic) infected individuals ( ), Quarantine ( ) and Recovered ( ) individuals. The parameters are described as: • * → Recruitment rate • * → Contact rate between the exposed and susceptible • * → The transmission rate of those individuals who come from exposed to asymptomatic infected class • * → The transmission rate of those individuals who come from exposed to quarantine class • * → Coronavirus related death rate in the class of symptomatic infected • * → The transmission rate of those individuals who come from quarantine to symptomatic infected class • * → The transmission rate of those individuals who join symptomatic infected class from exposed class • * → The transmission rate of those individuals who join the asymptomatic infected class from quarantine class • * → Mortality rate (natural) • * → The transmission rate of individuals who come to quarantine from susceptible class • * 1 → Recovery rate in asymptomatic infected class • * 2 → Recovery rate in symptomatic infected class The model (1) is built on the assumption of homogeneous population, which means that mixing of individuals is considered in such a way that there is no distinction between individuals in one place and individuals somewhere else in the domain. In genuine situations, the infection may spread quicker in one spot than in another, due to various conditions like diverse climate, etc. Consequently, it is fundamental for the model to rely upon space variable too. Hence by adding a diffusion term to (1), we get the following reaction-diffusion system: where represents time variable and X represents space variable. The initial conditions (ICs) for the above model are as follows: No-flux boundary conditions (BCs) for the model (5) where represent a smooth boundary with a bounded domain ⊂ R, and homogeneous Neumann boundary conditions mean that no individual crosses the boundary . The coefficients * 1 , . . . , * 6 are positive constants, which control movement or migration of susceptible, exposed, quarantine, asymptotically infected, symptomatic infected, and recovered classes of the population, respectively. In the present work, we always assume that * (X), * (X), * (X), * (X), * (X), * (X), * (X), * (X), * (X), * 1 (X) and * 2 (X) are all positive, continuous and bounded functions on .

Equilibrium points and threshold number
The spread and control of the infection are fundamentally connected with the threshold number 0 called the reproductive number. The number 0 predicts the possibilities of controlling the disease or otherwise. If 0 < 1, the infection vanishes from the community and the disease free equilibrium (DFE) exists, which is asymptotically stable both locally and globally. If 0 > 1, then the endemic equilibrium (EE) is stable both locally and globally under some specific conditions. In such a scenario, the infection dwells in the population for a long time. The disease free equilibrium (DFE) point for model (5) is given as follows: To find 0 for the proposed model (5), we use the next generation matrix technique [26,27]. Let us take the following classes of system (5) For simplicity we can write system (6) as; Linearizing the model (7) at steady state 00 , we get the following linear system: where 1 , 2 , 3 and 4 are the linearized form of , , and . Since 3 and 4 do not appear in first and second equations, we only need to consider the sub-system as follows We take .
From the expressions of and −1 , we obtain Thus, we have ) .

Stability analysis
In this section, we discuss asymptotic stability analysis of system (5) at endemic equilibrium point * as well as at disease-free point 0 .

Numerical stability of endemic equilibrium
Fourier series representation of (15) can be written as: where = 2 , ( = 1, 2, 3, …). Using these expressions in Eq. (15) we get: Box II.  The variational matrix 1 for the system (5) is given by (see Box II), where The determinant det where + ( * + * + * + * ), The eigenvalues of the characteristic equation are: Further, for 0 > 1, the quadratic equation [̄2 + 1̄+ 0 ] = 0 has positive coefficients and therefore, the real part of the reaming two eigenvalues must be negative, which implies that̄5 < 0 and̄6 < 0. Thus, for 0 > 1 all the eigenvalues of the variational matrix 1 are negative at the EE point * * , hence * * is stable. This completes the proof of the above theorem.   By evaluating the above determinant we obtained the following equation;

Stability of disease-free equilibrium
( The eigenvalues of the characteristic equation are: For 0 < 1, it is clear from Eq. (26) that the real part of each eigenvalue is negative. Thus, for 0 < 1 all the eigenvalues of the variational matrix are negative at the DFE point 00 , hence 00 is stable this completes the proof of the above theorem.  Fig. 2, confirming the fact that the solution variables are approaching to the equilibrium point 00 = ( 0 , 0 , 0 , 0 , 0 , 0 ) = (2.3491, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000).

Discretization of space and time variables
In the present section, we utilize MEM, FDOSM and MOSM for the numerical solutions of the system (5). In all cases, time step = 0.024 days and spatial step X = 1 meter are used.

Explicit discretization of time based on forward difference operator
Forward difference approximation of the first-order is applied on the system (5), which gives:

Discretization of time based on operator splitting
In this case the time discretization is carried out in two phases. In the first phase the first-order time derivative is discretized from to + * as: In the second phase, the time discretization is performed from + * to +1 as follows; where * = 2 .

Discretization of the space variable based on finite difference
In this case a second-order finite difference scheme is used for second order space derivative present in Eq. (29)

Discretization of space variable based on MQ RBF
Considering centers X 1 , X 2 , … , X , ∈ , we have Where ( ) is an RBF and = √ ( X − X ) 2 , for , = 1, 2, … , . The unknowns , = 1, … , are obtained from the interpolation condition at the set of nodal points X , = 1, 2, … , is R. Zarin et al.  As the centers and the collocation points coincide, hence The coefficient vector = ( 1 , 2 , … , ) is an × 1 vector and the matrix is an × interpolation matrix having entries of the MQ RBF, and The RBFs estimate for the derivatives of (X) can be represented as X, X ∈ (0, 5, ] R. Zarin et al.  Table 2.  Table 2.

Initial conditions
For the numerical simulation of the model, we have taken three sets of different types of initial conditions (ICs) to test the effectiveness of proposed numerical techniques. These conditions are given in Table 2, and the plots of different ICs are given in Figs. 3 and 4.

Numerical simulations
We consider the initial conditions in Table 2. For the numerical solution of (5), the parameters values given in Section 4 are considered. Analysis are conducted at different points within the domain [−10, 10] for 600 days. The shape parameter value in the MQ RBF is chosen as = 100 [29]. Value of the diffusion coefficients are taken 1 = 2 = 0.05, 3 = 6 = 0, 4 = 0.1, 5 = 0.02, unless stated otherwise. The spatial discretization points are taken 20, whereas the time steps are chosen 25000.

Initial condition-1
In this case we consider the IC-1 given in Table 2 at X = + 15 X. Fig. 5 shows the simulation results of the FDOSM with and without diffusion. The simulation results of the model with and without diffusion via the FDOSM are almost identical in the case of constant initial population. The MOSM and the MEM simulation results corresponding to the same initial conditions are depicted in Fig. 6. Accordingly, the methods MEM and MOSM yield identical numerical results.

Initial condition-2
The results produced by the FDOSM corresponding to IC-2 of Table 2 with and without diffusion are displayed in Fig. 7. The results of the MEM and the MOSM corresponding to the same initial condition are displayed in Fig. 8. It is clear from the figures that values of all the classes, especially, the class are higher in without diffusion case versus the diffusion case. This is due to the mobility of the population, R. Zarin et al.  Table 2, at X = + 15 X.  Table 2, at X = + 15 X.
which reduces the chances of concentration of the population in the case of diffusion. The mobility of the population allows to create social distancing among the concentrated population, which is the main factor in the strategy to contain the COVID-19 spread. The simulation results of the diffusion based model is different, showing steady-state behavior in the last 100 days. In the case of no diffusion model, considerable variation is recorded in that period of time, which seems unrealistic as only susceptible class of population is dominant for the first 500 day. This implies that there is a role of diffusion which has impacted the results over the passage of time. A more realistic understanding of the disease dynamics can be obtained if two-and three-dimensional diffusion models are incorporated for this purpose.

Initial condition-3 at X =
The simulation results of the FDOSM corresponding to initial condition 3 of Table 2 with and without diffusion models are shown in Fig. 9. Similarly, The numerical results of the MOSM and the MEM corresponding to the same initial condition are shown in Fig. 10. In this case, there is a higher concentration of population at X = . As shown in Figs. 9 and 10, it is clear that the value of the class is lower in the vicinity of = 600 days for the diffusion case due to the diffusion mobility of the population at X = . Also, it is evident from the same figures that there is a rapid increase in graphs of the classes , , , and for the same period. This may be due to the diffusion effects, which create mobility for the individuals, thus leading to a low number of susceptible individuals and a high number of individuals in other classes. However, a realistic picture of the epidemic disease can be captured once high-dimensional diffusion models are incorporated for this purpose.

+ 4 )
The simulation results of the FDOSM corresponding to initial condition 3 of Table 2 with and without diffusion models are shown in Fig. 11. Similarly, The numerical results of the MOSM and the MEM R. Zarin et al.  Table 2, at X = .  Table 2, at X = .
corresponding to the same initial condition are shown in Fig. 12. In this case, the simulation results have different patterns for with and without diffusion cases. The simulation results of the without diffusion case are not realistic as the infected class of population stay near zero even after the passage of 600 days. In contrast, the diffusion based model is more realistic as the peaks in the infected population results a basin in the susceptible of the class population. It can also be seen from Figs. 11 and 12 that there is slight variation in the simulation results of the asymptotic class of population.

Initial condition-3 at
Numerical solutions of the models with and without diffusion corresponding to the IC-3 of Table 2 are displayed in Fig. 13. While the numerical solutions of the diffusion model corresponding to the same initial condition produced by the MEM and the MOSM are shown in Fig. 14. In this case, the population is less dense in the vicinity of X = ( ) . However, when the diffusion term is added, the values of the classes , , , and decrease, and there is little increase in the value of the exposed class due to mixing of populations.

Conclusions
The purpose of epidemiological modeling is to predict the progression of the disease and to plan control measures. Forecasting based on epidemiological models of infection spread can help plan various public health interventions by the government such as lockdowns and administration of vaccines to the population. In this paper, we reported the SEQI I R model to predict the growth of the COVID-19 infected population, based on some theoratical initial population. The model takes into account both the dynamic behavior and spatial mobility of susceptible, exposed, infected, quarantined, and recovered classes of the population. The spatial mobility of the population is restricted to a single spatial dimension. The focus of study is to check performance of the model with and without diffusion. In addition to finite difference method, meshless collocation method based on MQ RBF is used to Fig. 11. Outcomes of the FDOSM for the IC-3 of Table 2, at X = (
approximate spatial derivatives and forward difference operator is used for approximating the time derivative operator of the model. This paper presents both distributed parameter models and lumped parameter models to predict the spread of the epidemic diseases such as COVID-19. The highly non-linear SEQI I R models are simulated by the three numerical methods namely; the FDOSM, the MEM, and the MOSM. The performance of these methods has been mutually compared and theoretical analysis of the model in terms of equilibrium points and stability have been investigated. The take aways of the proposed work are listed below: • The diffusion factor bears much impact on the final simulation results on the distributed parameter model of the three methods in the case of piece-wise initial populations with an initial concentration in certain areas of the domain. • The simulation results of the diffusion-added model are more realistic than the numerical results of the without diffusion model.
• The scope of the meshless methods has been extended to COVID-19 mathematical model. • Comparison and agreement of the meshless methods with standard finite difference method further ascertains validity of the meshless method in the context of COVID-19 mathematical model.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
No data was used for the research described in the article.  Table 2, at X = ( ) .  Table 2, at X = ( ) .