A Comparative Study between Discrete Stochastic Arithmetic and Floating-Point Arithmetic to Validate the Results of Fractional Order Model of Malaria Infection

to Validate the Results of Model Abstract: The researchers aimed to study the nonlinear fractional order model of malaria infection based on the Caputo-Fabrizio fractional derivative. The homotopy analysis transform method (HATM) is applied based on the ﬂoating-point arithmetic (FPA) and the discrete stochastic arithmetic (DSA). In the FPA, to show the accuracy of the method we use the absolute error which depends on the exact solution and a positive value ε . Because in real life problems we do not have the exact solution and the optimal value of ε , we need to introduce a new condition and arithmetic to show the efﬁciency of the method. Thus the CESTAC (Controle et Estimation Stochastique des Arrondis de Calculs) method and the CADNA (Control of Accuracy and Debugging for Numerical Applications) library are applied. The CESTAC method is based on the DSA. Also, a new termination criterion is used which is based on two successive approximations. Using the CESTAC method we can ﬁnd the optimal approximation, the optimal error and the optimal iteration of the method. The main theorem of the CESTAC method is proved to show that the number of common signiﬁcant digits (NCSDs) between two successive approximations are almost equal to the NCSDs of the exact and approximate solutions. Plotting several graphs, the regions of convergence are demonstrated for different number of iterations k = 5, 10. The numerical results based on the simulated data show the advantages of the DSA in comparison with the FPA.


Introduction
Malaria is a life-threatening disease caused by a parasite. The causative agent of this disease to humans is the bite of infected female Anopheles mosquitoes. But malaria is a preventable and treatable disease. According to studies, malaria can be caused by five types of parasites in humans, so that two of them can be a dangerous threat to humans. In 2019, almost half of the world's population was at risk for malaria. Sub-Saharan Africa has the highest mortality rate. In 2019, there were about 229 million malaria cases, and the estimated number of malaria deaths was 409,000 [1].
Given the above, it is easy to understand the importance of examining this disease. Accordingly, various mathematical models have been proposed to control and predict the behavior of this disease according to various human and environmental factors. The optimal control for a fractional order malaria transmission model has been studied by Sweilam et al. [2]. The comparison of mathematical frameworks for modeling erythropoiesis in the context of malaria infection and the numerical simulation of a reactiondiffusion malaria infection model using B-splines collocation have been illustrated by Fonseca et al. [3] and Mittal et al. [4]. The effects of climate change on Plasmodium vivax malaria transmission dynamics based on mathematical models and the endectocide-treated cattle for malaria control are studied by Kim et al. [5] and Yakob [6]. Recently, many authors have focused on solving the mathematical models. Naik et al. [7] have discussed the fractional model of human immunodeficiency virus (HIV) infection under immune control. The HIV transmission with memory and the optimal control on fractional model of HIV infection has been illustrated by Naik et al. [8,9]. The modified epidemiological model of computer viruses and the homotopy analysis method (HAM) has been applied to solve this model by Noeiaghdam et al. [10,11]. The fractional model of energy systems to control the supply and demand of energy has been explained by Noeiaghdam et al. [12] and solving the fractional system of partial differential equations has been demonstrated by Suleman et al. [13].
There are many mathematical and engineering problems that can be solved using the HAM. Liao for the first time has presented the HAM to solve the linear and nonlinear problems [14][15][16]. As we know, in this method we are free to choose the auxiliary functions and parameters that it is one of abilities of this method than other methods. The HAM for solving the singular Volterra integral equations [17], solving integral equations using the homotopy-regularization method [18] and solving integral equations based on finding the optimal convergence control parameter of the HAM [19] are only some of the applications of the HAM and HATM for solving linear and nonlinear problems.
We should know that all mathematical methods that we apply to solve the problems are based on the FPA and we apply the mathematical applications such as Mathematica, Maple and Matlab to provide the accurate results. In the FPA, when we want to show the accuracy of the method we apply the traditional absolute error. Also in order to control the accuracy, this condition depends on a small positive value ε. But we should know that in the real life problems we do not have the exact solution. On the other hand, we do not know the optimal value of ε. Thus we will have one or two iterations for large values of ε without producing the accurate results and large number of iterations for small values of ε without improving the accuracy. This is the big fault of the mathematical methods based on the FPA.
Because of these problems we introduce the DSA instead of the FPA. In the DSA, we should apply the CESTAC method to validate the numerical results. The mathematical methods can be implemented on the DSA using the CADNA library. In this method, we do not need to apply the absolute error. The termination criterion of the CESTAC method is based on two successive approximations and in the right hand side instead of ε we have @.0. This sign is the informatical zero. It can be produced only in the CESTAC method by the CADNA library and it shows that the NCSDs between two successive approximations are almost equal to zero [20][21][22][23]. The CADNA library should be done on Linux operating system and all CADNA codes should be written using C, C++, FORTRAN or ADA codes [24]. Recently, many methods have been validated using the CESTAC method and the CADNA library to solve various problems. This method has been applied to validate the numerical results of the homotopy analysis method, the homotopy perturbation method and the Admoian decomposition method for solving different kinds of integral equations [19,25,26]. Dynamical validation of fuzzy integral equations [27], control of accuracy on the load leveling problem [28], solving the reverse osmosis system [29,30], validation of collocation method for various problems [31], validation of mathematical models [32] and finding the optimal regularized parameter of the regularization method [18] are some applications of the CESTAC method and the CADNA library.
In this paper the fractional order model of malaria infection is introduced. The HATM is applied to solve the model. For this aim, the Laplace transformations as the linear operator of the HATM is used. Also the CESTAC method and the CADNA library are introduced and the advantages of the CESTAC method are discussed. The main theorem of the CESTAC method is proved. Using the method, the approximate solution of the problem is obtained and applying the CADNA library the numerical results are validated. Plotting severalh-curves the convergence regions are found. The numerical results are presented based on the FPA and the DSA. Comparing the results, the advantages of the DSA than the FPA to solve the fractional order of malaria infection can be found.

