A Pell–Lucas Collocation Approach for an SIR Model on the Spread of the Novel Coronavirus (SARS CoV-2) Pandemic: The Case of Turkey

Abstract

In this article, we present a study about the evolution of the COVID-19 pandemic in Turkey. The modelling of a new virus named SARS-CoV-2 is considered by an SIR model consisting of a nonlinear system of differential equations. A collocation approach based on the Pell–Lucas polynomials is studied to get the approximate solutions of this model. First, the approximate solution in forms of the truncated Pell–Lucas polynomials are written in matrix forms. By utilizing the collocation points and the matrix relations, the considered model is converted to a system of the nonlinear algebraic equations. By solving this system, the unknown coefficients of the assumed Pell–Lucas polynomial solutions are determined, and so the approximate solutions are obtained. Secondly, two theorems about the error analysis are given and proved. The applications of the methods are made by using a code written in MATLAB. The parameters and the initial conditions of the model are determined according to the reported data from the Turkey Ministry of Health. Finally, the approximate solutions and the absolute error functions are visualized. To demonstrate the effectiveness of the method, our approximate solutions are compared with the approximate solutions obtained by the Runge–Kutta method. The reliable results are obtained from numerical results and comparisons. Thanks to this study, the tendencies of the pandemic can be estimated. In addition, the method can be applied to other countries after some necessary arrangements.

1. Introduction

In December 2019, an epidemic first appeared in Wuhan, China’s Hubei province. The cause of this epidemic was not clear, and the epidemic quickly spread to other countries. Not long after, this infectious disease of unknown cause was identified as a new coronavirus (nCoV) and this virus was named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The World Health Organization (WHO) named this infectious disease as coronavirus disease 2019 (COVID-19) and the SARS-CoV-2 epidemic was declared a pandemic on 11 March 2020. According to worldometer data, as of 25 December 2022, worldwide, there have been a total of 661,711,220 cases, 6,685,775 deaths, and 634,178,985 recoveries.

To address COVID-19, measures such as the mutual stoppage of countries’ flights, border closings, taking quarantine decisions for infected people, curfews, education suspension, and the beginning of distance education were taken. In addition, all kinds of cultural, scientific, artistic, and similar meetings and events were postponed. Places such as theatres, cinemas, massage parlors, gyms, cafes, concert halls and wedding halls were temporarily closed. Simultaneously, scientists started the vaccine studies and soon after, people tried to immunize the population by vaccinating them. Thus, the normalization process was begun. However, the number of cases and deaths is still increasing significantly. For this reason, all studies related to the pandemic are of great importance for science and humanity.

On the other hand, studies were also started in the field of mathematics for this pandemic with the help of the models related to infectious diseases. The most important of these model problems are the continuous population models [1,2,3,4,5,6,7], the Lotka–Volterra population model [2,5,6,7,8,9,10,11,12,13,14], the Hantavirus infection model [15,16,17,18,19], the HIV infection models [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35], the SIR epidemic model [36,37,38,39,40,41], and the SIRD epidemic model [42,43,44,45].

Canto, Avila–Vales and Garcia–Almeida studied a SIRD-based COVID-19 models in Yucatan, Mexico in 2020 [46]. Canto and Avila–Vales worked on a parametric estimation of an SEIR and an SIRD models of COVID-19 pandemic in Mexico in 2020 [47]. Calafiore, Novara, and Possieri investigated a modified SIR model for the COVID-19 contagion in Italy in 2020 [48]. Calafiore and Novara studied a time-varying SIRD model for the COVID-19 contagion in Italy in 2020 [49]. Mohammadi, Rezapour, and Jajarmi worked the fractional SIRD mathematical model for the first and second waves of the disease in Iran and Japan in 2021 [50]. Pacheco and Lacerda made function estimation and regularization in an SIRD model applied to the COVID-19 pandemics in 2021 [51]. Faruk and Kar conducted a data-driven analysis and prediction of COVID-19 dynamics during the third wave by using an SIRD model in Bangladesh in 2021 [52]. Covid-19 epidemic data in Italy, using an adjusted time-dependent SIRD model, was modeled by Ferrari et al. in 2021 [53]. Kovalnogov, Simos, and Tsitouras studied Runge–Kutta pairs suited for SIR-type epidemic models in 2021 [54]. Martinez investigated a modified SIRD model to study the evolution of the COVID-19 pandemic in Spain in 2021 [55]. Pei and Zhang made long-term predictions of COVID-19 in some countries by a SIRD Model in 2021 [56]. The progress of the COVID-19 outbreak in India was worked by Chatterjee et al. in 2021 [57] by using a SIRD model. Fern$\acute{\mathrm{a}}$ndez–Villaverde and Jones estimated and simulated an SIRD model of COVID-19 for many countries, states, and cities in 2022 [58]. In addition, there are some studies in the literature regarding these models [59,60,61,62,63].

In 2020, a novel parametric model of the COVID-19 to estimate the casualties in Turkey was studied by Tutsoy et al. [64]. In 2020, the progress of COVID-19 in Turkey was estimated by Özdinç et al. [65]. Three mathematical models for forecasting the COVID-19 outbreak in Iran and Turkey were assessed by Niazkar et al. in 2020 [66]. The forecasting epidemic size for Turkey and Iraq using the logistic model was made by Ahmed et al. in 2020 [67]. Atangana and Araz studied the mathematical model of COVID-19 spread in Turkey and South Africa in 2020 [68]. Djilali and Ghanbari estimated analysis of the peak outbreak epidemic in South Africa, Turkey, and Brazil in 2020 [69]. The dynamics of the outbreak in Hubei and Turkey were predicted and analyzed by Aslan et al. in 2020 [70]. Atangana and Araz modeled third waves of COVID-19 spread with piecewise differential and integral operators for Turkey, Spain, and Czechia in 2021 [71].

On the other hand, various numerical methods based on the Pell–Lucas polynomials were studied to obtain the approximate solutions of some differential equations and integro-differential equations [7,72,73,74,75,76,77]. Accordingly, it is concluded that effective results are obtained with the help of the Pell–Lucas polynomials. To date, there is still no the collocation method based on the Pell–Lucas polynomials among the studies on the approximate solutions of the SIR model problem. Therefore, in this study, the parameters of the SIR model problem are determined according to Covid-19 data in Turkey and the Pell–Lucas collocation method is applied to this model.

In this study, the SIR epidemic model is considered in [47,52,57]

$$\begin{array}{c}\frac{dS\left(t\right)}{dt}=-\frac{\beta }{P}S\left(t\right)I\left(t\right),\hfill \\ \frac{dI\left(t\right)}{dt}=\frac{\beta }{P}S\left(t\right)I\left(t\right)-\gamma I\left(t\right),0\le t\le b,\hfill \\ \frac{dR\left(t\right)}{dt}=\gamma I\left(t\right),\hfill \end{array}$$

(1)

with the initial conditions

$$\begin{array}{c}S\left(0\right)=S_0=P-I_0-R_0,\hfill \\ I\left(0\right)=I_0,\hfill \\ R\left(0\right)=R_0,\hfill \end{array}$$

(2)

where $P=S\left(t\right)+I\left(t\right)+R\left(t\right)$. That is, population size P is constant.

The descriptions of the parameters and the variables in the model (1) and (2) are given in

Table 1. Representations of the parameters and the variables in the model (1) and (2).

Parameter/Variable	Explanation
t	The independent variable in units of days
$S\left(t\right)$	The dependent variable showing the number of the susceptible individuals at time t
$I\left(t\right)$	The dependent variable showing the number of individuals infected with COVID-19 at time t
$R\left(t\right)$	The dependent variable showing the number of individuals removed (recovered and died) from COVID-19 at time t
$\beta$	The rate of contact or transmission
$\gamma$	The rate of recovery

. Additionally, the arrows in Figure 1 indicate the flow between the populations of susceptible $\left(S\right)$, infected $\left(I\right)$, removed $\left(R\right)$. Note that the individuals $R\left(t\right)$ in the model represents the number of individuals who both recovered and died.

Our aim is to find the Pell–Lucas polynomial solutions of the model (1) and (2) as follows:

