A Novel Technique to Control the Accuracy of a Nonlinear Fractional Order Model of COVID-19: Application of the CESTAC Method and the CADNA Library

mathematics Article A Novel Technique to Control the Accuracy of a Nonlinear Fractional Order Model of COVID-19: Application of the CESTAC Method and the CADNA Library Samad Noeiaghdam 1,2,* , Sanda Micula 3and Juan J. Nieto 4   Citation: Noeiaghdam, S.; Micula, S.; Nieto, J.J. A Novel Technique to Control the Accuracy of a Nonlinear Fractional Order Model of COVID-19: Application of the CESTAC Method and the CADNA Library. Mathematics 2021,9, 1321. https://doi.org/ 10.3390/math9121321 Academic Editor: Mikhail Kolev Received: 1 May 2021 Accepted: 7 June 2021 Published: 8 June 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Industrial Mathematics Laboratory, Baikal School of BRICS, Irkutsk National Research Technical University, 664074 Irkutsk, Russia 2Department of Applied Mathematics and Programming, South Ural State University, Lenin Prospect 76, 454080 Chelyabinsk, Russia 3 Department of Mathematics and Computer Science, Babe¸s-Bolyai University, 400084 Cluj-Napoca, Romania; [email protected] 4 Instituto de Matemáticas, Departamento de Estatística, Análise Matemática e Optimización, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain; [email protected] *Correspondence: [email protected] or [email protected] Abstract: In this paper, a nonlinear fractional order model of COVID-19 is approximated. For this aim, at first we apply the Caputo–Fabrizio fractional derivative to model the usual form of the phenomenon. In order to show the existence of a solution, the Banach fixed point theorem and the Picard–Lindelof approach are used. Additionally, the stability analysis is discussed using the fixed point theorem. The model is approximated based on Indian data and using the homotopy analysis transform method (HATM), which is among the most famous, flexible and applicable semi-analytical methods. After that, the CESTAC (Controle et Estimation Stochastique des Arrondis de Calculs) method and the CADNA (Control of Accuracy and Debugging for Numerical Applications) library, which are based on discrete stochastic arithmetic (DSA), are applied to validate the numerical results of the HATM. Additionally, the stopping condition in the numerical algorithm is based on two successive approximations and the main theorem of the CESTAC method can aid us analytically to apply the new terminations criterion instead of the usual absolute error that we use in the floatingpoint arithmetic (FPA). Finding the optimal approximations and the optimal iteration of the HATM to solve the nonlinear fractional order model of COVID-19 are the main novelties of this study. Keywords: COVID-19; Caputo–Fabrizio derivative; homotopy analysis transform method; discrete stochastic arithmetic; CESTAC method; CADNA library 1. Introduction Corona virus infection (COVID-19) is an infectious disease caused by the newly discovered Corona virus. Most people with COVID-19 disease experience mild to moderate symptoms and recover without special treatment. The virus that causes COVID-19 is mainly transmitted through particles produced when a person coughs, sneezes, or exhales. These particles do not stay suspended in the air due to their weight and fall quickly on the ground or surfaces. If you are close to a person with COVID-19, you may become infected by inhaling the virus or touching the infected surface and then touching your eyes, nose, or mouth. According to WHO reports until 2 March 2021, more than 115 million infected, more than 2.5 million dead and more than 90 million recovered people from more than 200 countries have been identified [1]. Thus, providing mathematical models to control and predict the behavior of this pandemic is urgent [ 2 ]. Recently, many mathematical models related to COVID-19 have been studied [ 3 , 4 ]. In [ 5 ], an analysis model of diagnosis and treatment for the COVID- 19 pandemic based on medical information fusion was discussed, in [ 6 ] a model of the Mathematics 2021,9, 1321. https://doi.org/10.3390/math9121321 https://www.mdpi.com/journal/mathematics Mathematics 2021,9, 1321 2 of 26 COVID-19 epidemic and predicting its future situation in Ethiopia was studied and in [ 7 ] a fractional order COVID-19 model with real-life data application was analyzed. The COVID- 19 epidemic and clinical risk factors of patients under an epidemiological Markov model were discussed in [ 8 ] and the incorporating social determinants of health into modeling of COVID-19 and other infectious diseases was illustrated in [9]. Recently, due to their importance and more accurate results, fractional order models are among the most studied problems in the world. In [ 10 ], the authors focused on a Caputo–Fabrizio fractional model of COVID-19. In [ 11 ], a fractional order model of HIV infection was discussed, in [ 12 ] the application of the Atangana–Baleanu fractional derivative in mathematical biology was studied and in [ 13 ], a fractional order HIV model with Caputo and constant proportional Caputo operators was illustrated. An application of the Caputo–Fabrizio fractional derivative to solve a fractional order model of an energy supply–demand system can be found in [ 14 ]. For other applications of the fractional order models and also their usual forms please see [15–20]. The homotopy analysis method (HAM) presented by Liao [ 21 – 23 ] is one of the most flexible and accurate methods for solving various linear and nonlinear problems. This method has been applied for solving many problems, such as solving integral equations [24,25] , solving singular problems [ 26 , 27 ] and a modified non-linear epidemiological model of computer viruses [28]. We should note that all of the mentioned researches are based on the FPA, and in order to show the accuracy of the methods, we apply the traditional absolute error or residual error in the form of |g(x)−gm(x)|<ε,or |gm(x)−gm−1(x)|<ε, (1) where the exact and approximate solutions are denoted by g(x) and gm(x) , respectively and ε is a small positive value. Thus, in order to apply Condition (1), not only should we know the exact solution g(x) but we should also have the optimal value of ε , which are disadvantages of methods based on the FPA. Concerning ε , for small values we will have extra iterations without improving the accuracy of the results and for large values we will have only one or two iterations without providing accurate results. Thus, we apply the CESTAC method and the CADNA library, which are based on the DSA. This method was presented by Laporte and Vignes [ 29 ] for the first time and after that other researchers improved the method and library [ 30 ]. Using this method we do not need to have the exact solutions and εand the termination criterion will be in the following form: |gm(x)−gm−1(x)|[email protected], (2) where gm(x) and gm−1(x) are two successive approximations and @.0 denotes the informatical zero. Applying Condition (2), we do not need to have the exact solution and it is important that it does not depend on a positive small value such as ε . Additionally, in this method, instead of applying the usual softwares, we apply the CADNA library. This library should be applied on a Linux operating system and all codes should be written using C, C++, FORTRAN or ADA codes [ 31 – 33 ]. The informatical zero sign can be produced only using the DSA and the CADNA library. This shows that the number of common significant digits for two successive approximations are almost equal to the number of common significant digits of the exact and approximate solutions [ 33 – 35 ]. Thus, we will be able to apply Condition (2) instead of (1). Recently, we focused on validating the numerical results of various problems using the CESTAC method. In [ 36 ], the results of the Adomian decomposition method for solving Volterra integral equation with a discontinuous kernel was studied and in [ 37 , 38 ] the CESTAC method was used to validate the results on the reverse osmosis system. Dynamical control on the homotopy perturbation method and the Taylor-collocation method to solve Volterra integral equations with piecewise smooth kernels was discussed in [ 39 , 40 ]. In [ 41 ], we illustrated the Sinc-collocation method for Mathematics 2021,9, 1321 3 of 26 solving various kinds of crisp and fuzzy integral equations and finally, we applied the CADNA library to control the accuracy of the load leveling problem in [42]. In this paper, a nonlinear fractional order model of COVID-19 will be discussed based on the Caputo–Fabrizio fractional derivative. The Picard–Lindelof approach and the Banach fixed point theorem are used to show the existence of the solution of the model. Additionally, the stability of the model is illustrated applying the fixed point theorem. The HATM is applied to solve the model and the numerical results are validated applying the CESTAC method and the CADNA library. Additionally, the main theorem of the CESTAC method is proven, which guarantees the accuracy of the numerical results obtained using Condition (2) instead of (1). Using this method, the optimal approximations and the optimal iteration to solve the model using the HATM are obtained, which are the main novelties of this study. 2. Preliminaries In this section, we present some definitions and details of fractional derivatives [ 43 , 44 ]. Definition 1 ([ 43 ]) . Let g∈Cn . The ω -th order Caputo fractional derivative can be defined via integrable differentiations as follows: CDωg(t) = 1 Γ(n−ω)Zt 0 g(n)(s) (t−s)ω−n+1ds,n= [ω] + 1, where g(n)∈L1. Definition 2 ([ 43 ]) . The ω -th order Caputo–Fabrizio derivative for b>a , g∈H1(a , b) , and ω∈(0, 1)is defined in the following form: CFDωg(t) = M(ω) 1−ωZt aexp−ω 1−ω(t−s)g0(s)ds,t≥0, where M(ω) depending on ω denotes a normalization function and M( 0 ) = M( 1 ) = 1. For a function g that is not in H1(a,b), this derivative can be presented for g ∈L1(˘∞,b)as follows: CFDωg(t) = ωM(ω) 1−ωZb −∞(g(t)−g(s))exp(−ω 1−ω(t−s))ds, 0 <ω<1. In addition, the ω+n -th order fractional derivative CFDω+n for n≥ 1and ω∈( 0, 1 ) can be defined as CFDω+ng(t):=CF Dω(Dng(t)). Definition 3 ([ 44 ]) . For 0 <ω≤ 1and M(ω) = 1, the Laplace transform for the Caputo– Fabrizio derivative can be presented in the following form: L[CFDω+ng(t)](s) = sn+1L[g(t)] −sng(0)−sn−1g0(0)−. . . −f(n)(0) s+ω(1−s). Definition 4 ([ 44 ]) . The ω -order Riemann–Liouville fractional integral for Re(ω)> 0is specified as: Iωg(t) = 1 Γ(ω)Zt 0(t−s)ω−1g(s)ds. Definition 5 ([44]).The fractional integral of Caputo–Fabrizio is defined by CF Iωg(t) = 2(1−ω) (2−ω)M(ω)g(t) + 2(ω) (2−ω)M(ω)Zt 0g(s)ds, 0 <ω<1. Mathematics 2021,9, 1321 4 of 26 Additionally, the left and right fractional integrals of ( CF aDω ) are defined [ 45 ] respectively by: (CF Iω ag)(t) = 1−ω B(ω)g(t) + ω B(ω)Zt ag(s)ds, (CF Iω bg)(t) = 1−ω B(ω)g(t) + ω B(ω)Zb tg(s)ds. We have the following definitions for the Sumudu transform, x which is based on the classical Fourier integral [46–48]. Definition 6 ([ 47 ]) . Let A={G:∃λ , k1 , k2≥ 0, |G(t)|<λexp(t kj) , t∈( ˘1 )j×[ 0, ∞)} . Then ST[g(t) ; u] = G(u) denotes the Sumudu transform of a function g(t)∈A that is defined as G(u) = ST[g(t);u] = 1 uZ∞ 0exp(−t/u)g(t)dt,u∈(−k1,k2), for all t ≥0. Additionally, g(t) = ST˘1[G(u)] is the inverse Sumudu transform of G(u). Definition 7 ([ 46 ]) . For (m ˘1 <ω≤m) we can define the Sumudu transform for the Caputo derivative as ST[cDω tg(t);u] = u−ω"G(u)− m ∑ i=0 uω−i[cDω−ig(t)]t=0#. Definition 8 ([ 49 ]) . Assume that G is a function and it has a Caputo–Fabrizio fractional derivative. We can define the Sumudu transform of G with the Caputo–Fabrizio fractional derivative as ST(CF 0Dω t)(G(t)) = M(ω) 1−ω+ωu[ST(G(t)) −G(0)]. Definition 9. [ 50 ] For a metric space (X , d) , a map g:X→X is called a Picard operator whenever there exists x∗∈X such that Fix(g) = {x∗} and the sequence (gn(x0))n∈N converges to x∗for all x0∈X. Definition 10. Consider the Banach space (Y , k . k) , a self-map G on Y , and the recursive relation Pn+1=φ(G , Pn) . Let Ω(G) be the fixed point set of G , which includes Ω(G)6=∅ and limn→∞Pn=p∈Ω(G) . Suppose that {gn} ⊂ Ω and en=gn+1 ˘ φ(G , gn) , if limn→∞en= 0 implies that limn→∞gn=p , then the recursive procedure Pn+1=φ(G , Pn) is G -stable. Suppose that our sequence gn has an upper boundary. If Picard’s iteration Pn+1=GPn is satisfied in all these conditions, then Pn+1=GPnis G-stable. Theorem 1 ([50]).Assume that (Y,k.k)is a Banach space and G is a self-map of Y such that kGx−Gyk ≤ Rkx−Gxk+rkx−yk, for all x,y∈Y, where R ≥0and 0≤r<1. Then G is Picard G-stable. Mathematics 2021,9, 1321 5 of 26 3. Nonlinear Fractional Order Model of COVID-19 Consider the following non-linear bio-mathematical model [51,52]:                                                      dS dt =−βIS, dE dt =βIS −εE, dI dt =εE−γI−dI −qI, dQ dt =qI −qtQ−dQ, dR dt =γI+qtQ, dD dt =dI +dQ, (3) where S , E , I , Q , R and D are susceptible, exposed, infected, quarantined, recovered and dead populations, respectively, and we denote by N=S+E+I+Q+R+D the total size of the population. Susceptible individuals become infected at a rate, β , by infectious individuals, and the rate can be obtained using β=R0γ where R0 is the basic reproduction number and γ=1 in f ectious period . ε=1 incubation period . We should note that D is a fraction of I+Q . The death rate is denoted by d , the fraction of active cases quarantined is denoted by q and the time period of quarantine is indicated by qt . The graphical form of the model is shown in Figure 1[ 51 ]. The values and parameters of this model are presented in Table 1 , based on Indian data [51]. Model (3) can be written using the ω -th order Caputo–Fabrizio fractional derivative as follows:                                      CF 0Dω tS=−βIS, CF 0Dω tE=βIS −εE, CF 0Dω tI=εE−γI−dI −qI, CF 0Dω tQ=qI −qtQ−dQ, CF 0Dω tR=γI+qtQ, CF 0Dω tD=dI +dQ. (4) Figure 1. The model of COVID-19 [51]. Mathematics 2021,9, 1321 6 of 26 Table 1. List of parameters and functions [51]. Parameters Values Incubation Period 1 εMean= 5.1 days Infectious Period 1 γMean= 7 days Basic Reproduction Number R0Mean= 2.28 Mean Death Rate 0.43% Active Cases Quarantined q0.01 Period Quarantined qt14 days Quarantined Cases S020 Recovered Cases R012 Deaths D05 Exposed Cases E020 Infected Cases I015 Quarantined Cases Q015 4. Existence of Solution Applying Definition 5to both sides of Equation (4), we obtain:                                                            S(t)−u1(t) = 2(1−ω) (2−ω)M(ω){−βI(t)S(t)}+2ω (2−ω)M(ω)Zt 0[−βI(s)S(s)]ds, E(t)−u2(t) = 2(1−ω) (2−ω)M(ω){βI(t)S(t)−εE(t)}+2ω (2−ω)M(ω)Zt 0[βI(s)S(s)−εE(s)]ds, I(t)−u3(t) = 2(1−ω) (2−ω)M(ω){εE(t)−γI(t)−dI(t)−qI(t)}+2ω (2−ω)M(ω)Zt 0[εE(s)−γI(s)−dI(s)−qI(s)]ds, Q(t)−u4(t) = 2(1−ω) (2−ω)M(ω){qI(t)−qtQ(t)−dQ(t)}+2ω (2−ω)M(ω)Zt 0[qI(s)−qtQ(s)−dQ(s)]ds, R(t)−u5(t) = 2(1−ω) (2−ω)M(ω){γI(t) + qtQ(t)}+2ω (2−ω)M(ω)Zt 0[γI(s) + qtQ(s)]ds, D(t)−u6(t) = 2(1−ω) (2−ω)M(ω){dI(t) + dQ(t)}+2ω (2−ω)M(ω)Zt 0[dI(s) + dQ(s)]ds. (5) By assuming S0(t) = u1(t) , E0(t) = u2(t) , I0(t) = u3(t) , Q0(t) = u4(t) , R0(t) = u5(t) and D0(t) = u6(t), these equations can be written in the following form:                                                            Sn+1(t) = 2(1−ω) (2−ω)M(ω){−βIn(t)Sn(t)}+2ω (2−ω)M(ω)Zt 0[−βIn(s)Sn(s)]ds, En+1(t) = 2(1−ω) (2−ω)M(ω){βIn(t)Sn(t)−εEn(t)}+2ω (2−ω)M(ω)Zt 0[βIn(s)Sn(s)−εEn(s)]ds, In+1(t) = 2(1−ω) (2−ω)M(ω){εEn(t)−γIn(t)−dIn(t)−qIn(t)}+2ω (2−ω)M(ω)Zt 0[εEn(s)−γIn(s)−dIn(s)−qIn(s)]ds, Qn+1(t) = 2(1−ω) (2−ω)M(ω){qIn(t)−qtQn(t)−dQn(t)}+2ω (2−ω)M(ω)Zt 0[qIn(s)−qtQn(s)−dQn(s)]ds, Rn+1(t) = 2(1−ω) (2−ω)M(ω){γIn(t) + qtQn(t)}+2ω (2−ω)M(ω)Zt 0[γIn(s) + qtQn(s)]ds, Dn+1(t) = 2(1−ω) (2−ω)M(ω){dIn(t) + dQn(t)}+2ω (2−ω)M(ω)Zt 0[dIn(s) + dQn(s)]ds, (6) and by letting n approach infinity in Equation (6), the approximate solutions can be obtained as: limn→∞Sn(t) = S(t), limn→∞Qn(t) = Q(t), limn→∞En(t) = E(t), limn→∞Rn(t) = R(t), limn→∞In(t) = I(t), limn→∞Dn(t) = D(t). (7) Mathematics 2021,9, 1321 7 of 26 Now we can apply the Picard–Lindelof approach and the Banach fixed point theorem to show the existence of a solution to Equation (4). First we define the following operators as:                                      u1(t,S) = −βIS, u2(t,E) = βIS −εE, u3(t,I) = εE−γI−dI −qI, u4(t,Q) = qI −qtQ−dQ, u5(t,R) = γI+qtQ, u6(t,D) = dI +dQ. (8) Let `1=sup C[a,ρ1]ku1(t,S)k,`2=sup C[a,ρ2]ku2(t,E)k,`3=sup C[a,ρ3]ku3(t,I)k, `4=sup C[a,ρ4]ku4(t,Q)k,`5=sup C[a,ρ5]ku5(t,R)k,`6=sup C[a,ρ6]ku6(t,D)k, (9) where                                      C[a,ρ1] = |t−a,t+a|×|S−ρ1,S+ρ1|=B×C1, C[a,ρ2] = |t−a,t+a|×|E−ρ2,E+ρ2|=B×C2, C[a,ρ3] = |t−a,t+a|×|I−ρ3,I+ρ3|=B×C3, C[a,ρ4] = |t−a,t+a|×|Q−ρ4,Q+ρ4|=B×C4, C[a,ρ5] = |t−a,t+a|×|R−ρ5,R+ρ5|=B×C5, C[a,ρ6] = |t−a,t+a|×|D−ρ6,D+ρ6|=B×C6. (10) Assuming a uniform norm kZ(t)k∞=sup t∈[t−a,t+a]|Z(t)|, (11) on C[a,ρi],i=1, . . . , 5, and defining the Picard operator as O:C(B,C1, . . . , C6)→C(B,C1, . . . , C6), (12) we can write O(Z(t)) = Z0(t) + 2(1−ω) (2−ω)M(ω)G(t,Z(t)) + 2ω (2−ω)M(ω)Zt 0G(s,Z(s))ds, (13) where Z(t) = {S(t) , E(t) , I(t) , Q(t) , R(t) , D(t)} and Z0(t) = {S( 0 ) , E( 0 ) , I( 0 ) , Q( 0 ) , R( 0 ) , D(0)}and G(t,Z(t)) = {u1(t,S),u2(t,E),u3(t,I),u4(t,Q),u5(t,R),u6(t,D)}, (14) and thus kZ(t)k∞≤max{ρ1, . . . , ρ6}=C. (15) Mathematics 2021,9, 1321 8 of 26 Assume that L=max{L1, . . . , L6}and there is t0such that t0≥tand kO(Z(t)) −Z0(t)k=k2(1−ω) (2−ω)M(ω)G(t,Z(t)) + 2ω (2−v)M(v)Zt 0G(s,Z(s))dsk ≤2(1−ω) (2−ω)M(ω)kG(t,Z(t))k+2ω (2−ω)M(ω)Zt 0kG(s,Z(s))kds ≤2(1−ω) (2−ω)M(ω)+2ωt (2−ω)M(ω)L (16) ≤2(1−ω) (2−ω)M(ω)+2ωt0 (2−v)M(v)L≤µL≤C, (17) where µ<C L. (18) Now, we can write kOZ1−OZ2k=sup t∈B|Z1(t)−Z2(t)|, (19) and using definition of the Picard operator we obtain kOZ1−OZ2k=k2(1−ω) (2−ω)M(ω){G(t,Z1(t)) −G(t,Z2(t))} +2ω (2−ω)M(ω)Zt 0{G(s,Z1(s)) −G(s,Z2(s))}dsk ≤2(1−ω) (2−ω)M(ω)kG(t,Z1(t)) −G(t,Z2(t))k +2ω (2−ω)M(ω)Zt 0kG(s,Z1(s)) −G(s,Z2(s))kds ≤2(1−ω) (2−ω)M(ω)λkZ1(t)−Z2(t)k +2ωλ (2−ω)M(ω)Zt 0kZ1(s)−Z2(s)kds ≤2(1−ω)λ (2−ω)M(ω)+2ωλt0 (2−ω)M(ω)kZ1(t)−Z2(t)k ≤µλkZ1(t)−Z2(t)k, (20) where λ< 1. We know that G is a contraction so µλ < 1 and O is a contraction, showing that the proof is complete. Mathematics 2021,9, 1321 9 of 26 5. Special Solution via Iteration Approach In this section, using the Sumudu transform, we can provide a special solution. Applying this transformation on both sides of (4) we obtain:                                      ST(CF 0Dω tS(t)) = ST(−βI(t)S(t)), ST(CF 0Dω tE(t)) = ST(βI(t)S(t)−εE(t)), ST(CF 0Dω tI(t)) = ST(εE(t)−γI(t)−dI(t)−qI(t)), ST(CF 0Dω tQ(t)) = ST(qI(t)−qtQ(t)−dQ(t)), ST(CF 0Dω tR(t)) = ST(γI(t) + qtQ(t)), ST(CF 0Dω tD(t)) = ST(dI(t) + dQ(t)), (21) and based on the definition of the Sumudu transform we can write                                                        M(ω) 1−ω+ωu(ST(S(t)) −S(0)) = ST(−βI(t)S(t)), M(ω) 1−ω+ωu(ST(E(t)) −E(0)) = ST(βI(t)S(t)−εE(t)), M(ω) 1−ω+ωu(ST(I(t)) −I(0)) = ST(εE(t)−γI(t)−dI(t)−qI(t)), M(ω) 1−ω+ωu(ST(Q(t)) −Q(0)) = ST(qI(t)−qtQ(t)−dQ(t)), M(ω) 1−ω+ωu(ST(R(t)) −R(0)) = ST(γI(t) + qtQ(t)), M(ω) 1−ω+ωu(ST(D(t)) −D(0)) = ST(dI(t) + dQ(t)), (22) and                                                          ST(S(t)) = S(0) + 1−ω+ωu M(ω)ST(−βI(t)S(t)), ST(E(t)) = E(0) + 1−ω+ωu M(ω)ST(βI(t)S(t)−εE(t)), ST(I(t)) = I(0) + 1−ω+ωu M(ω)ST(εE(t)−γI(t)−dI(t)−qI(t)), ST(Q(t)) = Q(0) + 1−ω+ωu M(ω)ST(qI(t)−qtQ(t)−dQ(t)), ST(R(t)) = R(0) + 1−ω+ωu M(ω)ST(γI(t) + qtQ(t)), ST(D(t)) = D(0) + 1−ω+ωu M(ω)ST(dI(t) + dQ(t)). (23) Mathematics 2021,9, 1321 16 of 26 where                                                                                                            R1,m(~ Sm−1(t),~ Em−1(t),~ Im−1(t),~ Qm−1(t),~ Rm−1(t),~ Dm−1(t)) =L(S)−S0 s−s+ω(1−s) sL[−βIm−1Sm−1], R2,m(~ Sm−1(t),~ Em−1(t),~ Im−1(t),~ Qm−1(t),~ Rm−1(t),~ Dm−1(t)) =L(E)−E0 s−s+ω(1−s) sL[βIm−1Sm−1−εEm−1], R3,m(~ Sm−1(t),~ Em−1(t),~ Im−1(t),~ Qm−1(t),~ Rm−1(t),~ Dm−1(t)) =L(I)−I0 s−s+ω(1−s) sL[εEm−1−γIm−1−dIm−1−qIm−1], R4,m(~ Sm−1(t),~ Em−1(t),~ Im−1(t),~ Qm−1(t),~ Rm−1(t),~ Dm−1(t)) =L(Q)−Q0 s−s+ω(1−s) sL[qIm−1−qtQm−1−dQm−1], R5,m(~ Sm−1(t),~ Em−1(t),~ Im−1(t),~ Qm−1(t),~ Rm−1(t),~ Dm−1(t)) =L(R)−R0 s−s+ω(1−s) sL[γIm−1+qtQm−1], R6,m(~ Sm−1(t),~ Em−1(t),~ Im−1(t),~ Qm−1(t),~ Rm−1(t),~ Dm−1(t)) =L(D)−D0 s−s+ω(1−s) sL[dIm−1+dQm−1], (46) and using the inverse Laplace transforms we can write                                      Sm(t) = χmSm−1(t) + ¯hHL−1[R1,m(Sm−1)], Em(t) = χmEm−1(t) + ¯hHL−1[R2,m(Em−1)], Im(t) = χmIm−1(t) + ¯hHL−1[R3,m(Im−1)], Qm(t) = χmQm−1(t) + ¯hHL−1[R4,m(Qm−1)], Rm(t) = χmRm−1(t) + ¯hHL−1[R5,m(Rm−1)], Dm(t) = χmDm−1(t) + ¯hHL−1[R6,m(Dm−1)], (47) where we will apply this relation to find the successive iterations of the HATM. 7. CESTAC Method with CADNA Library Because of some advantages of the DSA in comparison to the FPA, we apply the mathematical methods based on the DSA instead of the methods based on the FPA. Thus, the CESTAC method and the CADNA library should be applied to validate the numerical results [35,53]. By collecting all representable values that are produced by the computer in set B , we can write V∗∈Bfor v∗∈Rwith αmantissa bits of the binary FPA as V∗=v∗−ψ2E−αη, (48) Mathematics 2021,9, 1321 17 of 26 where the sign, missing segment of the mantissa and the binary exponent of the result are denoted by ψ , 2 −αη and E , respectively [ 31 – 34 ]. In order to find the results with single and double precisions, the value α can be changed to 24 and 53, respectively. Let η be a random variable uniformly distributed on [− 1, 1 ] , constructing perturbation on the last mantissa bit of v∗ ; the mean (µ) and the standard deviation (σ) values can be produced for results of V∗ . If we repeat the process r times, then r samples of V∗ can be found as Φ=V∗ 1,V∗ 2, . . . , V∗ r and we should make perturbation on the bit of mantissa. For finding σ2 we need to calculate the mean value as ˜ V∗=∑r i=1V∗ i r . Thus we will have equality between the mean value and the exact v∗ . After that, the number of common significant digits between V∗and ˜ V∗, can be found by C˜ V∗,V∗=log10 √r˜ V∗ τδσ, where τδ is the value of T distribution as the confidence interval is 1 −δ , with r− 1 degrees of freedom [29,35,53]. This process will be stopped if ˜ V∗=0, or C˜ V∗,V∗≤0. When we want to apply the CESTAC method, we do not need to use the method directly. For using the method we should apply the CADNA library. This library can implement the algorithm automatically. The CADNA library should be used on a Linux operating system and all codes should be written using C, C++, FORTRAN or ADA codes. Thus, in this method we do not need to apply the usual mathematical softwares such as Mathematica, Maple and MATLAB. Applying the CESTAC method and DSA we have some advantages in comparison with the methods based on the FPA. In order to apply the termination criterion (1), which is based on the FPA, we need to have the exact solution, but in the DSA we do not need the exact solution and stopping Condition (2) is based on two successive approximations. In the FPA, we do not know the optimal ε and in the DSA we do not have the value ε definitely. In the FPA, the extra iterations can be produced without improving the accuracy, but in the DSA we can find the optimal number of iterations. In the FPA, the algorithm can be stopped in the first step without producing the accurate results but in the DSA, the optimal approximation can be identified. In the CESTAC method, we can produce @.0, which shows that the number of common significant digits between two successive approximations are zero, but in the FPA we cannot produce this sign. Definition 11 ([ 32 ]) . For two real numbers ˜ a1 and ˜ a2 , the number of common significant digits can be defined as        C˜ a1,˜ a2=log10 ˜ a1+˜ a2 2(˜ a1−˜ a2) =log10 ˜ a1 ˜ a1−˜ a2−1 2 ,˜ a16=˜ a2, C˜ a1,˜ a1= +∞. (49) Theorem 3. Let Sm(t) = ∑m j=0Sj(t),Em(t) = ∑m j=0Ej(t),Im(t) = ∑m j=0Ij(t), Qm(t) = ∑m j=0Qj(t),Rm(t) = ∑m j=0Rj(t),Dm(t) = ∑m j=0Dj(t), be the approximate solution of the mathematical model of COVID-19 (3) that is produced by the HATM; then CSm,Sm+1=CSm,S+O(1 m),CQm,Qm+1=CQm,Q+O(1 m), CEm,Em+1=CEm,E+O(1 m),CRm,Rm+1=CRm,R+O(1 m), CIm,Im+1=CIm,I+O(1 m),CDm,Dm+1=CDm,D+O(1 m). (50) Mathematics 2021,9, 1321 18 of 26 Proof. Using the mentioned definition we can write CSm,Sm+1−CSm,S=log10 Sm+Sm+1 2(Sm−Sm+1)−log10 Sm+S 2(Sm−S) =log10 Sm+Sm+1 2(S¯h m+1)−log10 Sm+S 2(Sm−S) =log10 Sm+Sm+1 Sm+S +log10 Sm−S S¯h m+1 . (51) It is obvious that in the first term of Equation (51), by increasing the number of iterations m , the approximate and exact solutions are almost equal. Thus we can neglect that term and we have log10 Sm−S Sm+1 =log10 ∑m i=0Si−∑∞ i=0Si Sm+1 =log10 ∑∞ i=m+1Si Sm+1 =log10 1+∑∞ i=m+2¯hL−1H(r)<i(ψi−1) ¯hL−1H(r)<m+1(ψm) . According to the traditional HAM [21–23], we can write ∑∞ i=m+2¯hL−1H(r)<i(ψi−1) ¯hL−1H(r)<m+1(ψm)= ∑∞ i=m+2¯hL−1H(r)1 (i−1)! ∂i−1N[Ψ(r;v)] ∂vi−1v=0 ¯hL−1H(r)1 m! ∂mN[Ψ(r;v)] ∂vmv=0 = ¯hL−1H(r)1 (m+1)! ∂m+1N[Ψ(r;v)] ∂vm+1v=0 ¯hL−1H(r)1 m! ∂mN[Ψ(r;v)] ∂vmv=0+ ¯hL−1H(r)1 (m+2)! ∂m+2N[Ψ(r;v)] ∂vm+2v=0 ¯hL−1H(r)1 m! ∂mN[Ψ(r;v)] ∂vmv=0+··· =O(1 m). Therefore, log10 Sm−S Sm+1 =log10 1+O(1 m) , log10 Qm−Q Qm+1 =log10 1+O(1 m) , log10 Em−E Em+1 =log10 1+O(1 m) , log10 Rm−R Rm+1 =log10 1+O(1 m) , log10 Im−I Im+1 =log10 1+O(1 m) , log10 Dm−D Dm+1 =log10 1+O(1 m) . It is clear that O(1 m)<< 1, and thus the right-hand sides of the above relations decrease as mincreases. Mathematics 2021,9, 1321 19 of 26 8. Numerical Results In this section, the numerical results of the HATM for solving the nonlinear fractional order model (4) are presented. As known in the HATM, we have some auxiliary functions and parameters that can help us find more accurate solutions with high-speed of convergence. In this method, the obtained solutions are based on t and ¯h such that the parameter ¯h can help identify and control the convergence region. To this aim, by plotting some ¯h -curves we can find these regions. In Figure 2, the ¯h -curves are presented for t= 0.2 and m= 15. The convergence region for functions S , E , I , Q , R and D is − 1 ≤¯h≤ − 0.6,, which is the parallel part of the ¯h -curve with the axis x . Additionally, by increasing the number of iterations of the HATM, we can plot more accurate ¯h -curves. In Figure 3, the ¯h -curves are plotted for m=25 and we obtain: −0.7 ≤¯hS≤ −0.5, −0.55 ≤¯hQ≤ −0.45, −0.7 ≤¯hE≤ −0.5, −0.65 ≤¯hR≤ −0.45, −0.7 ≤¯hI≤ −0.5, −0.7 ≤¯hD≤ −0.5. In Table 2, the numerical results obtained based on the DSA and using the CESTAC method are shown for t= 0.2, ρ= 1 and ¯h=− 0.8. However, it is important to describe why we need to apply the DSA instead of the FPA. In the methods based on the FPA, generally we apply the traditional absolute error or residual error to show the accuracy of the method, but in these cases we need to know the exact solution and also in order to stop the numerical algorithm we need to apply a suitable termination criterion. Thus, instead of applying the methods based on the FPA, we propose methods based on the DSA. In this case, we apply the CESTAC method and the CADNA library to validate the numerical results. In the CESTAC method we use the difference between two successive approximations instead of the traditional absolute error. Additionally, we do not need to have positive small value such as ε and an important note is that we do not need to have the exact solution. Using this method, we can find the optimal number of iteration and the optimal approximation of the HATM for solving the mentioned model. According to Table 4, the optimal iteration is mopt =7 and the optimal approximations are: Sopt =0.732556E+001, Qopt =0.91948E+000, Eopt =0.3162149E+002, Ropt =0.2654801E+002, Iopt =0.1575674E+002, Dopt =0.501759E+001. Sign @.0 in step 7 shows that the number of common significant digits between two successive approximations is zero. Thus, we can apply the difference between two successive approximations instead of traditional absolute error. Additionally, 24 numerical instabilities are reported by the CADNA library. In Table 3, the numerical results of the FPA based on Condition (1) for t= 0.2 and ε= 10 −4 are presented. It is obvious that the stopping condition depends on the value of ε . In Table 4, the numbers of iterations for different values of ε are obtained. As we described before, the algorithm is stopped at the first steps for large values of ε . Additionally, for small values of ε we have large number of iterations. In Table 5, the numerical results are presented for t= 0.2, ρ= 0.9 and ¯h=− 0.6 Mathematics 2021,9, 1321 20 of 26 based on the CESTAC method. Using this table we can find the optimal iteration mopt = 10 and the following optimal approximations: Sopt =0.761582E+001, Qopt =0.1916541E+002, Eopt =0.3076886E+002, Ropt =0.4576194E+002, Iopt =0.144270E+002, Dopt =0.528554E+001. In Figure 4, the approximate solution of the model using the HATM is demonstrated. We can see that by decreasing the number of susceptible people (red line), the quarantined people (dotted line) decreases, the infected people (dashed line) slowly decreases, the rate of exposed people (orange line) increases and after that slowly decreases, the death rate (green line) slowly increases, and the number of recovered people increases and then slowly increases. Thus, we have acceptable and logical results for the model. -1.0 -0.8 -0.6 -0.4 h 7.32546 7.32548 7.32550 7.32552 7.32554 7.32556 7.32558 S(t) -1.0 -0.8 -0.6 -0.4 h 31.6215 31.6216 31.6216 31.6217 E(t) -1.0 -0.8 -0.6 -0.4 h 15.7567 15.7567 15.7567 15.7567 15.7568 I(t) -1.0 -0.8 -0.6 -0.4 h 0.919 0.920 0.921 0.922 Q(t) -1.0 -0.8 -0.6 -0.4 h 26.5450 26.5455 26.5460 26.5465 26.5470 26.5475 26.5480 26.5485 R(t) -1.0 -0.8 -0.6 -0.4 h 5.01759 5.01759 5.01759 5.01759 5.01759 D(t) Figure 2. The ¯h-curves of the approximate solutions obtained form HATM for t=0.2 and m=15. Mathematics 2021,9, 1321 21 of 26 -1.0 -0.8 -0.6 -0.4 h 7.32551 7.32551 7.32551 7.32551 7.32551 7.32551 S(t) -1.0 -0.8 -0.6 -0.4 h 31.6215 31.6215 31.6215 31.6215 31.6215 E(t) -1.0 -0.8 -0.6 -0.4 h 15.7568 15.7568 15.7568 15.7568 15.7568 I(t) -1.0 -0.8 -0.6 -0.4 h 0.92174 0.92176 0.92178 0.92180 0.92182 0.92184 Q(t) -1.0 -0.8 -0.6 -0.4 h 26.5452 26.5453 26.5453 R(t) -1.0 -0.8 -0.6 -0.4 h 5.01759 5.01759 5.01759 5.01759 5.01759 5.01759 5.01759 D(t) Figure 3. The ¯h-curves of the approximate solutions obtained form HATM for t=0.2 and m=25. 0.2 0.4 0.6 0.8 1.0 t 5 10 15 20 25 30 35 D15 R15 Q15 I15 E15 S15 Figure 4. The approximate solution for m=15, ¯h=−0.8. Mathematics 2021,9, 1321 22 of 26 Table 2. Numerical results based on the DSA and using the CESTAC method for t= 0.2, ρ= 1 and ¯h=−0.8. mApproximate Solutions Difference between Two Iterations 1 4.365721 4.365721 35.00683 35.00683 15.39427 15.39427 18.58639 18.58639 45.94279 45.94279 5.020640 5.020640 2 7.081621 2.71589 31.93003 3.07679 15.70567 0.31139 12.3250 30.9114 15.1206 30.82215 5.013390 0.7249 ×10−2 3 7.329236 0.24761 31.62606 0.30397 15.75156 0.45893 ×10−1 2.0715 14.3966 29.53187 14.4112 5.01834 0.495 ×10−2 4 7.33754 0.831 ×10−2 31.61051 0.155 ×10−1 15.75638 0.482 ×10−2 1.0558 3.1274 26.40971 3.1221 5.017514 0.834 ×10−3 5 7.32929 0.825 ×10−2 31.61790 0.738 ×10−2 15.75678 0.39 ×10−3 1.00038 0.5541 ×10−1 26.46585 0.561 ×10−1 5.01755 0.4 ×10−4 6 7.32617 0.311 ×10−2 31.62090 0.2998352 ×10−2 15.75674 0.3 ×10−4 0.91942 0.8095 ×10−1 26.5475 0.8165 ×10−1 5.01759 0.3 ×10−4 7 7.32556 0.6 ×10−3 31.62149 0.59 ×10−3 15.75674 @.0 91948 0.6 ×10−4 26.54801 0.50 ×10−3 5.01759 @.0 Mathematics 2021,9, 1321 23 of 26 Table 3. Numerical results based on the FPA and Condition (1) for t= 0.2, ρ= 1, ¯h=− 0.8 and ε=10−4. mApproximate Solutions Difference between Two Iterations 1 4.365721 4.365721 35.00683 35.00683 15.39427 15.39427 18.58639 18.58639 45.94279 45.94279 5.020640 5.020640 2 7.081621 2.71589 31.93003 3.07679 15.70567 0.31139 12.3250 30.9114 15.1206 30.82215 5.013390 0.7249 ×10−2 3 7.329236 0.24761 31.62606 30397 15.75156 45893 ×10−1 2.0715 14.3966 29.53187 14.4112 5.01834 0.495 ×10−2 4 7.33754 0.831 ×10−2 31.61051 0.155 ×10−1 15.75638 0.482 ×10−2 1.0558 3.1274 26.40971 3.1221 5.017514 0.834 ×10−3 5 0.732929 0.825 ×10−2 31.61790 0.738 ×10−2 15.75678 0.39 ×10−3 1.00038 0.5541 ×10−1 26.46585 0.561 ×10−1 5.01755 0.4 ×10−4 Table 4. Number of iterations based on the FPA and Condition (1) for various ε and t= 0.2, ρ= 1, ¯h=−0.8. εSmall Values 10−510−410−10.5 Large Values m>> 9 9 5 3 1 1 Mathematics 2021,9, 1321 24 of 26 Table 5. Numerical results based on the DSA and using the CESTAC method for t= 0.2, ρ= 0.9 and ¯h=−0.6. mApproximate Solutions Difference between Two Iterations 1 0.6416038 0.6416040 35.75721 35.75721 14.30519 14.30519 21.33839 21.33840 47.63999 47.63999 5.294120 5.294120 2 8.290508 1.87447 30.37535 5.38186 14.48619 0.18099 34.01543 12.67704 5.92429 41.757 5.18415 0.10996 3 7.42754 0.86296 31.0735 0.69816 14.44400 0.4218 ×10−1 34.8834 0.86797 ×10−2 60.98040 0.17375 ×10−1 5.373754 0.18960 . . .. . .. . . 9 7.61575 0.32 ×10−3 30.7691 0.1 ×10−3 14.42713 0.19 ×10−3 19.1588 0.43 ×10−3 45.76194 0.2 ×10−4 5.33333 0.6258 ×10−1 10 7.61582 0.7 ×10−4 30.76886 @.0 14.4270 0.7 ×10−4 19.16541 0.6 ×10−4 45.76194 @.0 5.28554 0.4779 ×10−1 9. Conclusions Given the importance of modeling and controlling the COVID-19 pandemic, we focused on a nonlinear fractional order model of COVID-19. We applied the Caputo–Fabrizio fractional derivative to present a novel fractional order model. Some theorems were proven to show the existence of a solution and stability analysis. The HATM was used to solve the model numerically and the CESTAC method and the CADNA library were applied to validate the results. In order to show the accuracy of the results, instead of applying the absolute error we used the difference between two successive approximations. The main theorem of the CESTAC method enabled us to use the new termination criterion. This method is based on the DSA and using this method we can find the optimal iteration and the optimal approximation of the HATM for solving the fractional order model. A comparative study between the FPA and the DSA was presented to show the abilities of the CESTAC method and the CADNA library. Author Contributions: Conceptualization, S.N., S.M. and J.J.N.; methodology, S.N.; software, S.N.; validation, S.N., S.M. and J.J.N.; formal analysis, S.N., S.M. and J.J.N.; investigation, S.N., S.M. and J.J.N.; resources, S.N.; data curation, S.N.; writing—original draft preparation, S.N.; writing—review and editing, S.N., S.M. and J.J.N.; visualization, S.N., S.M. and J.J.N.; supervision, S.N.; project Mathematics 2021,9, 1321 25 of 26 administration, S.N.; funding acquisition, S.N. All authors have read and agreed to the published version of the manuscript. Funding: The work of J.J.N. has been partially supported by the Xunta de Galicia under grant ED431C 2019/02, as well as by Instituto de Salud Carlos III and the Ministerio de Ciencia e Innovación of Spain, research grant COV20/00617. The work of S. Noeiaghdam has been supported by a grant from the Academic Council in the direction of the scientific school of Irkutsk National Research Technical University No. 14-NSH-RAN-2020. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no conflict of interest. References 1. Coronavirus Disease (COVID-19) Pandemic. Available online: https://www.who.int/emergencies/diseases/novel-coronavirus- 2019# (accessed on 2 March 2021). 2. Srivastava, H.M.; Area, I.; Nieto, J.J. Power-series solutions of compartmental epidemiological models. Math. Biosci. Eng. 2021 ,18. [CrossRef] 3. Ndaïrou, F.; Area, I.; Nieto, J.J.; Silva, C.J.; Torres, D.F.M. Fractional model of COVID-19 applied to Galicia, Spain and Portugal. Chaos Solitons Fractals 2021,144, 110652. [CrossRef] 4. Bushnaq, S.; Saeed, T.; Torres, D.F.M.; Zeb, A. Control of COVID-19 dynamics through a fractional-order model. Alex. Eng. J. 2021,60, 3587–3592. [CrossRef] 5. Hu, F.; Huang, M.; Sun, J.; Zhang, X.; Liu, J. An analysis model of diagnosis and treatment for COVID-19 pandemic based on medical information fusion. Inf. Fusion 2021,73, 11–21. [CrossRef] 6. Gebremeskel, A.A.; Berhe, H.W.; Atsbaha, H.A. Mathematical modeling and analysis of COVID-19 epidemic and predicting its future situation in Ethiopia. Results Phys. 2021,22, 103853 [CrossRef] 7. Inc, M.; Acay, B.; Berhe, H.W.; Yusuf, A.; Khan, A.; Yao, S.-W. Analysis of novel fractional COVID-19 model with real-life data application. Results Phys. 2021,23, 103968. [CrossRef] 8. Zhang, W.; Zhang, C.; Bi, Y.; Yuan, L.; Jiang, Y.; Hasi, C.; Zhang, X.; Kong, X. Analysis of COVID-19 epidemic and clinical risk factors of patients under epidemiological Markov model. Results Phys. 2021,22, 103881. [CrossRef] 9. Galanis, G.; Hanieh, A. Incorporating Social Determinants of Health into Modelling of COVID-19 and other Infectious Diseases: A Baseline Socioeconomic Compartmental Model. Soc. Sci. Med. 2021,274, 113794 [CrossRef] 10. Panwar, V.S.; Sheik Uduman, P.S.; Gómez-Aguilar, J.F. Mathematical modeling of coronavirus disease COVID-19 dynamics using CF and ABC non-singular fractional derivatives. Chaos Solitons Fractals 2021,145, 110757. [CrossRef] 11. Naik, P.A.; Zu, J.; Owolabi, K.M. Modeling the mechanics of viral kinetics under immune control during primary infection of HIV-1 with treatment in fractional order. Phys. A Stat. Mech. Appl. 2020,545, 123816. [CrossRef] 12. Ghanbari, B.; Günerhan, H.; Srivastava, H.M. An application of the Atangana–Baleanu fractional derivative in mathematical biology: A three-species predator-prey model. Chaos Solitons Fractals 2020,138, 109910. [CrossRef] 13. Günerhan, H.; Dutta, H.; Dokuyucu, M.A.; Adel, W. Analysis of a fractional HIV model with Caputo and constant proportional Caputo operators. Chaos Solitons Fractals 2020,139, 110053. [CrossRef] 14. Noeiaghdam, S.; Sidorov, D. Caputo–Fabrizio Fractional Derivative to Solve the Fractional Model of Energy Supply-Demand System. Math. Model. Eng. Probl. 2020,7, 359–367. [CrossRef] 15. Naik, P.A.; Yavuz, M.; Zu, J. The role of prostitution on HIV transmission with memory: A modeling approach. Alex. Eng. J. 2020 , 59, 2513–2531. [CrossRef] 16. Naik, P.A.; Zu, J.; Owolabi, K.M. Global dynamics of a fractional order model for the transmission of HIV epidemic with optimal control. Chaos Solitons Fractals 2020,138, 109826. [CrossRef] 17. Suleman, M.; Lu, D.; He, J.H.; Farooq, U.; Noeiaghdam, S.; Chandio, F.A. Elzaki Projected Differential Transform method for Fractional order System of Linear and Nonlinear Fractional Partial Differential Equation. Fractals 2018,26, 1850041. 18. Noeiaghdam, S. A novel technique to solve the modified epidemiological model of computer viruses. SeMA J. 2019 ,76, 97–108. [CrossRef] 19. Hattaf, K.; Mohsen, A.A.; Harraq, J.; Achtaich, N. Modeling the dynamics of COVID-19 with carrier effect and environmental contamination. Int. J. Model. Simul. Sci. Comput. 2021. [CrossRef] 20. Abboubakar, H.; Kumar, P.; Rangaig, N.A.; Kumar, S. A malaria model with Caputo–Fabrizio and Atangana–Baleanu derivatives. Int. J. Model. Simul. Sci. Comput. 2021 12, 2150013. [CrossRef] 21. Liao, S.J. Beyond Perturbation: Introduction to Homotopy Analysis Method; Chapman & Hall/CRC Press: Boca Raton, FL, USA, 2003. 22. Liao, S.J. Homotopy Analysis Method in Nonlinear Differential Equations; Higher Education Press: Beijing, China; Springer: Berlin/Heidelberg, Germany, 2012.