Preliminaries
In this section, we present some definitions and details of the fractional derivatives [33][34][35].
Definition 1. Let g be a continuous function then the ρ-th order Caputo fractional derivative can be defined via integrable differentiations as follows

Definition 2.
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(ρ) shows a normalization function and M(0) = M(1) = 1. For function g which is not member of H 1 (a, b), this derivative can be presented for g ∈ L 1 (˘∞, b) as follows: Also, ρ + 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 3.
For 0 < ρ ≤ 1 and M(ρ) = 1, the Laplace transform for the Caputo-Fabrizio derivative can be presented in the following form

Definition 5. The fractional integral of Caputo-Fabrizio is defined by
Also, the left and right fractional integrals of ( CF a D ρ ) are defined [36] respectively by

Fractional Order Model of Malaria Infection
In this section, we consider the following fractional order malaria infection model [6] where 0 < ρ ≤ 1 is the order of Caputo-Fabrizio fractional derivative and Ψ Z shows the infectious mosquitoes, Ψ S is the bited susceptible hosts and Ψ I is the infected hosts where the infected hosts before reversion to fully susceptible at rate τ can have level of temporary immunity (1 − θ, where 0 ≤ θ ≤ 1). The hosts can become asymptomatically infected Ψ A during the recovered state Ψ R . Symptomatic infection is a temporary condition that may have different transmission potentials to (σ) carriers, after which people return to the improved group at a rate of κ. In this model, considering that a fixed human population is assumed, human mortality is compensated by births and we get Ψ S = 1 − (Ψ I + Ψ R + Ψ A ). Sensitive (Ψ X ) vectors begin to become infected after being bitten by an infected infectious host, but they are not yet infected (Ψ Y ), but after the extrinsic incubation period ( 1 ξ ), become infectious (Ψ Z ). We apply m to show the ratio of vector-to-hosts which vectors typically outnumber hosts. Here, a stable vector population is assumed whereby births are set to balance deaths (µ v ) and so Ψ X = 1 − (Ψ Y + Ψ Z ). All parameters and values are detailed in Table 1 based on the mentioned data in [6].

