Error Estimation of the Homotopy Perturbation Method to Solve Second Kind Volterra Integral Equations with Piecewise Smooth Kernels: Application of the CADNA Library

This paper studies the second kind linear Volterra integral equations (IEs) with a discontinuous kernel obtained from the load leveling and energy system problems. For solving this problem, we propose the homotopy perturbation method (HPM). We then discuss the convergence theorem and the error analysis of the formulation to validate the accuracy of the obtained solutions. In this study, the Controle et Estimation Stochastique des Arrondis de Calculs method (CESTAC) and the Control of Accuracy and Debugging for Numerical Applications (CADNA) library are used to control the rounding error estimation. We also take advantage of the discrete stochastic arithmetic (DSA) to find the optimal iteration, optimal error and optimal approximation of the HPM. The comparative graphs between exact and approximate solutions show the accuracy and efficiency of the method.


Introduction
The problem of finding approximate solution for linear Volterra IEs is one of the oldest problems in the applied mathematics researches. Specially, this problem with discontinuous kernel has many applications in the load leveling problems, energy storage with renewable and diesel generation, charge/discharge storages control and others [1][2][3].
There are various methods for solving linear and nonlinear problems [4][5][6][7][8][9][10] specially the Volterra IEs with discontinuous kernel. Muftahov et al. in [11] applied the Lavrentiev regularization and direct quadrature method, Sidorov in [12] used the successive approximations and Noeiaghdam et al. studied the Taylor-collocation method for solving Volterra IEs with discontinuous kernel [1,13]. Also, the nonlinear system of Volterra IE with applications was studied in [14,15]. Furthermore, the existence of a continuous solution depending on free parameters and sufficient conditions for the existence of a unique continuous solution of the system of Volterra IE with discontinuous kernels were derived in [16]. The class of integral operator equations of Volterra type with applications to p-Laplacian equations was illustrated in [17]. The problem of generalized solution (in the Sobolev-Schwartz sense) to the Volterra equations with piecewise continuous kernel was illustrated in [18]. Belbas and Bulka in [19] considered the multiple Volterra IEs. The problem of global solution's existence and blow-up of nonlinear Volterra IEs were discussed in [20]. For systematic study of the qualitative theory of Volterra IE with discontinuous kernels readers may refer to monograph [21] and part 1 in monograph [22].
The parametric continuation method for the first time was justified by Bernstein [23] for partial differential equations. Here readers may also refer to excellent review by Lusternik [24]. In community of numerical analysts the parametric continuation method is known as the HPM. This method is among of semi-analytical methods that was popularized by J.H. He [25][26][27]. Then, this method has been extended by many other researchers for solving different problems. The HPM was applied to find the approximate analytical solution of the Allen-Cahn equation in [28], to study the maximum power extraction from fractional order doubly fed induction generator based wind turbines in [29], dissipative nonplanar solitons in an electronegative complex plasma in [30] and others [31][32][33]. Convergence of the parameter continuation method in the homotopy method based on the theorem of V.A. Trenogin (see [34], Section 14, p. 146) will be global with respect to a parameter if there is an a priori estimate of the solution for all values of the parameter (this condition can be replaced with a more stringent requirement for the existence of a unique solution bounded for all values of the parameter). If there is no a priori estimate of the solution, then on the basis of the inverse operator theorem (see [34] p. 135), at least local convergence in the homotopy method can be guaranteed. Due to the models complexity, we addressed only some classes of the results in this field. Many other interesting results concerning nonlinear equations with discontinuous symmetric kernels with application of group symmetry have remained beyond this paper. Results of present paper in combination with methods of representation theory and group analysis in the bifurcation theory [35,36] make it possible to construct solutions of nonlinear models with discontinuous kernels using the HPM.
In the mentioned studies and many other researches, the numerical results have been obtained from the floating point arithmetic (FPA) and the accuracy of the method has been discussed using the traditional absolute error as follows where w(t) and w n (t) are the exact and approximate solutions. This condition depends on the existence of the exact solution and optimal value of ε. Also, based on condition (1) we will not be able to find the more accurate approximation because we do not have information about optimal ε and in some cases we do not know the exact solution. For small values of ε, the numerical algorithm can not be stopped and extra iterations will be produced without improving the accuracy. For large values of ε, the numerical algorithm will be stopped in initial steps without producing enough iterations. Moreover, in condition (1), researchers do not have any idea about optimal approximations, optimal errors or numerical instabilities. The aim of this study is to apply the HPM to solve the second kind linear Volterra IEs with jumping kernel and validate the numerical results using the CESTAC method [37][38][39][40].
In this method, instead of applying the condition (1), we need to produce other and better condition without having the disadvantages of (1) as follows: where @.0 is the informatical zero sign [41] and w n (t) and w n+1 (t) are two successive approximations. Condition (2) is based on the DSA and Theorem 2 can support us to apply this condition theoretically. In this condition, not only we do not need to have the exact solution but also we would be able to identify the optimal approximation, optimal iteration and optimal error of numerical procedure. Also, the CADNA library is applied as an important software for this validation. The CESTAC method and the CADNA library have been introduced and developed during decades by researchers from LIP6, the computer science laboratory in Sorbonne University in Paris, France (https://wwwpequan.lip6.fr/). This principle was introduced in [38] and it was extended to various quadrature rules in [42][43][44] and others [45,46]. The CADNA library should be done on the LINUX operating system and its codes should be written using C, C++ or ADA codes [40,[47][48][49]. The CESTAC method is based on the DSA and instead of applying the absolute error to show the precision of method, a termination criterion is applied based on two successive approximations [50][51][52][53]. Thus in this technique we do not need to have the exact solution to compare the results. Also, we will prove that number of common significant digits (NCSD) of two successive approximations are almost equal to the NCSD of exact and approximate solutions. So the new theorem gives the license to apply the new stopping condition instead of previous one. This technique has some advantages than other methods based on the FPA [37,39,50,52,53]. Due to the advantages of the CESTAC method we can find the optimal iteration of iterative and numerical methods, optimal approximation and optimal error. Furthermore, the extra iterations can be neglected and some of numerical instabilities can be detected too [13,[54][55][56].
This paper is arranged as follows: In the next section, the preliminaries are described regarding to the HPM. In third section, the DSA and the CESTAC method are discussed. Also, algorithm of the CESTAC method and sample code of the CADNA library are presented. In forth section the main idea is described. Then using the HPM we solve the second kind linear Volterra IEs with jumping kernel. Furthermore, the convergence theorem is proved. Also, a theorem is presented which proves that instead of traditional absolute error which depends on the exact and approximate solutions, a termination criterion can be applied which depends on two successive approximations. Section five includes some examples. Also, several tables are presented to show the efficiency of method. The last section is conclusion.

Preliminaries
For operator F, given function g and prepared function x we get the following operator equation as We can write the operator F in the following form where the remain part of F showed by N and L is the linear operator. Now, Equation (3) can be presented as According to the traditional homotopy [25][26][27], for parameterâ ∈ [0, 1], the homotopy operator H can be presented as where v(z,â) is defined on Γ × [0, 1] → R and x 0 is the initial guess of Equation (3). Now, by applying Equation (4) we get can be applied to find the solution of H(v,â) = 0. Then comparing the same powers of parameterâ we can find the successive functions v j , j = 0, · · · , n. Finally, applying the solution of Equation (3) can be found and the n-th order approximation is in the following form

Stochastic Arithmetic and the CESTAC Method
In this section, the CESTAC method is described and the algorithm of this method is presented. Also, a sample program of the CADNA library is demonstrated and finally advantages of presented method based on the DSA are investigated in comparison with the traditional FPA [37][38][39][40]50].
Assume that some representable values are produced by computer and they are collected in set A. Then W ∈ A can be produced for w ∈ R with R mantissa bits of the binary FPA in the following form where sign of w showed by χ, missing segment of the mantissa presented by 2 −R ξ and the binary exponent of the result characterized by E. Moreover, there are single and double precisions by choosing R = 24, 53 [40,[50][51][52][53]. Assume ξ is the casual variable that uniformly distributed on [−1, 1]. After making perturbation on final mantissa bit of w we will have (µ) and (σ) as mean and standard deviation for results of W which they have important rule in precision of W. Repeating this process J times for W i , i = 1, . . . , J we will have quasi Gaussian distribution for results. It means that µ for these data equals to the exact w.
It is clear that we should find µ and σ based on W i , i = 1, . . . , J. For more consideration, the following Algorithm 1 is presented where τ δ is the value of T distribution as the confidence interval is 1 − δ with J − 1 freedom degree [52].
Usually, in order to find the numerical results we need to apply the usual packages like Mathematica and Matlab. Here, instead of them we introduce the CADNA library and the CESTAC method to validate the numerical results [1,55,56,62].
This library should run on LINUX operating system and all commands should be written by C, C++, FORTRAN or ADA codes [13,54,59,60,63].
We have many advantages to apply the CESTAC method and the CADNA library instead of traditional schemes using the FPA. In this method, a novel criterion independence of absolute error and tolerance value like ε is presented. Applying the CADNA library, we can find the optimal iteration, approximation and error of numerical methods. Moreover, the numerical instabilities can be identified [13,[54][55][56]. A sample program of the CADNA library is presented as Algorithm 1: Step 1-Make the perturbation of the last bit of mantissa to produce J samples of W as Φ = W 1 , W 2 , ..., W J . Step Step 3-Compute Step 4-Find the NCSDs of w and W ave applying C W ave ,w = log 10 Step Write the main codes here; printf(" %s ",strp(Parameter)); } while(u[n]-u[n-1]!=0); cadna − end(); }