$$\begin{array}{c}{S}_N\left(t\right)=\sum _{n=0}^N{a}_n{Q}_n\left(t\right),\hfill \\ I_N\left(t\right)=\sum _{n=0}^N{b}_n{Q}_n\left(t\right),\hfill \\ R_N\left(t\right)=\sum _{n=0}^N{c}_n{Q}_n\left(t\right),\hfill \end{array}$$

(3)

where N is any positive integer, and $a_n$, $b_n$, $c_n$, $d_n$ are the Pell–Lucas coefficients. In addition, $Q_n\left(t\right)$$(n=0,\dots ,N)$ are the Pell–Lucas polynomials defined by [78,79]

$$Q_n\left(t\right)=\sum _{k=0}^{〚n/2〛}{2}^{n-2k}\frac{n}{n-k}\left(\genfrac{}{}{0pt}{}{n-k}{k}\right)t^{n-2k}.$$

(4)

Here, $〚n/2〛$ shows the integer value of $n/2$. For features about the Pell–Lucas polynomials, please see [78,79].

2. Fundamental Matrix Relations

In this section, the Pell–Lucas polynomial solutions of the SIR model (1) and (2) are written in matrix forms.

Lemma 1.

The Pell–Lucas polynomials $Q_n\left(t\right)$ in (4) are expressed in the following matrix form [77],

$${\mathbf{Q}}_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N,$$

(5)

where ${\mathbf{T}}_N\left(t\right)=\left[\begin{array}{ccccc}1& t& t^2& \cdots & t^N\end{array}\right]$, ${\mathbf{Q}}_N\left(t\right)=\left[\begin{array}{ccccc}{Q}_0\left(t\right)& Q_1\left(t\right)& Q_2\left(t\right)& \cdots & Q_N\left(t\right)\end{array}\right]$ and if N is even

$${\mathbf{D}}_N^T=\left[\begin{array}{ccccc}2& 0& 0& \cdots & 0\\ 0& 2^1\frac{1}{1}\left(\genfrac{}{}{0pt}{}{1}{0}\right)& 0& \cdots & 0\\ 2^0\frac{2}{1}\left(\genfrac{}{}{0pt}{}{1}{1}\right)& 0& 2^2\frac{2}{2}\left(\genfrac{}{}{0pt}{}{2}{0}\right)& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 2^0\frac{N}{\frac{N}{2}}\left(\genfrac{}{}{0pt}{}{\frac{N}{2}}{\frac{N}{2}}\right)& 0& 2^2\frac{N}{\frac{N+2}{2}}\left(\genfrac{}{}{0pt}{}{\frac{N+2}{2}}{\frac{N-2}{2}}\right)& \cdots & 2^N\frac{N}{N}\left(\genfrac{}{}{0pt}{}{N}{0}\right)\end{array}\right],$$

and if N is odd

$${\mathbf{D}}_N^T=\left[\begin{array}{ccccc}2& 0& 0& \cdots & 0\\ 0& 2^1\frac{1}{1}\left(\genfrac{}{}{0pt}{}{1}{0}\right)& 0& \cdots & 0\\ 2^0\frac{2}{1}\left(\genfrac{}{}{0pt}{}{1}{1}\right)& 0& 2^2\frac{2}{2}\left(\genfrac{}{}{0pt}{}{2}{0}\right)& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 2^1\frac{N}{\frac{N+1}{2}}\left(\genfrac{}{}{0pt}{}{\frac{N+1}{2}}{\frac{N-1}{2}}\right)& 0& \cdots & 2^N\frac{N}{N}\left(\genfrac{}{}{0pt}{}{N}{0}\right)\end{array}\right].$$

Proof. 

When the vector ${\mathbf{T}}_N\left(t\right)$ is multiplied by the matrix ${\mathbf{D}}_N$ from the right side, we have the vector ${\mathbf{T}}_N\left(t\right){\mathbf{D}}_N$, which is ${\mathbf{Q}}_N\left(t\right)$. □

Lemma 2.

The Pell–Lucas polynomial solutions (3) of the SIR model (1) and (2) for any selected value of N are written in following forms:

$$\begin{array}{c}S\left(t\right)\approxeq S_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N,\hfill \\ I\left(t\right)\approxeq I_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ R\left(t\right)\approxeq R_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{C}}_N,\hfill \end{array}$$

(6)

where

$${\mathbf{A}}_N={\left[\begin{array}{cccc}{a}_0\hfill & a_1\hfill & \cdots \hfill & a_N\hfill \end{array}\right]}^T,{\mathbf{B}}_N={\left[\begin{array}{cccc}{b}_0\hfill & b_1\hfill & \cdots \hfill & b_N\hfill \end{array}\right]}^T,{\mathbf{C}}_N={\left[\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \cdots \hfill & c_N\hfill \end{array}\right]}^T.$$

Here, the matrices ${\mathbf{T}}_N\left(t\right)$ and ${\mathbf{D}}_N$ are as in Lemma 1.

Proof. 

If the vector ${\mathbf{T}}_N\left(t\right){\mathbf{D}}_N$ is multiplied by ${\mathbf{A}}_N$ from the right, we get $S_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N$. Similarly, the vector ${\mathbf{T}}_N\left(t\right){\mathbf{D}}_N$ is multiplied by ${\mathbf{B}}_N$ from the right, we have $I_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N$. Finally, when the vector ${\mathbf{T}}_N\left(t\right){\mathbf{D}}_N$ is multiplied from the right by ${\mathbf{C}}_N$, the approximate solution $R_N\left(t\right)$ is obtained in matrix form as $R_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{C}}_N$. □

Lemma 3.

The matrix relations for the derivatives of the Pell–Lucas polynomial solutions (3) are as follows:

$$\begin{array}{c}{S}^{{}^{\prime }}\left(t\right)\approxeq S_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{A}}_N,\hfill \\ I^{{}^{\prime }}\left(t\right)\approxeq I_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ R^{{}^{\prime }}\left(t\right)\approxeq R_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{C}}_N,\hfill \end{array}$$

(7)

where

$${\mathbf{H}}_N=\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 2& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & N\\ 0& 0& 0& \cdots & 0\end{array}\right].$$

Here, the matrices ${\mathbf{T}}_N\left(t\right)$, ${\mathbf{D}}_N$, ${\mathbf{A}}_N$, ${\mathbf{B}}_N$ and ${\mathbf{C}}_N$ are as in Lemma 2.

Proof. 

By taking the derivatives of the solutions in matrix forms (6), the following matrix forms are obtained:

$$\begin{array}{c}{S}^{{}^{\prime }}\left(t\right)\approxeq S_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N^{{}^{\prime }}\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N,\hfill \\ I^{{}^{\prime }}\left(t\right)\approxeq I_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N^{{}^{\prime }}\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ R^{{}^{\prime }}\left(t\right)\approxeq R_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N^{{}^{\prime }}\left(t\right){\mathbf{D}}_N{\mathbf{C}}_N.\hfill \end{array}$$

(8)

Now, the derivative of the matrix ${\mathbf{T}}_N\left(t\right)$ is taken and so the term ${\mathbf{T}}_N^{{}^{\prime }}\left(t\right)$ is converted to the form [77]

$${\mathbf{T}}_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{H}}_N.$$

(9)

Hence, the relation (9) is substituted in (8) and then the approximate solutions are written in the next forms

$$\begin{array}{c}{S}^{{}^{\prime }}\left(t\right)\approxeq S_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{A}}_N,\hfill \\ I^{{}^{\prime }}\left(t\right)\approxeq I_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ R^{{}^{\prime }}\left(t\right)\approxeq R_N^{{}^{\prime }}\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{C}}_N.\hfill \end{array}$$

□

Lemma 4.

The matrix representation of the nonlinear term in the SIR model (1) for any selected value of N is written as

$$S\left(t\right)I\left(t\right)\approxeq S_N\left(t\right)I_N\left(t\right)=\left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N\right).$$

(10)

Here, the matrices ${\mathbf{T}}_N\left(t\right)$, ${\mathbf{D}}_N$, ${\mathbf{A}}_N$ and ${\mathbf{B}}_N$ are as in Lemma 2.

Proof. 

If we use the matrix representations of $S\left(t\right)\approxeq S_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N$ and $I\left(t\right)\approxeq I_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N$ in the Lemma 2, then we have

$$S\left(t\right)I\left(t\right)\approxeq S_N\left(t\right)I_N\left(t\right)=\left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N\right).$$