Homotopy Analysis Transform Method
Applying the Laplace transform for both sides of Equation (1) we get and we can write Now, we have We define the nonlinear operator as follows and using Equation (5), the zero order deformation equations can be defined as where p ∈ [0, 1] is the embedding parameter,h = 0 and H(t) = 0 are the auxiliary parameter and function and L shows an auxiliary linear operator. Also, ϕ j (t; p) for j = 1, 2, 3, 4, 5 are unknown functions and thus clearly we can see that by increasing p from 0 to 1, the solution can be produced from initial guesses where . By choosing suitable values for auxiliary functions and parameters, series (8) will be convergent at p = 1 and we have The k-th order deformation equation can be constructed as follows where Applying inverse Laplace transform for both sides of the k-th order deformation Equation (10) we get and by substituting the suitable values of χ k as the approximate solutions will be obtained using the following series

Discrete Stochastic Arithmetic-CESTAC Method
In the mathematical methods based on the FPA, researchers have applied the absolute error or difference between to discuss the accuracy of the methods as follows where Ψ I (t) and Ψ I k (t) are the exact and approximate solutions of the problem and ε is a small positive value. But this condition can not be acceptable. In the real life problems we do not know the exact solution. Also, generally the researchers do not know the optimal value of ε. If we choose a small value as ε we will have huge number of iterations without improving the accuracy of the results and we will have so many extra iterations. Also, if we choose a large value as ε, the algorithm will be stopped at the first step without producing accurate results. Thus we propose the CESTAC method and the CADNA library which are based on the SA and instead of the stopping condition (14) we have new termination criterion where Ψ I k (t) and Ψ I k−1 (t) are two successive approximations and @.0 denotes the informatical zero [24,[37][38][39].
In order to apply the CESTAC method, we collect all representable values by computer in set B. For w * ∈ R with α mantissa bits we can write the value W * ∈ B based on the binary FPA as where ρ, 2 −α φ and E show the sign, missing segment of the mantissa and the binary exponent of the result respectively [20,21,24,37]. By changing value α from 24 to 53, the accuracy of the results can be changed from single to double precisions. Let φ be a casual variable which is uniformly distributed on [−1, 1]. Making perturbation on the last mantissa bit of W * , the mean (µ) and the standard deviation (σ) can be found for the results of W * . It means that if we have k samples of W * as Φ = W * 1 , W * 2 , ..., W * k and make perturbation on the last mantissa bit then we will haveW Now, we will be able to find the NCSDs of W * andW * as CW * ,W * = log 10 where τ δ is the value of T distribution as the confidence interval is 1 − δ, with k − 1 freedom degree. In the CESTAC method, if we haveW * = 0 or CW * ,W * ≤ 0, the algorithm will be stopped by showing the informatical zero sign as W * = @.0. The mentioned process of the CESTAC method will be done using the CADNA library. Thus we do not need to produce the samples or calculate the mean and standard deviation values. We should apply the CADNA codes to write the solving process and after that all steps of the CESTAC method will be done using this library. For more details of this method and the CADNA library please see [20,21,24,37].
We have the following highlights on the CESTAC method and the CADNA library [22,23,38]: Applying this method the optimal number of iterations, optimal approximation and error can be found. • The CADNA library can find some of the numerical instabilities. • The CADNA library should be run on LINUX operating system. • The CADNA codes should be written using C, C++, FORTRAN or ADA codes.

Definition 6.
For two real numbers ω 1 and ω 2 , the NCSDs can be defined [37,38] as C ω 1 ,ω 1 = +∞. (17) Theorem 1. Assume that the series (13) are the approximate solutions of the mathematical model of Malaria infection (1) which are produced by the HATM, then the NCSDs between two successive approximations are almost equal to the NCSDs between the exact and approximate solutions.
Proof. Using Definition 6 we have By increasing the iterations number k, we can see in the first term of Equation (18), the approximate and exact solutions Ψ I k , Ψ I are close together and we can neglect that. For the second term we have Based on the traditional HAM we can write Therefore, we get By repeating the process for other terms we have Since O( 1 k ) << 1, then the right hand side of above relation decreases as k increases. Thus, the proof is complete.

Results and Discussions
In this section, the numerical results of validation of the HATM based on both the SA and the FPA are presented using the simulated data. We know that because of existence of auxiliary parameters and functions in the HATM, this method is one of the flexible methods among semi-analytical methods. It can help us to make accurate solutions with fast convergence rate. Applying the HATM to solve the fractional model (1), the obtained solution will be based on t andh. Thus for various values of t we will be able to find the solutions based onh and plot theh-curves. Applying these graphs we can find the convergence region of the solution. The parallel part ofh-curve with axis x will be the convergence region of the problem. We apply the presented values in Table 1 based on the mentioned data in [6]. The following approximate solutions are obtained for k = 5: .59242 × 10 −6h5 t 5 , Ψ Y 5 (t) = 100 + 138.5ht + 277h 2 t + 277h 3 t + · · · + 0.0615969h 5 t 4 + 0.000524847h 5 t 5 , Ψ Z 5 (t) = 55 − 22.5ht − 45h 2 t − 45h 3 t + · · · − 0.0453403h 5 t 4 − 0.000534686h 5 t 5 .
Now the convergence regions can be obtained. In Figures 1 and 2 the convergence regions based on the approximate solutions for k = 5, 10 andh are plotted. According to these figures the regions for k = 5, ρ = 1 and t = 1 are In Figure 3, the residual error functions based onh for t = 1 and k = 5, 10 are plotted. Using this graph, obviously we can see the regions of convergence for different functions. For k = 5 and −1.2 h −0.7 we will have more accurate results than other regions. Also, for k = 10 we have −1.3 h −0.5. Figure 4, shows the residual error functions based on t forh = −1 and k = 5, 10. Completely clear that by increasing number of iterations k from 5 to 10 the HATM can produce more accurate results. In Table 2, the numerical results are obtained based on the FPA forh = −1 and ε = 10 −10 . But in this table we do not know the optimal value of ε. It means that we do not know for ε = 10 −10 we can provide the accurate results or no. In Table 3, obviously the number of iterations for various ε are demonstrated. For large values of ε we can not provide accurate results and the algorithm is stopped soon. On the other hand, for small values of ε we have a large number of iterations without improving the accuracy. Because of the mentioned problems, we use the CESTAC method and the CADNA library which are based on the SA. Table 4 shows the results applying the CESTAC method. Thus instead of applying the absolute error or residual error, by finding difference of two successive approximations we can show the accuracy of the results and the main Theorem 1 supports us to do it. Using this technique, not only we do not need to have the exact solution but also we can find the optimal iteration and the optimal approximations. Thus, the optimal iteration is k opt = 9 and the optimal approximations are The sign @.0 in the ninth step shows that the NCSDs for two successive approximations is almost equal to the NCSDs for exact and approximate solutions. But this sign can not be produced in the FPA and in the mathematical methods based on the FPA we need to provide a small positive value ε. Since we do not know the optimal value of ε, it can make problem in computations.  Table 2. The numerical results based on the FPA for t = 0.5, ρ = 1,h = −1 and ε = 10 −10 .   Table 4. The numerical results using the CESTAC method for t = 0.5, ρ = 1 andh = −1.

Conclusions
The main theorem of the CESTAC method has been proved which can support us to use the new termination criterion. In the numerical results by plotting severalh-curves we showed the convergence regions. According to the obtained results, we can see the advantages of the DSA than the FPA. In the DSA, we do not need to make more and extra iterations. Also, we can find the optimal iteration and the optimal approximations. But in the FPA, we do not have these abilities. More over in the DSA, the informatical zero sign @.0 can be produced which shows that the NCSDs for two successive approximations are almost equal to the number of common significant digits for exact and approximate solutions and then the numerical algorithm can be stopped. But in the FPA, the stopping condition depends on a small positive value ε that we do not have information about its optimality. Applying the method on the fractional models of HIV infection and COVID-19 are among our future plans.