Main idea
Consider the following second kind linear IE where k(t, s) is discontinuous along continuous curves γ i , i = 0, 1, · · · , m − 1 and it can be written in the following form and finally for brief form we get Indeed, the kernel is the principal part of the IE (14). One may think about considered Volterra IE as generalization of classic Duhamel integral. So, the kernel can be understood as instrumental response function (IF, or spectral sensitivity, transmission function, point spread function, frequency response), see e.g., [67]. In this study, we do not focus on specific physical problems, but more on numerical aspects of solutions only.
Based on the HPM and applying Equations (4) and (5) for solving Equation (14), operators L(v) and N (v) should be defined as follows For next step, using Equation (7) the homotopy map can be constructed as follows and we have Now, Equation (17) can be written in the following form where By disjointing the different powers ofâ in both sides of Equation (18) the following successive iterations can be obtained aŝ Applying Equation (10) and successive iterations (19), the approximate solution of Equation (14) can be obtained. (14)

Theorem 1. Assume that functions k i (t, s) and g(t) of Equation
then for initial approximation w 0 which is continuous in [a, b], the series solution (9) will be uniformly convergent to the exact solution for eachâ ∈ [0, 1].
From Equation (20), the following remark can be deduced: Remark 1. Based on the n-th order approximate solution (10), the error function E n = sup t∈[a,b] |w(t) − w n (t)| can be approximated as follows: Order of error E n can be obtained in the following form: where L is a positive real number.
Theorem 2. Let w(t) and w n (t) be the exact and numerical solutions of problem (12) which w n (t) is obtained by using the HPM and Equation (10). Based on assumptions of Theorem 1 and Remark 1 for n enough large we have C w n (t),w n+1 (t) C w n (t),w(t) , where C w n (t),w(t) shows the NCSDs of w n (t), w(t) and C w n (t),w n+1 (t) is the NCSDs of two successive iterations w n (t), w n+1 (t).
Proof. Using Definition 1 and Remark 1 we get Also, Applying Equations (23) and (24) we have From Remark 1, we can write Thus for n enough large we get O 1 n << 1 and consequently C w n (t),w n+1 (t) C w n (t),w(t) .
Theorem 2 shows that when n increases, the NCSDs between two sequential results obtained from the algorithm is almost equal to the common significant digits of the n-th iteration and the exact solution at the given point t which means that for an optimal index like n = n o pt, when w n (t) − w n+1 (t) = @.0 then w n (t) − w(t) = @.0.