(11)

□

Lemma 5.

The matrix relations of the initial conditions (2) for the solutions (3) are in forms

$$\left\{\begin{array}{c}{\mathbf{U}}_N{\mathbf{A}}_N=S_0,{\mathbf{U}}_N={\mathbf{T}}_N\left(0\right){\mathbf{D}}_N,\hfill \\ {\mathbf{U}}_N{\mathbf{B}}_N=I_0,{\mathbf{U}}_N={\mathbf{T}}_N\left(0\right){\mathbf{D}}_N,\hfill \\ {\mathbf{U}}_N{\mathbf{C}}_N=R_0,{\mathbf{U}}_N={\mathbf{T}}_N\left(0\right){\mathbf{D}}_N.\hfill \end{array}\right.$$

(12)

Here, the matrices ${\mathbf{T}}_N\left(t\right)$, ${\mathbf{D}}_N$, ${\mathbf{A}}_N$, ${\mathbf{B}}_N$ and ${\mathbf{C}}_N$ are as in Lemma 2.

Proof. 

By writting 0 instead of t in the equations in the system (6), we obtain the following matrix relations:

$$\begin{array}{c}S\left(0\right)\approxeq S_0\left(0\right)={\mathbf{T}}_N\left(0\right){\mathbf{D}}_N{\mathbf{A}}_N,\hfill \\ I\left(0\right)\approxeq I_0\left(0\right)={\mathbf{T}}_N\left(0\right){\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ R\left(0\right)\approxeq R_0\left(0\right)={\mathbf{T}}_N\left(0\right){\mathbf{D}}_N{\mathbf{C}}_N.\hfill \end{array}$$

(13)

Consequently, the matrix multiplication ${\mathbf{T}}_N\left(0\right){\mathbf{D}}_N$ is represented by ${\mathbf{U}}_N$, and thus we have the matrix relations in the Equation (12). □

Theorem 1.

It is supposed that the solutions of the model (1) and (2) are sought in the form (3). In that case, we get the following matrix relations:

$$\left\{\begin{array}{c}{\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{A}}_N=-\frac{\beta }{P}\left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N\right),\hfill \\ {\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{B}}_N=\frac{\beta }{P}\left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N\right)-\gamma {\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ {\mathbf{T}}_N\left(t\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{C}}_N=\gamma \left({\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N\right).\hfill \end{array}\right.$$

(14)

Here, the matrices ${\mathbf{T}}_N\left(t\right)$, ${\mathbf{H}}_N$, ${\mathbf{D}}_N$, ${\mathbf{A}}_N$, ${\mathbf{B}}_N$ and ${\mathbf{C}}_N$ are as in Lemmas 2 and 3.

Proof. 

If Lemma 3 is used for the terms $S^{{}^{\prime }}\left(t\right)$, $I^{{}^{\prime }}\left(t\right)$ and $R^{{}^{\prime }}\left(t\right)$ in Equation (1), Lemma 2 is used for the term $I\left(t\right)$ in Equation (1), and Lemma 4 is used for the term $S\left(t\right)I\left(t\right)$ in Equation (1), the proof is completed. □

3. The Method for the Solutions of the SIR Model

In this section, a collocation method based on the Pell–Lucas polynomials is presented for the SIR model. In application of the method, we use the evenly spaced collocation points.

Definition 1.

The evenly spaced collocation points in $[0,b]$ are defined by

$$t_i=\frac{b}{N}i,i=0,1,\dots ,N.$$

(15)

Theorem 2.

It is assumed that the approximate solutions of the Equation (1) under the conditions (2) can be represented in the form (3). In that case, the model (1) can be reduced to the system

$$\left\{\begin{array}{c}{\mathbf{W}}_0{\mathbf{A}}_N+{\mathbf{G}}_{1,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_0{\mathbf{B}}_N+{\mathbf{G}}_{2,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_0{\mathbf{C}}_N+{\mathbf{G}}_{3,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{A}}_N+{\mathbf{G}}_{1,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{B}}_N+{\mathbf{G}}_{2,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{C}}_N+{\mathbf{G}}_{3,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ ⋮\\ {\mathbf{W}}_N{\mathbf{A}}_N+{\mathbf{G}}_{1,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_N{\mathbf{B}}_N+{\mathbf{G}}_{2,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_N{\mathbf{C}}_N+{\mathbf{G}}_{3,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\end{array}\right.$$

(16)

where

$$\begin{array}{c}{\mathbf{W}}_i={\mathbf{T}}_N\left(t_i\right){\mathbf{H}}_N{\mathbf{D}}_N,\hfill \\ {\mathbf{G}}_{1,i}=\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N{\mathbf{A}}_N\right){\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N,\hfill \\ {\mathbf{G}}_{2,i}=-\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N{\mathbf{A}}_N\right){\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N+\gamma {\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N,\hfill \\ {\mathbf{G}}_{3,i}=-\gamma {\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N,\hfill \\ {\mathbf{0}}_{(N+1)\times 1}=zeros((N+1)\times 1).\hfill \end{array}$$

Here, the matrices ${\mathbf{A}}_N$, ${\mathbf{B}}_N$, ${\mathbf{C}}_N$, ${\mathbf{T}}_N\left(t_i\right)$, ${\mathbf{H}}_N$ and ${\mathbf{D}}_N$ are as in Lemmas 2 and 3.

Proof. 

By writing the collocation points (15) in the Equation (14), we get

$$\left\{\begin{array}{c}{\mathbf{T}}_N\left(t_0\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{A}}_N=-\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_0\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t_0\right){\mathbf{D}}_N{\mathbf{B}}_N\right),\hfill \\ {\mathbf{T}}_N\left(t_0\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{B}}_N=\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_0\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t_0\right){\mathbf{D}}_N{\mathbf{B}}_N\right)-\gamma {\mathbf{T}}_N\left(t_0\right){\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ {\mathbf{T}}_N\left(t_0\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{C}}_N=\gamma \left({\mathbf{T}}_N\left(t_0\right){\mathbf{D}}_N{\mathbf{B}}_N\right),\hfill \\ {\mathbf{T}}_N\left(t_1\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{A}}_N=-\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_1\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t_1\right){\mathbf{D}}_N{\mathbf{B}}_N\right),\hfill \\ {\mathbf{T}}_N\left(t_1\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{B}}_N=\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_1\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t_1\right){\mathbf{D}}_N{\mathbf{B}}_N\right)-\gamma {\mathbf{T}}_N\left(t_1\right){\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ {\mathbf{T}}_N\left(t_1\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{C}}_N=\gamma \left({\mathbf{T}}_N\left(t_1\right){\mathbf{D}}_N{\mathbf{B}}_N\right),\hfill \\ ⋮\hfill \\ {\mathbf{T}}_N\left(t_N\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{A}}_N=-\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_N\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t_N\right){\mathbf{D}}_N{\mathbf{B}}_N\right),\hfill \\ {\mathbf{T}}_N\left(t_N\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{B}}_N=\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_N\right){\mathbf{D}}_N{\mathbf{A}}_N\right)\left({\mathbf{T}}_N\left(t_N\right){\mathbf{D}}_N{\mathbf{B}}_N\right)-\gamma {\mathbf{T}}_N\left(t_N\right){\mathbf{D}}_N{\mathbf{B}}_N,\hfill \\ {\mathbf{T}}_N\left(t_N\right){\mathbf{H}}_N{\mathbf{D}}_N{\mathbf{C}}_N=\gamma \left({\mathbf{T}}_N\left(t_N\right){\mathbf{D}}_N{\mathbf{B}}_N\right),\hfill \end{array}\right.$$

(17)

or

$$\left\{\begin{array}{c}{\mathbf{W}}_0{\mathbf{A}}_N+{\mathbf{G}}_{1,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_0{\mathbf{B}}_N+{\mathbf{G}}_{2,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_0{\mathbf{C}}_N+{\mathbf{G}}_{3,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{A}}_N+{\mathbf{G}}_{1,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{B}}_N+{\mathbf{G}}_{2,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{C}}_N+{\mathbf{G}}_{3,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ ⋮\\ {\mathbf{W}}_N{\mathbf{A}}_N+{\mathbf{G}}_{1,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_N{\mathbf{B}}_N+{\mathbf{G}}_{2,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_N{\mathbf{C}}_N+{\mathbf{G}}_{3,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1}.\end{array}\right.$$

(18)

Consequently, by using the following equations,

$$\begin{array}{c}{\mathbf{W}}_i={\mathbf{T}}_N\left(t_i\right){\mathbf{H}}_N{\mathbf{D}}_N,\hfill \\ {\mathbf{G}}_{1,i}=\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N{\mathbf{A}}_N\right){\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N,\hfill \\ {\mathbf{G}}_{2,i}=-\frac{\beta }{P}\left({\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N{\mathbf{A}}_N\right){\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N+\gamma {\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N,\hfill \\ {\mathbf{G}}_{3,i}=-\gamma {\mathbf{T}}_N\left(t_i\right){\mathbf{D}}_N,\hfill \end{array}$$

we complete the proof. □

Theorem 3.

Supposing that the Pell–Lucas polynomial solutions of the SIR model (1) and (2) can be represented as in (3). Then, we have the following system:

$$\left\{\begin{array}{c}{\mathbf{W}}_0{\mathbf{A}}_N+{\mathbf{G}}_{1,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_0{\mathbf{B}}_N+{\mathbf{G}}_{2,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_0{\mathbf{C}}_N+{\mathbf{G}}_{3,0}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{A}}_N+{\mathbf{G}}_{1,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{B}}_N+{\mathbf{G}}_{2,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_1{\mathbf{C}}_N+{\mathbf{G}}_{3,1}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ ⋮\\ {\mathbf{W}}_N{\mathbf{A}}_N+{\mathbf{G}}_{1,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_N{\mathbf{B}}_N+{\mathbf{G}}_{2,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{W}}_N{\mathbf{C}}_N+{\mathbf{G}}_{3,N}{\mathbf{B}}_N={\mathbf{0}}_{(N+1)\times 1},\\ {\mathbf{U}}_N{\mathbf{A}}_N=S_0,\\ {\mathbf{U}}_N{\mathbf{B}}_N=I_0,\\ {\mathbf{U}}_N{\mathbf{C}}_N=R_0.\end{array}\right.$$

(19)

Here, the matrices ${\mathbf{W}}_i$, ${\mathbf{G}}_{i,N}$, ${\mathbf{A}}_N$, ${\mathbf{B}}_N$, ${\mathbf{C}}_N$, ${\mathbf{0}}_{(N+1)\times 1}$ and ${\mathbf{U}}_N$ are as in Theorem 2 and Lemma 5.

Proof. 

If the matrix systems (12) and (16) are written as a single system, then we obtain a new $3(N+2)\times 1$-dimensional matrix system. Hence, we get the desired result. □

Corollary 1.

By solving the obtained system (19) with the help of a program written in MATLAB, we have the coefficient matrices ${\mathbf{A}}_N$, ${\mathbf{B}}_N$ and ${\mathbf{C}}_N$ in (6). The calculated coefficient matrices ${\mathbf{A}}_N$, ${\mathbf{B}}_N$ and ${\mathbf{C}}_N$ are written in (6) and thus the approximate solutions of the model (1) and (2) are found.

4. Error Analysis

In this section, we give two important theorems. First, we determine the upper boundary of the errors for the method. Secondly, we present an error estimation method by using the residual function.

Theorem 4.

(Upper Boundary of Errors) Let $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$ be the exact solutions of the problem (1) and (2). It is supposed that $S_N\left(t\right)$, $I_N\left(t\right)$, $R_N\left(t\right)$ are the Pell–Lucas polynomial solutions (3) with $N-th$ degree of the problem (1) and (2). In addition, the expansions of the generalized Maclaurin series with $N-th$ degree of $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$ are $S_N^M\left(t\right)$, $I_N^M\left(t\right)$, $R_N^M\left(t\right)$. Then, the absolute errors of the Pell–Lucas polynomial solutions for $0\le t\le b$ are bounded by the inequality

$$\begin{array}{c}\parallel S\left(t\right)-S_N{\left(t\right)\parallel }_{\infty }\le k_N(\parallel {\tilde{\mathbf{A}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }\parallel {\mathbf{A}}_N{\parallel }_{\infty })+\frac{b^{N+1}}{(N+1)!}{\parallel S^{(N+1)}\left(c_t\right)\parallel }_{\infty },\hfill \\ \parallel I\left(t\right)-I_N{\left(t\right)\parallel }_{\infty }\le k_N(\parallel {\tilde{\mathbf{B}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }\parallel {\mathbf{B}}_N{\parallel }_{\infty })+\frac{b^{N+1}}{(N+1)!}{\parallel I^{(N+1)}\left(c_t\right)\parallel }_{\infty },\hfill \\ \parallel R\left(t\right)-R_N{\left(t\right)\parallel }_{\infty }\le k_N(\parallel {\tilde{\mathbf{C}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }\parallel {\mathbf{C}}_N{\parallel }_{\infty })+\frac{b^{N+1}}{(N+1)!}{\parallel R^{(N+1)}\left(c_t\right)\parallel }_{\infty },\hfill \end{array}$$

(20)

where $\parallel {\mathbf{T}}_N{\left(t\right)\parallel }_{\infty }\le max\{b^N,1\}:=k_N$, $▵{\mathbf{A}}_N=\parallel {\mathbf{A}}_{N+1}{\parallel }_{\infty }-{\parallel {\mathbf{A}}_N\parallel }_{\infty }$, $▵{\mathbf{B}}_N=\parallel {\mathbf{B}}_{N+1}{\parallel }_{\infty }-{\parallel {\mathbf{B}}_N\parallel }_{\infty }$, $▵{\mathbf{C}}_N=\parallel {\mathbf{C}}_{N+1}{\parallel }_{\infty }-{\parallel {\mathbf{C}}_N\parallel }_{\infty }$. Also, the coefficient matrix of $S_N^M\left(t\right)$, the coeffcient matrix of $I_N^M\left(t\right)$, the coeffcient matrix of $R_N^M\left(t\right)$ are represented, respectively, ${\tilde{\mathbf{A}}}_N$, ${\tilde{\mathbf{B}}}_N$ and ${\tilde{\mathbf{C}}}_N$.

Proof. 

First, we add and subtract the functions $S\left(t\right)-S_N\left(t\right)$, $I\left(t\right)-I_N\left(t\right)$, $R\left(t\right)-R_N\left(t\right)$ to the functions the Maclaurin expansions $S_N^M\left(t\right)$, $I_N^M\left(t\right)$, $R_N^M\left(t\right)$, respectively. Next, we use the triangle inequality and so we have

$$\begin{array}{c}\parallel S\left(t\right)-S_N{\left(t\right)\parallel }_{\infty }=\parallel S\left(t\right)-S_N^M\left(t\right)+S_N^M\left(t\right)-S_N{\left(t\right)\parallel }_{\infty }\le \parallel S\left(t\right)-S_N^M{\left(t\right)\parallel }_{\infty }+{\parallel S_N^M\left(t\right)-S_N\left(t\right)\parallel }_{\infty },\hfill \\ \parallel I\left(t\right)-I_N{\left(t\right)\parallel }_{\infty }=\parallel I\left(t\right)-I_N^M\left(t\right)+I_N^M\left(t\right)-I_N{\left(t\right)\parallel }_{\infty }\le \parallel I\left(t\right)-I_N^M{\left(t\right)\parallel }_{\infty }+{\parallel I_N^M\left(t\right)-I_N\left(t\right)\parallel }_{\infty },\hfill \\ \parallel R\left(t\right)-R_N{\left(t\right)\parallel }_{\infty }=\parallel R\left(t\right)-R_N^M\left(t\right)+R_N^M\left(t\right)-R_N{\left(t\right)\parallel }_{\infty }\le \parallel R\left(t\right)-R_N^M{\left(t\right)\parallel }_{\infty }+{\parallel R_N^M\left(t\right)-R_N\left(t\right)\parallel }_{\infty }.\hfill \end{array}$$

(21)

By examining the terms $\parallel S\left(t\right)-S_N^M{\left(t\right)\parallel }_{\infty }$, $\parallel I\left(t\right)-I_N^M{\left(t\right)\parallel }_{\infty }$, $\parallel R\left(t\right)-R_N^M{\left(t\right)\parallel }_{\infty }$, we write the remainder terms of the Maclaurin series $S_N^M\left(t\right)$, $I_N^M\left(t\right)$, $R_N^M\left(t\right)$ as follows:

$$\begin{array}{c}\frac{t^{N+1}}{(N+1)!}{S}^{(N+1)}\left(c_t\right),0\le t\le b,\hfill \\ \frac{t^{N+1}}{(N+1)!}{I}^{(N+1)}\left(c_t\right),0\le t\le b,\hfill \\ \frac{t^{N+1}}{(N+1)!}{R}^{(N+1)}\left(c_t\right),0\le t\le b,\hfill \end{array}$$

(22)

and thus we get

$$\begin{array}{c}\parallel S\left(t\right)-S_N^M{\left(t\right)\parallel }_{\infty }\le \frac{b^{N+1}}{(N+1)!}{\parallel S^{(N+1)}\left(c_t\right)\parallel }_{\infty },\hfill \\ \parallel I\left(t\right)-I_N^M{\left(t\right)\parallel }_{\infty }\le \frac{b^{N+1}}{(N+1)!}{\parallel I^{(N+1)}\left(c_t\right)\parallel }_{\infty },\hfill \\ \parallel R\left(t\right)-R_N^M{\left(t\right)\parallel }_{\infty }\le \frac{b^{N+1}}{(N+1)!}{\parallel R^{(N+1)}\left(c_t\right)\parallel }_{\infty }.\hfill \end{array}$$

(23)

According to Lemma 2, we know that $S_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{A}}_N$, $I_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{B}}_N$, $R_N\left(t\right)={\mathbf{T}}_N\left(t\right){\mathbf{D}}_N{\mathbf{C}}_N$ are the matrix forms of the Pell–Lucas polynomial solutions $S_N\left(t\right)$, $I_N\left(t\right)$, $R_N\left(t\right)$, respectively. In addition, we denote the expansions of the Maclaurin series of $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$ as $S_N^M\left(t\right)={\mathbf{T}}_N\left(t\right){\tilde{\mathbf{A}}}_N$, $I_N^M\left(t\right)={\mathbf{T}}_N\left(t\right){\tilde{\mathbf{B}}}_N$, $R_N^M\left(t\right)={\mathbf{T}}_N\left(t\right){\tilde{\mathbf{C}}}_N$. Hence, we can write the terms $\parallel S_N^M\left(t\right)-S_N{\left(t\right)\parallel }_{\infty }$, $\parallel I_N^M\left(t\right)-I_N{\left(t\right)\parallel }_{\infty }$, $\parallel R_N^M\left(t\right)-R_N{\left(t\right)\parallel }_{\infty }$ as the following forms:

$$\begin{array}{c}\parallel S_N^M\left(t\right)-S_N{\left(t\right)\parallel }_{\infty }=\parallel {\mathbf{T}}_N\left(t\right)({\tilde{\mathbf{A}}}_N-{\mathbf{D}}_N{\mathbf{A}}_N){\parallel }_{\infty }\le {\parallel {\mathbf{T}}_N\left(t\right)\parallel }_{\infty }\left(\parallel {\tilde{\mathbf{A}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }{\parallel {\mathbf{A}}_N\parallel }_{\infty }\right),\hfill \\ \parallel I_N^M\left(t\right)-I_N{\left(t\right)\parallel }_{\infty }=\parallel {\mathbf{T}}_N\left(t\right)({\tilde{\mathbf{B}}}_N-{\mathbf{D}}_N{\mathbf{B}}_N){\parallel }_{\infty }\le {\parallel {\mathbf{T}}_N\left(t\right)\parallel }_{\infty }\left(\parallel {\tilde{\mathbf{B}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }{\parallel {\mathbf{B}}_N\parallel }_{\infty }\right),\hfill \\ \parallel R_N^M\left(t\right)-R_N{\left(t\right)\parallel }_{\infty }=\parallel {\mathbf{T}}_N\left(t\right)({\tilde{\mathbf{C}}}_N-{\mathbf{D}}_N{\mathbf{C}}_N){\parallel }_{\infty }\le {\parallel {\mathbf{T}}_N\left(t\right)\parallel }_{\infty }\left(\parallel {\tilde{\mathbf{C}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }{\parallel {\mathbf{C}}_N\parallel }_{\infty }\right).\hfill \end{array}$$

(24)

Now, because of $0\le t\le b$, we express the term $\parallel {\mathbf{T}}_N{\left(t\right)\parallel }_{\infty }$ as follows:

$$\parallel {\mathbf{T}}_N{\left(t\right)\parallel }_{\infty }\le max\{b^N,1\}:=k_N.$$

(25)

By using the expression (25), we get the inequalities in (24) as

$$\begin{array}{c}\parallel S_N^M\left(t\right)-S_N{\left(t\right)\parallel }_{\infty }\le k_N\left(\parallel {\tilde{\mathbf{A}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }{\parallel {\mathbf{A}}_N\parallel }_{\infty }\right),\hfill \\ \parallel I_N^M\left(t\right)-I_N{\left(t\right)\parallel }_{\infty }\le k_N\left(\parallel {\tilde{\mathbf{B}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }{\parallel {\mathbf{B}}_N\parallel }_{\infty }\right),\hfill \\ \parallel R_N^M\left(t\right)-R_N{\left(t\right)\parallel }_{\infty }\le k_N\left(\parallel {\tilde{\mathbf{C}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }{\parallel {\mathbf{C}}_N\parallel }_{\infty }\right).\hfill \end{array}$$

(26)

By substituting the inequalities (23) and (26) in (21), we have

$$\begin{array}{c}\parallel S\left(t\right)-S_N{\left(t\right)\parallel }_{\infty }\le k_N(\parallel {\tilde{\mathbf{A}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }\parallel {\mathbf{A}}_N{\parallel }_{\infty })+\frac{b^{N+1}}{(N+1)!}{\parallel S^{(N+1)}\left(c_t\right)\parallel }_{\infty },\hfill \\ \parallel I\left(t\right)-I_N{\left(t\right)\parallel }_{\infty }\le k_N(\parallel {\tilde{\mathbf{B}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }\parallel {\mathbf{B}}_N{\parallel }_{\infty })+\frac{b^{N+1}}{(N+1)!}{\parallel I^{(N+1)}\left(c_t\right)\parallel }_{\infty },\hfill \\ \parallel R\left(t\right)-R_N{\left(t\right)\parallel }_{\infty }\le k_N(\parallel {\tilde{\mathbf{C}}}_N{\parallel }_{\infty }+\parallel {\mathbf{D}}_N{\parallel }_{\infty }\parallel {\mathbf{C}}_N{\parallel }_{\infty })+\frac{b^{N+1}}{(N+1)!}{\parallel R^{(N+1)}\left(c_t\right)\parallel }_{\infty }.\hfill \end{array}$$

(27)

As a result, the proof is completed. □

Theorem 5.

(Error Estimation) Let $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$ be the exact solutions of the model (1) and (2) and $S_N\left(t\right)$, $I_N\left(t\right)$, $R_N\left(t\right)$ be the Pell–Lucas polynomial solutions (3) with $N-th$ degree of the model (1) and (2). In this case, the following error problem is obtained:

$$\left\{\begin{array}{c}{e}_{S,N}^{{}^{\prime }}\left(t\right)+\frac{\beta }{P}\left(e_{S,N}\left(t\right)e_{I,N}\left(t\right)+I_N\left(t\right)e_{S,N}\left(t\right)+S_N\left(t\right)e_{I,N}\left(t\right)\right)=-R{e}_{1,N}\left(t\right),\hfill \\ e_{I,N}^{{}^{\prime }}\left(t\right)-\frac{\beta }{P}\left(e_{S,N}\left(t\right)e_{I,N}\left(t\right)+I_N\left(t\right)e_{S,N}\left(t\right)+S_N\left(t\right)e_{I,N}\left(t\right)\right)+\gamma e_{I,N}\left(t\right)=-R{e}_{2,N}\left(t\right),\hfill \\ e_{R,N}^{{}^{\prime }}\left(t\right)-\gamma e_{I,N}\left(t\right)=-R{e}_{3,N}\left(t\right),\hfill \\ e_{S,N}\left(0\right)=0,e_{I,N}\left(0\right)=0,e_{R,N}\left(0\right)=0.\hfill \end{array}\right.$$

(28)

Here, $e_{S,N}\left(t\right)=S\left(t\right)-S_N\left(t\right)$, $e_{I,N}\left(t\right)=I\left(t\right)-I_N\left(t\right)$, $e_{R,N}\left(t\right)=R\left(t\right)-R_N\left(t\right)$. In addition, $R{e}_{1,N}\left(t\right)$, $R{e}_{2,N}\left(t\right)$, $R{e}_{3,N}\left(t\right)$ are the residual functions of the model (1) and (2) for the the Pell–Lucas polynomial solutions (3).

Proof. 

Because the Pell–Lucas polynomial solutions in (3) provide the Equation (1) and initial conditions (2), we can write

$$\begin{array}{c}R{e}_{1,N}\left(t\right)=S_N^{{}^{\prime }}\left(t\right)+\frac{\beta }{P}{S}_N\left(t\right)I_N\left(t\right),\hfill \\ R{e}_{2,N}\left(t\right)=I_N^{{}^{\prime }}\left(t\right)-\frac{\beta }{P}{S}_N\left(t\right)I_N\left(t\right)+\gamma I_N\left(t\right),\hfill \\ R{e}_{3,N}\left(t\right)=R_N^{{}^{\prime }}\left(t\right)-\gamma I_N\left(t\right),\hfill \\ S_N\left(0\right)=S_0,I_N\left(0\right)=I_0,R_N\left(0\right)=R_0.\hfill \end{array}$$

(29)

The model (29) is subtracted from model (1) and (2) and thus we have the error problem

$$\left\{\begin{array}{c}{e}_{S,N}^{{}^{\prime }}\left(t\right)+\frac{\beta }{P}\left(e_{S,N}\left(t\right)e_{I,N}\left(t\right)+I_N\left(t\right)e_{S,N}\left(t\right)+S_N\left(t\right)e_{I,N}\left(t\right)\right)=-R{e}_{1,N}\left(t\right),\hfill \\ e_{I,N}^{{}^{\prime }}\left(t\right)-\frac{\beta }{P}\left(e_{S,N}\left(t\right)e_{I,N}\left(t\right)+I_N\left(t\right)e_{S,N}\left(t\right)+S_N\left(t\right)e_{I,N}\left(t\right)\right)+\gamma e_{I,N}\left(t\right)=-R{e}_{2,N}\left(t\right),\hfill \\ e_{R,N}^{{}^{\prime }}\left(t\right)-\gamma e_{I,N}\left(t\right)=-R{e}_{3,N}\left(t\right),\hfill \\ e_{S,N}\left(0\right)=0,e_{I,N}\left(0\right)=0,e_{R,N}\left(0\right)=0.\hfill \end{array}\right.$$

(30)

Here, $e_{S,N}\left(t\right)=S\left(t\right)-S_N\left(t\right)$, $e_{I,N}\left(t\right)=I\left(t\right)-I_N\left(t\right)$, $e_{R,N}\left(t\right)=R\left(t\right)-R_N\left(t\right)$. Consequently, we complete the proof of the theorem. □

Corollary 2.

By solving the problem (28) with the help of the method in the previous section, we obtain the estimated error functions $e_{S,N,M}\left(t\right)$, $e_{I,N,M}\left(t\right)$, $e_{R,N,M}\left(t\right)$.

5. Numerical Verification and Discussion

In this section, we make the applications of the methods presented in the Section 3 and Section 4 for the SIR model. First, we determine the parameters and the initial conditions in this model by using the COVID-19 data in Turkey [80]. Secondly, by using a program for the method in MATLAB, we get the Pell–Lucas polynomial solutions. In addition, we compare our approximate solutions with the approximate solutions of the Runge–Kutta method. Finally, we present application results in tables and graphs and discuss the numerical verification.

In order to determine the parameters $\beta ,\gamma$ and the initial conditions $S_0,I_0,R_0$ in the SIR model (1) and (2), the COVID-19 data in Turkey are used. Hence, the numbers of the susceptible individuals, the infected individuals, the removed individuals on April 4, 2020 are selected as the initial condition [80]. In addition, we give representations of the solutions and the errors in

Table 2. Representations of the solutions and the errors in the Section 5.

Data	Explanation
$S\left(t\right)$	The susceptible individuals at time t
$I\left(t\right)$	The individuals infected with COVID-19 at time t
$R\left(t\right)$	The individuals removed (recovered and died) from COVID-19 at time t
$S_N\left(t\right)$	The susceptible individuals at time t according to the method in Section 3
$I_N\left(t\right)$	The individuals infected with COVID-19 at time t according to the method in Section 3
$R_N\left(t\right)$	The individuals removed (recovered and died) from COVID-19 at time t according to the method in Section 3
$e_{S,N,M}\left(t\right)$	The estimated error function for the susceptible population according to the method in Section 4
$e_{I,N,M}\left(t\right)$	The estimated error function for the infected population according to the method in Section 4
$e_{R,N,M}\left(t\right)$	The estimated error function for the removed population (recovered and died) according to the method in Section 4

and we give the values of parameters $\beta ,\gamma$ and initial conditions $S_0,I_0,R_0$ in SIR model (1) and (2) in

Table 3. The parameters $\beta ,\gamma$ and the initial conditions $S_0,I_0,R_0$ in the SIR model (1) and (2).

Parameters	${\mathit{S}}_0$	${\mathit{I}}_0$	${\mathit{R}}_0$	$\mathit{\beta }$	$\mathit{\gamma }$
Values	83,996,609	3013	378	$1/14$	1287/23,934
				[1/day]	[Total Removed/Total Infected]
	[80]	[80]	[80]	Estimated [81,82]	Estimated [80,83]

.

We consider the SIR epidemic model together with the conditions according to the selected parameters for Covid-19 data in Turkey as follows:

$$\begin{array}{c}\frac{dS\left(t\right)}{dt}=-8.5034e-10S\left(t\right)I\left(t\right),\hfill \\ \frac{dI\left(t\right)}{dt}=8.5034e-10S\left(t\right)I\left(t\right)-0.0538I\left(t\right),\hfill \\ \frac{dR\left(t\right)}{dt}=0.0538I\left(t\right),\hfill \\ S\left(0\right)=83996609,I\left(0\right)=3013,R\left(0\right)=378.\hfill \end{array}$$

(31)

Now, let’s apply the Pell–Lucas collocation method in the range $[0,60]$. First, we write the Pell–Lucas polynomial solutions for $N=5$ as

$$\begin{array}{c}{S}_5\left(t\right)=\sum _{n=0}^5{a}_n{Q}_n\left(t\right),\hfill \\ I_5\left(t\right)=\sum _{n=0}^5{b}_n{Q}_n\left(t\right),\hfill \\ R_5\left(t\right)=\sum _{n=0}^5{c}_n{Q}_n\left(t\right).\hfill \end{array}$$

(32)

By using the Lemma 2, we express the Pell–Lucas polynomial solutions in (32) in matrix forms

$$\begin{array}{c}S\left(t\right)\approxeq S_5\left(t\right)={\mathbf{T}}_5\left(t\right){\mathbf{D}}_5{\mathbf{A}}_5,\hfill \\ I\left(t\right)\approxeq I_5\left(t\right)={\mathbf{T}}_5\left(t\right){\mathbf{D}}_5{\mathbf{B}}_5,\hfill \\ R\left(t\right)\approxeq R_5\left(t\right)={\mathbf{T}}_5\left(t\right){\mathbf{D}}_5{\mathbf{C}}_5.\hfill \end{array}$$

(33)

Here,

$${\mathbf{A}}_5={\left[\begin{array}{cccccc}{a}_0\hfill & a_1\hfill & a_2\hfill & a_3\hfill & a_4\hfill & a_5\hfill \end{array}\right]}^T,{\mathbf{B}}_5={\left[\begin{array}{cccccc}{b}_0\hfill & b_1\hfill & b_2\hfill & b_3\hfill & b_4\hfill & b_5\hfill \end{array}\right]}^T,$$

$${\mathbf{C}}_5={\left[\begin{array}{cccccc}{c}_0\hfill & c_1\hfill & c_2\hfill & c_3\hfill & c_4\hfill & c_5\hfill \end{array}\right]}^T,{\mathbf{T}}_5\left(t\right)=\left[\begin{array}{cccccc}1& t& t^2& t^3& t^4& t^5\end{array}\right],$$

$${\mathbf{D}}_5^T=\left[\begin{array}{cccccc}2& 0& 0& 0& 0& 0\\ 0& 2& 0& 0& 0& 0\\ 2& 0& 4& 0& 0& 0\\ 0& 6& 0& 8& 0& 0\\ 2& 0& 16& 0& 16& 0\\ 0& 10& 0& 40& 0& 32\end{array}\right].$$

Secondly, we determine the collocation points for the range $[0,60]$. Because $b=60$, the collocation points become $t_0=0,t_1=12,t_2=24,t_3=36,t_4=48,t_5=60.$ Thus, by using the system (16), we get

$$\left\{\begin{array}{c}{\mathbf{W}}_0{\mathbf{A}}_5+{\mathbf{G}}_{1,0}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\\ {\mathbf{W}}_0{\mathbf{B}}_5+{\mathbf{G}}_{2,0}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\\ {\mathbf{W}}_0{\mathbf{C}}_5+{\mathbf{G}}_{3,0}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\\ {\mathbf{W}}_1{\mathbf{A}}_5+{\mathbf{G}}_{1,1}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\\ {\mathbf{W}}_1{\mathbf{B}}_5+{\mathbf{G}}_{2,1}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\\ {\mathbf{W}}_1{\mathbf{C}}_5+{\mathbf{G}}_{3,1}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\\ ⋮\\ {\mathbf{W}}_5{\mathbf{A}}_5+{\mathbf{G}}_{1,5}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\\ {\mathbf{W}}_5{\mathbf{B}}_5+{\mathbf{G}}_{2,5}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\\ {\mathbf{W}}_5{\mathbf{C}}_5+{\mathbf{G}}_{3,5}{\mathbf{B}}_5={\mathbf{0}}_{6\times 1},\end{array}\right.$$

(34)

where

$$\begin{array}{c}{\mathbf{W}}_i={\mathbf{T}}_5\left(t_i\right){\mathbf{H}}_5{\mathbf{D}}_5,\hfill \\ {\mathbf{G}}_{1,i}=8.5034e-10\left({\mathbf{T}}_5\left(t_i\right){\mathbf{D}}_5{\mathbf{A}}_5\right){\mathbf{T}}_5\left(t_i\right){\mathbf{D}}_5,\hfill \\ {\mathbf{G}}_{2,i}=-8.5034e-10\left({\mathbf{T}}_5\left(t_i\right){\mathbf{D}}_5{\mathbf{A}}_5\right){\mathbf{T}}_5\left(t_i\right){\mathbf{D}}_5+0.0538{\mathbf{T}}_5\left(t_i\right){\mathbf{D}}_5,\hfill \\ {\mathbf{G}}_{3,i}=-0.0538{\mathbf{T}}_5\left(t_i\right){\mathbf{D}}_5,\hfill \\ \\ {\mathbf{D}}_5=\left[\begin{array}{cccccc}2& 0& 2& 0& 2& 0\\ 0& 2& 0& 6& 0& 10\\ 0& 0& 4& 0& 16& 0\\ 0& 0& 0& 8& 0& 40\\ 0& 0& 0& 0& 16& 0\\ 0& 0& 0& 0& 0& 32\end{array}\right],{\mathbf{H}}_5=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 2& 0& 0& 0\\ 0& 0& 0& 3& 0& 0\\ 0& 0& 0& 0& 4& 0\\ 0& 0& 0& 0& 0& 5\\ 0& 0& 0& 0& 0& 0\end{array}\right],\hfill \\ \\ {\mathbf{0}}_{6\times 1}=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\end{array}\right],{\mathbf{T}}_5\left(0\right)=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\end{array}\right],\hfill \\ \\ {\mathbf{T}}_5\left(12\right)=\left[\begin{array}{cccccc}1& 12& {12}^2& {12}^3& {12}^4& {12}^5\end{array}\right],{\mathbf{T}}_5\left(24\right)=\left[\begin{array}{cccccc}1& 24& {24}^2& {24}^3& {24}^4& {24}^5\end{array}\right],\hfill \end{array}$$

$$\begin{array}{c}{\mathbf{T}}_5\left(36\right)=\left[\begin{array}{cccccc}1& 36& {36}^2& {36}^3& {36}^4& {36}^5\end{array}\right],{\mathbf{T}}_5\left(48\right)=\left[\begin{array}{cccccc}1& 48& {48}^2& {48}^3& {48}^4& {48}^5\end{array}\right],\hfill \\ \\ {\mathbf{T}}_5\left(60\right)=\left[\begin{array}{cccccc}1& 60& {60}^2& {60}^3& {60}^4& {60}^5\end{array}\right].\hfill \end{array}$$

Subsequently, we express the matrix relations of the initial conditions by using (12) in the following matrix forms:

$$\left\{\begin{array}{c}{\mathbf{U}}_5{\mathbf{A}}_5=83996609,{\mathbf{U}}_5={\mathbf{T}}_5\left(0\right){\mathbf{D}}_5,\hfill \\ {\mathbf{U}}_5{\mathbf{B}}_5=3013,{\mathbf{U}}_5={\mathbf{T}}_5\left(0\right){\mathbf{D}}_5,\hfill \\ {\mathbf{U}}_5{\mathbf{C}}_5=378,{\mathbf{U}}_5={\mathbf{T}}_5\left(0\right){\mathbf{D}}_5.\hfill \end{array}\right.$$

(35)

Here, ${\mathbf{T}}_5\left(0\right)=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\end{array}\right].$

As the next step, we combine (34) and (35) and we solve the combined system with the help of MATLAB. The solution of this system determines the coefficients matrices ${\mathbf{A}}_5$, ${\mathbf{B}}_5$ and ${\mathbf{C}}_5$. By writing the determined coefficient matrices in (33), we get the approximate solutions of (1) and (2) as

$$\begin{array}{cc}{S}_5\left(t\right)=\hfill & 83996533-2.1537e+02t-1.8846t^2-0.0118t^3-3.4679e-05t^4\hfill \\ \hfill & -3.2610e-07t^5,\hfill \\ I_5\left(t\right)=\hfill & 3013+5.3228e+01t+0.4657t^2+0.0029t^3+8.6146e-06t^4\hfill \\ \hfill & +7.9518e-08t^5,\hfill \\ R_5\left(t\right)=\hfill & 378+1.6215e+02t+1.4190t^2+0.0089t^3+2.6064e-05t^4\hfill \\ \hfill & +2.4658e-07t^5.\hfill \end{array}$$

(36)

In Figure 2, Figure 3 and Figure 4, we show the Pell–Lucas polynomial solutions $S_N\left(t\right)$, $I_N\left(t\right)$, $R_N\left(t\right)$ of the SIR model (31) for $N=5$, $N=8$ and $N=10$. According to this, we interpret that although the susceptible population is decreasing, the infected population and the removed population are increasing. In Figure 5, we demonstrate that the Pell–Lucas polynomial solutions $I_N\left(t\right)$ and $R_N\left(t\right)$ of model (31) for $N=5$. From here, we said that the removed population is increased at a greater rate. Accordingly, the removed rate is quite high compared to the infected rate at 60 days. Also, we compare the Pell–Lucas polynomial solutions $S_N\left(t\right)$, $I_N\left(t\right)$, $R_N\left(t\right)$ of the SIR model (31) for $N=5$ with those of the Runge–Kutta method in Figure 6. According to Figure 6, it is said that the graphs of the presented method and the Runge–Kutta method are similar. That is, we observe that the method is accurate and effective.

In Figure 7, Figure 8 and Figure 9, we compare the residual absolute error functions of the SIR model (31) for $N=5$, $N=8$ and $N=10$. In addition, we compare the estimated absolute error functions of the SIR model (31) for $(N,M)=(5,6)$, $(N,M)=(8,9)$ and $(N,M)=(10,11)$ in Figure 10, Figure 11 and Figure 12. Accordingly, we observe that as the value of N increases, the errors decrease.

In

Table 4. Comparison of the residual absolute errors and the estimated absolute errors of the susceptible individuals.

	Residual Absolute Errors		Estimated Absolute Errors
${\mathit{t}}_{\mathit{i}}$	$|{\mathit{R}}_{\mathbf{1},\mathbf{8}}\left(\mathit{t}\right)|$	$|{\mathit{R}}_{\mathbf{1},\mathbf{10}}\left(\mathit{t}\right)|$	$|{\mathit{e}}_{\mathbf{1},\mathbf{8},\mathbf{9}}\left(\mathit{t}\right)|$	$|{\mathit{e}}_{\mathbf{1},\mathbf{10},\mathbf{11}}\left(\mathit{t}\right)|$
0	2.5055 × 10${}^{-5}$	1.7821 × 10${}^{-8}$	2.7755 × 10${}^{-22}$	6.9117 × 10${}^{-26}$
10	3.9551 × 10${}^{-7}$	1.0382 × 10${}^{-10}$	5.4301 × 10${}^{-5}$	3.0992 × 10${}^{-8}$
20	5.4811 × 10${}^{-8}$	1.0065 × 10${}^{-11}$	5.3574 × 10${}^{-5}$	3.1013 × 10${}^{-8}$
30	1.7356 × 10${}^{-14}$	1.7963 × 10${}^{-14}$	5.3476 × 10${}^{-5}$	3.1007 × 10${}^{-8}$
40	2.7206 × 10${}^{-8}$	5.4862 × 10${}^{-12}$	5.3568 × 10${}^{-5}$	3.1005 × 10${}^{-8}$
50	7.7954 × 10${}^{-8}$	2.3212 × 10${}^{-11}$	5.3686 × 10${}^{-5}$	3.0992 × 10${}^{-8}$
60	8.1501 × 10${}^{-15}$	6.2135 × 10${}^{-14}$	5.2078 × 10${}^{-5}$	3.03 × 10${}^{-8}$

,

Table 5. Comparison of the residual absolute errors and the estimated absolute errors of the infected individuals.

	Residual Absolute Errors		Estimated Absolute Errors
${\mathit{t}}_{\mathit{i}}$	$|{\mathit{R}}_{\mathbf{2},\mathbf{8}}\left(\mathit{t}\right)|$	$|{\mathit{R}}_{\mathbf{2},\mathbf{10}}\left(\mathit{t}\right)|$	$|{\mathit{e}}_{\mathbf{2},\mathbf{8},\mathbf{9}}\left(\mathit{t}\right)|$	$|{\mathit{e}}_{\mathbf{2},\mathbf{10},\mathbf{11}}\left(\mathit{t}\right)|$
0	4.5783 × 10${}^{-6}$	1.5858 × 10${}^{-8}$	1.2801 × 10${}^{-23}$	4.7937 × 10${}^{-26}$
10	7.1278 × 10${}^{-8}$	9.2552 × 10${}^{-11}$	5.4365 × 10${}^{-6}$	1.5551 × 10${}^{-8}$
20	9.7302 × 10${}^{-9}$	9.0502 × 10${}^{-12}$	3.107 × 10${}^{-6}$	9.1025 × 10${}^{-9}$
30	1.0093 × 10${}^{-13}$	5.9179 × 10${}^{-14}$	1.7995 × 10${}^{-6}$	5.3174 × 10${}^{-9}$
40	4.6659 × 10${}^{-9}$	4.9293 × 10${}^{-12}$	1.0621 × 10${}^{-6}$	3.1079 × 10${}^{-9}$
50	1.3106 × 10${}^{-8}$	2.1138 × 10${}^{-11}$	6.4465 × 10${}^{-7}$	1.8073 × 10${}^{-9}$
60	3.4531 × 10${}^{-14}$	3.2251 × 10${}^{-14}$	1.1109 × 10${}^{-7}$	4.2956 × 10${}^{-10}$

and

Table 6. Comparison of the residual absolute errors and the estimated absolute errors of the removed individuals.

	Residual Absolute Errors		Estimated Absolute Errors
${\mathit{t}}_{\mathit{i}}$	$|{\mathit{R}}_{\mathbf{3},\mathbf{8}}\left(\mathit{t}\right)|$	$|{\mathit{R}}_{\mathbf{3},\mathbf{10}}\left(\mathit{t}\right)|$	$|{\mathit{e}}_{\mathbf{3},\mathbf{8},\mathbf{9}}\left(\mathit{t}\right)|$	$|{\mathit{e}}_{\mathbf{3},\mathbf{10},\mathbf{11}}\left(\mathit{t}\right)|$
0	2.0477 × 10${}^{-5}$	1.9418 × 10${}^{-9}$	3.159 × 10${}^{-23}$	2.9812 × 10${}^{-26}$
10	3.2423 × 10${}^{-7}$	1.1024 × 10${}^{-11}$	4.8865 × 10${}^{-5}$	1.5403 × 10${}^{-8}$
20	4.5081 × 10${}^{-8}$	1.0284 × 10${}^{-12}$	5.0467 × 10${}^{-5}$	2.1872 × 10${}^{-8}$
30	1.4048 × 10${}^{-14}$	5.2996 × 10${}^{-14}$	5.1677 × 10${}^{-5}$	2.5652 × 10${}^{-8}$
40	2.254 × 10${}^{-8}$	5.4701 × 10${}^{-13}$	5.2506 × 10${}^{-5}$	2.786 × 10${}^{-8}$
50	6.4847 × 10${}^{-8}$	2.2316 × 10${}^{-12}$	5.3041 × 10${}^{-5}$	2.9147 × 10${}^{-8}$
60	2.9582 × 10${}^{-15}$	4.9546 × 10${}^{-14}$	5.1967 × 10${}^{-5}$	2.9831 × 10${}^{-8}$

, we tabulate the residual absolute errors and the estimated absolute errors of the SIR model (31) for $(N,M)=(8,9)$ and $(N,M)=(10,11)$. According to Table 4, Table 5 and Table 6, we observe that as the value of N increases, the error decreases. Although the residual absolute errors are better than the estimation absolute errors, the estimation absolute errors are not bad either. In other words, the error estimation method presented in Section 4 give very successful results.

6. Conclusions

This paper proposes a numerical method for an SIR model to investigate the present condition of COVID-19 disease contamination and to estimate its future improvements in Turkey. The parameters and the initial conditions of this model are determined by using real data. The presented method is a collocation approach based on the Pell–Lucas polynomials. According to the Pell–Lucas collocation method, the SIR model is reduced to a system of nonlinear algebraic equations. The solutions of this nonlinear algebraic system determine coefficients of the Pell–Lucas polynomial solutions of the SIR model. Additionally, two error analyses are made. According to Figure 2, Figure 3 and Figure 4, it is interpreted that although the susceptible population is decreasing, the infected population and the removed population are increasing. Also in Figure 5, it is observed that the removed population increases from 378 to 17,667 for the same value of N whereas the infected population increases from 3013 to 8685 for $N=5$. In the 60-day period from 4 April 2020, an increase in the number of the infected patients is observed. Nevertheless, a faster increase is observed in the number of the removed patients. In that case, we expect that the pandemic will diminish when enough isolation precautions are continued. In Figure 6, we compare the approximate solutions $S_N\left(t\right)$, $I_N\left(t\right)$, $R_N\left(t\right)$ for $N=5$ with those of the Runge–Kutta method. Accordingly, it is concluded that the graphs obtained from the presented method and the Runge–Kutta method are similar.

In Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 and Table 4, Table 5 and Table 6, we examine the residual absolute errors and the estimated absolute errors of the approximate solution functions. According to these, we deduce that as the value of N increases the error decreases. Even though, the residual absolute errors are better than the estimation absolute errors, the estimation absolute errors are not bad either. Accordingly, we comment that the Pell–Lucas collocation method is the effective method to get the approximate solutions of the SIR model. A limitation of the method is that the individuals $R\left(t\right)$ in the model represents the number of individuals who both recovered and died. However, the method can be improved by making necessary adjustments to the model. A more important advantage of the method than all these advantages is that the parameters in the model can be determined for different countries, and this method can be developed for other countries as well. Moreover, this method can be developed for similar infections. In the future, in similar epidemic situations, the method can be applied by determining the parameters of the model and the initial conditions in the model. Moreover, the results are obtained in a very short time thanks to the code written in MATLAB. Hereby, the cautious provisions can be made to minimize infections and to intercept an overloading of the health system.