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

: In this paper, a nonlinear fractional order model of COVID-19 is approximated. For this aim, at ﬁrst 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 ﬁxed point theorem and the Picard–Lindelof approach are used. Additionally, the stability analysis is discussed using the ﬁxed point theorem. The model is approximated based on Indian data and using the homotopy analysis transform method (HATM), which is among the most famous, ﬂexible 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 ﬂoating-point 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.


Introduction
Corona virus infection  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 where the exact and approximate solutions are denoted by g(x) and g m (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: where g m (x) and g m−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][32][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][34][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 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.

Preliminaries
In this section, we present some definitions and details of fractional derivatives [43,44]. Definition 1 ([43]). Let g ∈ C n . The ω-th order Caputo fractional derivative can be defined via integrable differentiations as follows:

Definition 2 ([43]).
The ω-th order Caputo-Fabrizio derivative for b > a, g ∈ H 1 (a, b), and ω ∈ (0, 1) is defined in the following form: where M(ω) depending on ω denotes a normalization function and M(0) = M(1) = 1. For a function g that is not in H 1 (a, b), this derivative can be presented for g ∈ L 1 (˘∞, b) as follows: In addition, the ω + n-th order fractional derivative CF D ω+n for n ≥ 1 and ω ∈ (0, 1) can be defined as CF D ω+n g(t) := CF D ω (D n g(t)).

Definition 4 ([44]).
The ω-order Riemann-Liouville fractional integral for Re(ω) > 0 is specified as: Definition 5 ([44]). The fractional integral of Caputo-Fabrizio is defined by Additionally, the left and right fractional integrals of ( CF a D ω ) are defined [45] respectively by: We have the following definitions for the Sumudu transform, x which is based on the classical Fourier integral [46][47][48].
Then ST[g(t); u] = G(u) denotes the Sumudu transform of a function g(t) ∈ A that is defined as 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

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 (g n (x 0 )) n∈N converges to x * for all x 0 ∈ X.
Definition 10. Consider the Banach space (Y, . ), a self-map G on Y, and the recursive relation P n+1 = φ(G, P n ). Let Ω(G) be the fixed point set of G, which includes Ω(G) = ∅ and lim n→∞ P n = p ∈ Ω(G). Suppose that {g n } ⊂ Ω and e n = g n+1˘φ (G, g n ), if lim n→∞ e n = 0 implies that lim n→∞ g n = p, then the recursive procedure P n+1 = φ(G, P n ) is G-stable. Suppose that our sequence g n has an upper boundary. If Picard's iteration P n+1 = GP n is satisfied in all these conditions, then P n+1 = GP n is G-stable.
for all x, y ∈ Y, where R ≥ 0 and 0 ≤ r < 1. Then G is Picard G-stable.

Nonlinear Fractional Order Model of COVID-19
Consider the following non-linear bio-mathematical model [51,52]: 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 β = R 0 γ where R 0 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 q t . 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:

Existence of Solution
Applying Definition 5 to both sides of Equation (4), we obtain: By assuming S 0 (t) = u 1 (t), E 0 (t) = u 2 (t), I 0 (t) = u 3 (t), Q 0 (t) = u 4 (t), R 0 (t) = u 5 (t) and D 0 (t) = u 6 (t), these equations can be written in the following form: [dI n (s) + dQ n (s)]ds, (6) and by letting n approach infinity in Equation (6), the approximate solutions can be obtained as: lim n→∞ S n (t) = S(t), lim n→∞ Q n (t) = Q(t), 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: where Assuming a uniform norm on C[a, ρ i ], i = 1, . . . , 5, and defining the Picard operator as we can write where and thus Z(t) ∞ ≤ max{ρ 1 , . . . , Assume that L = max{L 1 , . . . , L 6 } and there is t 0 such that t 0 ≥ t and where Now, we can write and using definition of the Picard operator we obtain where λ < 1. We know that G is a contraction so µλ < 1 and O is a contraction, showing that the proof is complete.

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: and based on the definition of the Sumudu transform we can write and Finally, the following recursive relation can be obtained as where the approximate solutions can be found using the following relations Fixed Point Theorem for Stability Analysis of the Iteration Method Theorem 2. Let F be a self-map that is defined in the following form: Formula (26) Proof. In order to prove the theorem, first the following relation should be computed for (n, m) ∈ N × N to show that F has a fixed point as and taking the norm on both sides of Equation (28) we obtain: Because of the same role of both solutions, we obtain: and applying (29) and (30) we can write We know that S n , E n , I n , Q n , R n , D n are bounded because they are convergent sequences. Thus for all t, there are values P 1 , P 2 , P 3 , P 4 , P 5 , P 6 such that S m < P 1 , E n < P 2 , I n < P 3 , Q n < P 4 , R n < P 5 , D n < P 6 , where (m, n) ∈ N × N. Now, based on Equations (31) and (32) we obtain: where k i are obtained from ST −1 1−ω+ωu M(ω) ST . . By repeating the process we have where Then, the F self-mapping has a fixed point. In addition, we show that F satisfies the conditions in Theorem 1. By assuming that (33) and (34) hold, we have R = (0, 0, 0, 0, 0) and Then the conditions of Theorem 1 are satisfied and the proof is concluded.

Application of the HATM to Solve the Model
In order to apply the HATM for solving the fractional order model (4), applying the Laplace transformation we obtain: Additionally, we can write and we have We define the nonlinear operators as follows: Based on the traditional HAM, the zero-order deformation equation can be defined as t; p), . . . , ϕ 6 (t; p)), (1 − p)L[ϕ 5 (t; p) − R 0 (t)] = phH(t)(ϕ 1 (t; p), . . . , ϕ 6 (t; p)), (1 − p)L[ϕ 6 (t; p) − D 0 (t)] = phH(t)(ϕ 1 (t; p), . . . , ϕ 6 (t; p)), where p ∈ [0, 1],h, H(t) and L are the embedding parameter, auxiliary convergence control parameter, auxiliary function and the linear operator, respectively. By assuming S 0 (t), E 0 (t), I 0 (t), Q 0 (t), R 0 (t) and D 0 (t) as initial guesses and increasing p from 0 to 1 as follows: the solutions of the problem can be found using the initial guesses in the form of a Taylor series where . Thus, the series (43) will be convergent to the exact solution by choosing the suitable values of auxiliary parameters and functions as The m-th order deformation equation can be written as follows: where and using the inverse Laplace transforms we can write where we will apply this relation to find the successive iterations of the HATM.

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 * ∈ B for v * ∈ R with α mantissa bits of the binary FPA as where the sign, missing segment of the mantissa and the binary exponent of the result are denoted by ψ, 2 −α η and E, respectively [31][32][33][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Ṽ * = ∑ r i=1 V * 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Ṽ * , can be found by 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Ṽ * = 0, or CṼ * ,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ã 1 andã 2 , the number of common significant digits can be defined as be the approximate solution of the mathematical model of COVID-19 (3) that is produced by the HATM; then Proof. Using the mentioned definition we can write 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 According to the traditional HAM [21][22][23], we can write Therefore, It is clear that O( 1 m ) << 1, and thus the right-hand sides of the above relations decrease as m increases.

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 andh such that the parameter h can help identify and control the convergence region. To this aim, by plotting somē h-curves we can find these regions. In Figure 2, theh-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 theh-curve with the axis x. Additionally, by increasing the number of iterations of the HATM, we can plot more accurateh-curves. In Figure 3 In Table 2, the numerical results obtained based on the DSA and using the CESTAC method are shown for t = 0.2, ρ = 1 andh = −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 m opt = 7 and the optimal approximations are: S opt = 0.732556E + 001, Q opt = 0.91948E + 000, E opt = 0.3162149E + 002, R opt = 0.2654801E + 002, 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 andh = −0.6 based on the CESTAC method. Using this table we can find the optimal iteration m opt = 10 and the following optimal approximations: S opt = 0.761582E + 001, Q opt = 0.1916541E + 002, E opt = 0.3076886E + 002, R opt = 0.4576194E + 002, 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.

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. 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.