Numerical Results
In this section, several examples of second kind linear Volterra IEs with discontinuous kernel are presented. The numerical process is based on the HPM that we discussed in previous sections. Also, using the CESTAC method and the CADNA library for all examples we will arrange some numerical procedures based on the following algorithm to find the optimal approximation, optimal error and optimal step of the HPM for solving linear Volterra IEs with jumping kernels. Having the exact solution in the examples is only to compare the numerical results based on both conditions (1) and (2). Some comparative graphs between exact and approximate solutions are plotted to show the accuracy and efficiency of the method.
Step 2-3-n = n + 1. } Example 1. Consider the following second kind Volterra IE with discontinuous kernel with exact solution w(t) = t 2 − t. Applying the homotopy map (16) and relations (17), (18) and successive iterations (19) we get and finally using series solution (10), the approximate solution for n = 5 can be obtained as follows In this example, in order to show the accuracy of method, the CESTAC method and the CADNA library are applied according to Algorithm 2. Also, instead of applying the termination criterion (1) and using the traditional absolute error, the stopping condition (2) is applied. This condition is based on two successive approximations w n (t) and w n+1 (t). When the difference of these terms is @.0 the CESTAC algorithm will be stopped. It shows that the NCSDs of the difference between two successive iterates is zero. The numerical results using the DSA are presented in Table 1 for t = 0.2 in double precision. According to this table the optimal step of iterations for the HPM is n opt = 10, the optimal approximation is −0.16 and the optimal absolute error is 0.231 × 10 −13 . Figure 1, shows the comparison between the exact and approximate solutions for optimal value n opt = 10 obtained from the CESTAC method.
where the exact solution is w(t) = t 3 8 . Applying the homotopy map (16) and relations (17), (18) and successive iterations (19) we can find the approximate solution in the following form In this example, the DSA and the CADNA library are applied to validate the numerical approximations. Also, using the stopping condition (2) we do not need to have the exact solution to show that accuracy of presented method. The numerical results are presented in Table 2 for t = 0.3 by applying Algorithm 2. Using these results, the optimal approximation is 0.337499999999999 × 10 −2 and the optimal absolute error is 0.26 × 10 −15 and n opt = 11 is the optimal step of iteration for HPM method for solving Example 2. Comparison between the exact and approximate solutions for n opt = 11 is demonstrated in Figure 2.
where the exact solution is w(t) = t 5 .
The numerical results are presented in Table 3. The optimal iteration of the HPM for solving this example is n opt = 10, the optimal approximation is 0.1 × 10 −4 and the optimal error is 0.3 × 10 −19 . To validate the results, the CADNA library is applied based on the termination criterion (2). Theorem 2 is able to permit us to apply the stopping condition instead of the traditional condition (1). In Figure 3, the graph of exact and approximate solutions for optimal value n opt = 10 is studied.

Conclusions
Volterra IEs with discontinuous kernel are among applicable problems in power engineering and especially in load leveling problems. In this study, we applied the HPM while the CESTAC method and the CADNA library used to examine the numerical results. Applying this method not only the optimal iteration of the HPM, the optimal approximation and the optimal error can be found but also some of numerical instabilities can be detected. Furthermore, the substantial theorem is provided which approves the appropriateness of the termination criterion (2) instead of traditional absolute error. We will focus on validating the nonlinear Volterra IEs with discontinuous kernel in fuzzy and crisp forms using the CESTAC method for our future works.