A Novel Algorithm to Evaluate Definite Integrals by the Gauss-Legendre Integration Rule Based on the Stochastic Arithmetic: Application in the Model of Osmosis System

A Novel Algorithm to Evaluate Definite Integrals by the Gauss-Legendre Integration Rule Based on the Stochastic Arithmetic: Application in the Model of Osmosis System

Samad NoeiaghdamMohammad Ali Fariborzi Araghi 

Department of Mathematics, Central Tehran Branch, Islamic Azad University, Tehran 1955847881, Iran

South Ural State University, Lenin prospect 76, Chelyabinsk 454080, Russian Federation

Corresponding Author Email: 
noiagdams@susu.ru
Page: 
577-586
|
DOI: 
https://doi.org/10.18280/mmep.070410
Received: 
20 August 2020
|
Accepted: 
1 December 2020
|
Published: 
18 December 2020
| Citation

OPEN ACCESS

Abstract: 

Finding the optimal iteration of Gaussian quadrature rule is one of the important problems in the computational methods. In this study, we apply the CESTAC (Controle et Estimation Stochastique des Arrondis de Calculs) method and the CADNA (Control of Accuracy and Debugging for Numerical Applications) library to find the optimal iteration and optimal approximation of the Gauss-Legendre integration rule (G-LIR). A theorem is proved to show the validation of the presented method based on the concept of the common significant digits. Applying this method, an improper integral in the solution of the model of the osmosis system is evaluated and the optimal results are obtained. Moreover, the accuracy of method is demonstrated by evaluating other definite integrals. The results of examples illustrate the importance of using the stochastic arithmetic in discrete case in comparison with the common computer arithmetic.

Keywords: 

stochastic arithmetic, CESTAC method, CADNA library, model of osmosis system, Gauss-Legendre integration rule

1. Introduction

Numerical solution of mathematical models are one of the important subjects of applied mathematics. There are many mathematical models that their importance in the real life can not be hidden. Some of these models are the model of computer viruses [1-3], model of HIV infection [4-7], model of smoking habit [8], model of energy supply-demand system [9, 10], nervous stomach mathematical model [11], susceptible-infected-recovered epidemic model [12] and many other models [13].

The reverse osmosis system is a famous technique of water purification. Based on this model, semi-penetrable covers can be applied to the process of isolating and eliminating dissolved particles from water [14]. In converse osmosis, this is done by applying high pressure to the centralized cover so that by applying pressure, pure water flows through the semi-penetrable membrane to the other side at low condensation [15]. This phenomenon can be formulated as an advection-diffusion equation and we can apply it to simulate water and engineering problems such as the model of converse osmosis system [16-18]. In this paper, one of numerical integration rules is applied to find the approximate solution of the osmosis model.

Numerical integrations or quadrature rules are well known schemes to approximate the definite or improper integrals. The general form of a quadrature rule is formed as:

$Q=\int_{\alpha}^{\beta} v(r) S(r) d r \simeq Q_{J}=\sum_{\ell=0}^{J} w_{\ell} S\left(r_{\ell}\right)$,    (1)

where, ν(r) is a weight function, $w_{\ell}, \ell=0,1, \ldots, J$ are coefficients (weights) , $r_{\ell}$ are points of the rule and S is a given function integrable on the interval [α,β] with the weight function v.

The Gaussian integration rule has been applied to estimate the different kinds of functions [19-22]. Deloff [23] applied the G-LIR to approximate the singular integrals and Verlinden [24] used this rule for a definite integral with an endpoint singularity. Also, Lian and Guo [25] worked on interpolation by combining the Gauss-Legendre and the Lobatto-Birkhoff integration rules. Furthermore Elliott and Johnston [26] applied this rule to evaluate the integrals involving the Hankel function.

In order to apply the G-LIR, the floating-point arithmetic (FPA) and the stochastic arithmetic (SA) can be applied. But what are the differences between the FPA and the SA? What are the advantages of the SA? Why do we want to apply the SA instead of the FPA?

In all of the mentioned papers, the authors used the FPA to find the approximate solution without considering optimal iteration, optimal approximation and the validation of numerical results. In these works, the numerical results are obtained for fixed and determined step size. In the FPA, the results may be false or with low accuracy without awareness of the user. The usual termination criterion based on the FPA is in form:

$\left|Q-Q_{J}\right| \leq \varepsilon, o r\left|Q_{J}-Q_{J-1}\right| \leq \varepsilon$,    (2)

where, ε is an arbitrary positive value and $Q_{J}$ is the approximate of Q. This criterion may not be acceptable. Because for large value of ε the process will be stopped quickly and the suitable approximation can not be produced and for small ε many iterations can be generated without progressing the accuracy of the results. Also, the first criterion in (2) is traditional absolute error to be less or equal to ε that depends on the exact solution.

In the SA, in place of (2), the following criterion is applied,

$\left|Q_{J}-Q_{J-1}\right|=@ .0$,    (3)

which depends on two successive approximations $Q_{J}$ and $Q_{J-1}$. Also, sign $@ .0$ is called the informatical zero as it includes the mathematical zero. It means that the difference between these two numerical results has not any significant digits and therefore the number of common significant digits (NCSDs) of $Q_{J}$ and $Q_{J-1}$ is the same.

The aim of this work is to find the optimal iteration, error and approximation of G-LIR for solving the reverse osmosis model using the CESTAC method. La Porte and Vignes [27, 28] described about applications of the SA for reliable scientific computation and error analysis in computing, Chesneaux [29, 30] presented the equality relations in scientific computing and also properties of the SA and he [31] described the CADNA library as an ADA tool for round-off errors analysis and for numerical debugging, Lamotte et al. [32] presented a version of the CADNA for use with C or C++ programs and Alt et al. applied this technique to validate the results of collocation methods for ODEs [33]. Jézéquel et al. studied the infinite integrals using Romberg’s method [34] and a dynamical strategy for approximation methods [35]. Also, they illustrated the CADNA library for estimating round-off error propagation in the paper [36]. Scott et al. [37] presented the numerical study of ’health check’ for scientific codes using the CADNA approach. Graillat et al. [38, 39] studied the SA in multi precision, dynamical control of Newton’s method for multiple roots of polynomials and numerical validation of compensated summation algorithms with SA. Eberhart et al. [40] presented a high performance numerical validation using the SA. In recent years, the CESTAC method was applied to validate the results of several numerical methods such as estimating the value of interpolation polynomials [41], improper and definite integrals [42, 43], dynamical control on Gauss-Laguerre integration rule [44], evaluating the fuzzy definite integrals [45], numerical validation of the Sinc-collocation method to solve linear integral equations [46], finding optimal iteration of the power and inverse iteration methods, solving linear systems, fuzzy Newton-Cotes integration rules [47], finding the optimal parameter of the homotopy analysis method for solving integral equations [48], solving fuzzy integral equations [49], validating Taylor-collocation method for solving Volterra integral equations with discontinuous kernels [50] and load leveling problem [51].

The CESTAC method applies the SA in its computations and it is more useful than the FPA. Also, in order to implement the CESTAC method, the CADNA library (http://cadna.lip6.fr) must be used instead of the common mathematical packages such as Matlab or Maple. In this library we prepare the programs using C/C++ or FORTRAN and they are performed on the LINUX operating system. This method is able to eliminate the unessential iterations but the FPA does not have this ability.

2. Preleminaries

2.1 Mathematical osmosis model

The following model is presented to forecast the condensation of salt solutions in semi-penetrable covers in the converse osmosis model as [17]:

$y \frac{\partial \Phi}{\partial x}=\alpha \frac{\partial^{2} \Phi}{\partial^{2} y}$,    (4)

with $\alpha=\frac{D h}{v_{0}}$ and the boundary conditions are:

$\Phi(0, y)=\phi_{0}, \Phi(x, \infty)=\phi_{1}$,    (5)

and

$-D \frac{\partial \Phi}{\partial y}(x, 0)=q \Phi(x, 0)$,    (6)

where, space variables are demonstrated by x and y and $\Phi=\Phi(x, y)$ shows the condensation of salt solutions in semi-penetrable covers at point (x, y). q is the velocity of water flow in semi-permeable distribution, D is the salt diffusion in water, h is the distance from semi-permeable boundary to canal center and finally $\phi_{0}$ and $v_{0}$ are the concentration away from semi-permeable membranes and horizontal velocity at distance h from semi-boundary. Also, we should note that q, D, h and $v_{0}$ are constant values.

2.2 Gauss-Legendre integration rule

The quadrature integration formula (1) for $v(r)=1$ and $[\alpha, \beta]=[-1,1]$ is called the G-LIR as follows:

$Q=\int_{-1}^{1} S(r) d r \simeq Q_{J}=\sum_{\ell=0}^{J} w_{\ell} S\left(r_{\ell}\right)$,    (7)

where, the nodes $\left\{r_{\ell}\right\}_{\ell=0}^{n}$, are zeros of the Legendre polynomial.

$P_{J}(r)=\frac{1}{2^{J} J !} \frac{d^{J}}{d r^{J}}\left[\left(r^{2}-1\right)^{J}\right]$.

The polynomials $P_{J}(r)$ satisfy in the following recurrence relation:

$\left\{\begin{array}{l}P_{0}(r)=1, \\ P_{1}(r)=r, \\ J P_{J}(r)=(2 J-1) r P_{J-1}(r)-(J-1) P_{I-2}(r), J \geq 2.\end{array}\right.$    (8)

Theorem 1. The solution of G-LIR on [-1,1] with J points for polynomials of degree $\leq 2 J-1$ is exact.

Proof: see the papers [52, 53].

Theorem 2. If $S \in C^{2 J}[-1,1]$, the J-the order error term of the G-LIR is in the following form:

$\begin{aligned} E_{J}(S)=Q_{J}-Q &=\frac{2^{2 J+1}(J !)^{4}}{(2 J+1)((2 J) !)^{2}} \frac{S^{(2 J)}(\eta)}{(2 J) !} \\ &=e_{J} \frac{S^{(2 J)}(\eta)}{(2 J) !} \end{aligned}$,    (9)

where, $e_{J}=\frac{2^{2 J+1}(J !)^{4}}{(2 J+1)((2 J) !)^{2}}$ and $-1<\eta<1$.

Proof: see the references [52, 53].

Now, for finding the order of error term in Eq. (9) at first we define:

$M_{\hbar}=\max _{-1 \leq r \leq 1} \frac{\left|S^{(\hbar)}(r)\right|}{\hbar !}, \hbar \geq 0$    (10)

For a large class of infinitely differentiable functions S on [-1,1], we have $M_{\hbar}<\infty, \hbar \geq 0$. By combining Eqns. (9) and (10), we obtain:

$\left|E_{J}(S)\right| \leq e_{J} M_{2 J}, J \geq 1$    (11)

and the size of $e_{J}$ is essential in examining the rate of convergence. The term $e_{J}$ can be made more understandable by estimating it using Stirling’s formula [54],

$J ! \simeq e^{-J} J^{J} \sqrt{2 \pi J}$,    (12)

which is true in a relative error sense as $J \rightarrow \infty$. Then we obtain:

$e_{J} \simeq \frac{\pi}{4^{J}} a s J \rightarrow \infty$.    (13)

Now, we can write

$E_{J}(S)=\mathcal{O}\left(4^{-J}\right)$    (14)

The details of the G-LIR can be found in the papers [52, 53, 55, 56].

For integrals on other finite intervals, the following formula can be applied

$\int_{\alpha}^{\beta} S(t) d t=\left(\frac{\beta-\alpha}{2}\right) \int_{-1}^{1} S\left(\frac{\alpha+\beta+r(\beta-\alpha)}{2}\right) d r$,    (15)

reducing the integral to the standard interval [-1,1].

2.3 CESTAC method-CADNA library

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 the presented method based on the SA in discrete case are investigated in comparison with the traditional FPA [32, 34-37, 40].

Assume that some representable values are produced by computer and they are collected in set A. Then $W \in A$ can be produced for $w \in \mathbb{R}$ with $\mathcal{R}$ mantissa bits of the binary FPA in the following form:

$W=w-\chi 2^{E-\mathcal{R}} \xi$,    (16)

where, sign of w showed by x, missing segment of the mantissa presented by $2^{-\mathcal{R}} \xi$ and the binary exponent of the result characterized by E. Moreover, in single and double precisions $\mathcal{R}=24,53$ respectively.

Assume ξ is the casual variable that uniformly distributed on [-1,1]. After making perturbation on the final mantissa bit of w we will have (μ) and (σ) as mean and standard deviation for results of W which have important role in accuracy of W. Repeating this process J times for $W_{i}, i=1, \ldots, 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, \ldots, J$. For more consideration, the following algorithm [28] is presented where $\tau_{\delta}$ is the value of T distribution as the confidence interval is 1-δ with J-1 freedom degree.

Algorithm 1:

Step 1- Make the perturbation of the last bit of mantissa to produce J samples of W as $\Phi=\left\{W_{1}, W_{2}, \ldots, W_{J}\right\}$.

Step 2- Find $W_{a v e}=\frac{\sum_{i=1}^{J} W_{i}}{J}$.

Step 3- Compute $\sigma^{2}=\frac{\sum_{i=1}^{J}\left(W_{i}-W_{a v e}\right)^{2}}{J-1}$.

Step 4- Find the estimation of NCSDs of w and $W_{\text {ave}}$ applying $C_{W_{\text {ave}}, w}=\log _{10} \frac{\sqrt{J}\left|W_{\text {ave}}\right|}{\tau_{\delta} \sigma}$.

Step 5- Print W=@.0 if $W_{\text {ave}}=0$, or $C_{W_{a v e}, W} \leq 0$.

In general form, 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 as a new tool to implement the CESTAC method [31, 43].

There are important advantages to applying the CESTAC method and the CADNA library instead of traditional packages which are based on the FPA. In this method, instead of applying Eq. (2) we present a criterion independence of absolute error and tolerance value like ε [34, 35, 37]. Applying the CADNA library, we can find the optimal iteration, best approximation in the point of computational view and estimation of the accuracy of numerical methods. Moreover, the numerical instabilities can be identified. A sample program of the CADNA library is presented as:

# include <cadna.h>

cadna _ init(-1);

 main()

 {

double _ st Parameter;

do

{

Write the main program here;

printf(" %s ",strp(Parameter));

}

 while(Q[J]-Q[J-1]!=0);

 cadna _ end();

 }

The details of the SA, CESTAC method and CADNA library can be found in the paper [40].

3. Main Idea

In order to estimate the roots $r_{\ell}$ and weights $w_{\ell}, \ell=0,1, \ldots, J$ of G-LIR the following Newton-Raphson iteration formula [52] is applied.

$r_{\ell+1}=r_{\ell}-\frac{S\left(r_{\ell}\right)}{S^{\prime}\left(r_{\ell}\right)}=r_{\ell}-\frac{P_{J}\left(r_{\ell}\right)}{P_{I}^{\prime}\left(r_{\ell}\right)}$,    (17)

where, the following relation is used to find the first guess $r_{0}$ for the $\ell$-th root of $P_{J}$.

$r_{0, \ell}=\cos \left(\frac{\pi\left(\ell-\frac{1}{4}\right)}{J+\frac{1}{2}}\right)$,    (18)

and the recursive relation for derivative of $P_{J}$ is in the following form:

$P_{J}^{\prime}(r)=\frac{J}{r^{2}-1}\left(r P_{J}(r)-P_{J-1}(r)\right)$.    (19)

Also, the appropriate weights $w_{\ell}, \ell=0,1, \ldots, J$ is given by:

$w_{\ell}=\frac{2}{\left(1-r_{\ell}^{2}\right)\left[P_{I}^{\prime}\left(r_{\ell}\right)\right]^{2}}.$    (20)

Now, by substituting the obtained roots and weights in Eq. (7), the approximate solution can be estimated. For example, the nodes and weights for J=1,2,3,4,5 are presented in Table 1.

Table 1. The nodes and weights of the G-LIR

J

$\boldsymbol{r}_{\ell}$

$\boldsymbol{w}_{\boldsymbol{\ell}}$

1

±0.5773502692

1.0000000000

2

0.0000000000

±0.7745966692

0.8888888889

0.5555555556

3

±0.3399810436

±0.8611363116

0.6521451549

0.3478548451

4

±0.0000000000

±0.5384693101

±0.9061798459

0.5688888889

0.4786286705

0.2369268851

5

±0.2386191861

±0.6612093865

±0.9324695142

0.4679139346

0.3607615730

0.1713244924

Definition 1. [57] For $\rho_{1^{\prime}} \rho_{2} \in \mathbb{R}$, the NCSDs is estimated as:

(1) for $\rho_{1} \neq \rho_{2}$,

$C_{\rho_{1}, \rho_{2}}=\log _{10}\left|\frac{\rho_{1}+\rho_{2}}{2\left(\rho_{1}-\rho_{2}\right)}\right|=\log _{10}\left|\frac{\rho_{1}}{\rho_{1}-\rho_{2}}-\frac{1}{2}\right|$.    (21)

(2) $C_{\rho_{1}, \rho_{1}}=+\infty$.

Theorem 3. Let $S(r) \in C^{2 J+2}[-1,1]$ and $Q_{J}$ be the numerical solution of Q by using the G-LIR (7), then

$C_{Q_{I}, Q_{I+1}} \simeq C_{Q_{I}, Q}$.    (22)

Proof: Using Definition 1 we get:

$C_{Q_{n}, Q_{n+1}}=\log _{10}\left|\frac{Q_{n}}{Q_{n}-Q_{n+1}}-\frac{1}{2}\right|$

$=\log _{10}\left|\frac{Q_{n}}{Q_{n}-Q_{n+1}}\right|+\log _{10}\left|1-\frac{1}{2 Q_{n}}\left(Q_{n}-Q_{n+1}\right)\right|$

$=\log _{10}\left|\frac{Q_{n}}{Q_{n}-Q_{n+1}}\right|+\mathcal{O}\left(Q_{n}-Q_{n+1}\right)$

$=\log _{10}\left|\frac{Q_{n}}{\left(Q_{n}-Q\right)-\left(Q_{n+1}-Q\right)}\right|+\mathcal{O}\left[\left(Q_{n}-Q\right)-\left(Q_{n+1}-Q\right)\right]$

$=\log _{10}\left|\frac{Q_{n}}{\left(Q_{n}-Q\right)\left[1-\frac{Q_{n+1}-Q}{Q_{n}-Q}\right.}\right|+\mathcal{O}\left(E_{n}\right)+\mathcal{O}\left(E_{n+1}\right)$,   (23)

and using Eq. (14) we get:

$C_{Q_{J}, Q_{j+1}}=\log _{10}\left|\frac{Q_{J}}{Q_{I}-Q}\right|-\log _{10}\left|1-\frac{Q_{J+1}-Q}{Q_{I}-Q}\right|+\mathcal{O}\left(4^{-J}\right)$.    (24)

Also,

$\begin{aligned} C_{Q_{J}, Q} &=\log _{10}\left|\frac{Q_{J}}{Q_{J}-Q}-\frac{1}{2}\right| \\ &=\log _{10}\left|\frac{Q_{J}}{Q_{J}-Q}\right|+\mathcal{O}\left(E_{J}\right) \\ &=\log _{10}\left|\frac{Q_{J}}{Q_{J}-Q}\right|+\mathcal{O}\left(4^{-J}\right) \end{aligned}$.    (25)

Applying Eqns. (24) and (25) we have:

$C_{Q_{J}, Q_{J+1}}=C_{Q_{J}, Q}-\log _{10}\left|1-\frac{Q_{J+1}-Q}{Q_{J}-Q}\right|+\mathcal{O}\left(4^{-J}\right)$.

Consequently, based on Eq. (9) we have $\frac{Q_{J+1}-Q}{Q_{J}-Q}=\frac{\mathcal{O}((2 J+3)(2 J+1))}{\mathcal{O}\left(2^{J}(J+1)^{2}\right)}$ and we can find that for J enough large this term tends to zero and $O\left(4^{-} J\right)<<1$ thus we get:

$C_{Q_{J}, Q_{J+1}} \simeq C_{Q_{J}, Q}$.

Theorem 3 shows that, for enough large value of J, the NCSDs between $Q_{J}$ and $Q_{J+1}$ are almost equal to the NCSDs between Q and $Q_{J}$. Therefore, if we use the CESTAC method, the computations of $Q_{J}$’s sequence will be stopped when $Q_{J+1}-Q_{J}=@ .0$. In this case, the values J+1 and $Q_{J+1}$ are the optimal iteration and approximation of the G-LIR respectively which are shown by by $J_{o p t}$ and $Q_{J o p t}$.

4. Numerical Illustrations

In this section, the G-LIR is validated by applying the CESTAC method to estimate some definite integrals. The approximate solution is produced until $\left|Q_{J+1}-Q_{J}\right| \neq @ .0$. If this condition equals to @.0, the presented algorithm will be stopped. It means that the values $Q_{J}$ and $Q_{J+1}$ are equal stochastically and the number of iteration will be the optimal iteration of G-LIR and shows by $J_{o p t}$. In some of examples, the numerical results are obtained based on two arithmetics; the FPA and the SA. By comparing the results, it can be shown the SA is more applicable and reliable than the FPA.

In order to perform the proposed algorithm, the CADNA library is presented. All codes are provided by C++ and implement on Linux with double precision.

Algorithm 2:

Step 1- Let J=0 and $Q_{0}=0$.

Step 2- Do the following steps while $\left|Q_{J+1}-Q_{J}\right| \neq @ .0$

{

Step 2-1- Construct $P_{J}(r)$ and $P_{J}^{\prime}(r)$ by using Eqns. (8) and (19).

Step 2-2- Calculate the first guess $r_{0, \ell}$ based on Eq. (18).

Step 2-3- Produce roots $r_{\ell}$ by using Newton-Raphson iteration formula (17).

Step 2-4- Compute the weights $w_{\ell}$ by applying Eq. (20).

Step 2-5- Calculate the approximate solution $Q_{J+1}$ by means of Eq. (4).

Step 2-6- Print $Q_{J+1}$, $\left|Q_{J+1}-Q_{J}\right|$ and $\left|Q_{J+1}-Q\right|$.

Step 2-7-$J=J+1$.

 }

Table 2. Numerical results of Example 1

J

$Q_{J+1}$

$\left|Q_{J+1}-Q_{J}\right|$

$\left|Q_{J+1}-Q\right|$

0

0.84200591794852E-003

0.84200591794852E-003

0.450530640557518E+000

1

0.748058884321424E+000

0.747216878403476E+000

0.29668623784595E+000

2

0.86409516325659E+000

0.11603627893516E+000

0.41272251678112E+000

3

0.50123463565037E+000

0.36286052760621E+000

0.4986198917490E-001

4

0.30198192671927E+000

0.19925270893109E+000

0.14939071975619E+000

5

0.37081909274923E+000

0.68837166029966E-001

0.80553553726227E-001

6

0.49496778905600E+000

0.12414869630676E+000

0.4359514258053E-001

7

0.493059958255680E+000

0.190783080032E-002

0.41687311780213E-001

8

0.44280099467176E+000

0.5025896358391E-001

0.8571651803703E-002

9

0.43418145795721E+000

0.8619536714545E-002

0.1719118851824E-001

10

0.45248079593407E+000

0.182993379768E-001

0.11081494586E-002

11

0.457637149034123E+000

0.515635310004E-002

0.6264502558657E-002

12

0.45135438319577E+000

0.628276583834E-002

0.1826327968E-004

13

0.44922625662707E+000

0.212812656870E-002

0.214638984838E-002

14

0.45136405761603E+000

0.213780098895E-002

0.8588859428E-005

15

0.452078451409432E+000

0.71439379339E-003

0.705804933966E-003

16

0.45134806353781E+000

0.73038787161E-003

0.2458293764E-004

17

0.45115062558794E+000

0.19743794987E-003

0.22202088752E-003

18

0.451397794645076E+000

0.24716905713E-003

0.2514816960E-004

19

0.451438116342501E+000

0.40321697425E-004

0.65469867034E-004

20

0.451357582657320E+000

0.80533685181E-004

0.15063818146E-004

21

0.451355225400271E+000

0.235725704E-005

0.1742107519E-004

22

0.45137959305346E+000

0.2436765318E-004

0.6946577994E-005

23

0.45137653311710E+000

0.305993635E-005

0.388664163E-005

24

0.451369996731362E+000

0.6536385742E-005

0.26497441042E-005

25

0.451372061091343E+000

0.2064359980E-005

0.585384123E-006

26

0.451373494695151E+000

0.143360380E-005

0.848219684E-006

27

0.45137262420657E+000

0.870488577E-006

0.2226889E-007

28

0.451372423846007E+000

0.20036056E-006

0.222629458E-006

29

0.451372705102738E+000

0.28125673E-006

0.58627271E-007

30

0.45137269022487E+000

0.1487786E-007

0.4374940E-007

31

0.45137262003629E+000

0.7018858E-007

0.2643917E-007

32

0.45137264237891E+000

0.2234262E-007

0.409655E-008

33

0.451372654286528E+000

0.1190761E-007

0.78110617E-008

34

0.451372645260884E+000

0.9025644E-008

0.121458E-008

35

0.451372644906413E+000

0.35447E-009

0.156905E-008

36

0.45137264724825E+000

0.234184E-008

0.772790E-009

37

0.451372646613299E+000

0.63495E-009

0.13783E-009

38

0.451372646246873E+000

0.366426E-009

0.22859E-009

39

0.451372646517676E+000

0.270803E-009

0.42209E-010

40

0.45137264651480E+000

0.2875E-011

0.3933E-010

41

0.451372646452266E+000

0.6253E-010

0.2320E-010

42

0.451372646474531E+000

0.22265E-010

0.934E-012

43

0.451372646481058E+000

0.6526E-011

0.5592E-011

44

0.451372646473655E+000

0.7403E-011

0.181E-011

45

0.45137264647488E+000

0.123E-011

0.580E-012

46

0.45137264647608E+000

0.1194E-011

0.613E-012

47

0.45137264647535E+000

0.72E-012

0.10E-012

48

0.45137264647537E+000

0.1E-013

0.93E-013

49

0.451372646475524E+000

0.15E-012

0.57E-013

50

0.45137264647546E+000

0.6E-013

@.0

51

0.451372646475454E+000

@.0

@.0

Table 3. Applying the CESTAC method for solving Example 2

J

$Q_{J+1}$

$\left|Q_{J+1}-Q_{J}\right|$

$\left|Q_{J+1}-Q\right|$

0

0.558607885129995E+000

0.558607885129995E+000

0.80340631273229E-001

1

0.476468795302816E+000

0.82139089827178E-001

0.179845855394E-002

2

0.478282274977675E+000

0.181395447393E-002

0.1549591998E-004

$\vdots$

$\vdots$

$\vdots$

$\vdots$

5

0.47826725385638E+000

0.19885E-009

0.37E-012

6

0.478267253856765E+000

0.7E-012

@.0

7

0.478267253856765E+000

@.0

@.0

Example 1. According to the Ref. [17], solution of osmosis model (4) is in the following form:

$\Phi(x, 0)=\frac{3^{-\frac{1}{3}} q \phi_{0}}{D I}\left(\frac{D h}{v_{0}}\right)^{\frac{1}{3}} x^{\frac{1}{3}}+\phi_{0}$,    (26)

where, $I=\int_{0}^{\infty} e^{-v^{3}} v d v$ and Eq. (26) is used to predict the condensation of salt solutions in semi-penetrable covers in the converse osmosis model. Solving the improper integral I by using analytical methods is difficult, so we need to solve this integral numerically.

In order to evaluate I, we consider interval [0, m] where m is a number enough large such that $\left|\int_{m}^{\infty} e^{-v^{3}} v d v\right|=@ .0$ which means this value has not any significant digits. By choosing m=10, we have:

$Q=\int_{0}^{10} e^{-v^{3}} v d v=0.4513726464754668 \ldots$,

and using Eq. (15), we can transform this integral to the following form:

$25 \int_{-1}^{1} e^{-125(1+z)^{3}}(1+z) d z$.

Based on the presented results in Table 2, the optimal step of iteration is $J_{o p t}=51$ and the optimal result is $Q_51=0.451372646475454E+000$ with optimal error 0.6E-013. Thus the optimal solution of the osmosis model is:

$\frac{\Phi(x, 0)}{\phi_{0}}=1.5361171751\left(\frac{q}{D}\right)\left(\frac{D h}{v_{0}}\right)^{\frac{1}{3}} x^{\frac{1}{3}}+1$,

where, it is a mathematical formula to forecast the condensation of salt solutions in semi-penetrable covers in the converse osmosis model.

Example 2. Let [53]

$Q=\int_{-1}^{1} r^{2} \cos r d r=0.47826725385676605 \ldots$.    (27)

According to obtained results which are presented in Table 3 the optimal iteration of G-LIR to approximate integral (27) is $J_{o p t}=7$ and the approximate solution is $Q_{I+1}=0.478267253856765 E+000$.

Example 3. In this example, the integral

$Q=\int_{0}^{1} \frac{d r}{1+r}=0.6931471805599453 \ldots$,    (28)

is evaluated [53]. In order to apply the G-LIR the interval [0,1] is transformed to [-1,1] as:

$\int_{-1}^{1} \frac{d z}{z+3}$.    (29)

In Table 4, the optimal step of G-LIR, approximate solution, absolute error and difference between two successive approximations are presented. In this case, $J_{o p t}=9$ is the optimal step of this rule. According to Table 4 the optimal value of integral is $Q_{J+1}=0.693147180559945 E+000$.

Example 4. Consider the following integral [53]

$Q=\int_{0}^{2} \frac{r^{2}+2 r+1}{r^{2}+2} d r=2.42310142981207 \ldots$.    (30)

By using the transformation formula (15), the interval [0,2] should be changed to [-1,1] as:

$\int_{-1}^{1} \frac{(z+2)^{2}}{z^{2}+2 z+3} d z$.    (31)

In Table 5, the numerical results and the accuracy of this method are shown. The optimal step of G-LIR is $J_{o p t}=13$ and the optimal approximation is $Q_{J+1}=0.242310142981206 E+001$.

Table 4. Numerical results of Example 3

J

$Q_{J+1}$

$\left|Q_{J+1}-Q_{J}\right|$

$\left|Q_{J+1}-Q\right|$

0

0.692307692307692E+000

0.692307692307692E+000

0.839488252253E-003

1

0.693121693121692E+000

0.814000814000E-003

0.25487438252E-004

2

0.693146417445482E+000

0.2472432378E-004

0.763114462E-006

$\vdots$

$\vdots$

$\vdots$

$\vdots$

6

0.693147180559355E+000

0.1934E-010

0.589E-012

7

0.693147180559927E+000

0.571E-012

0.1E-013

8

0.693147180559944E+000

0.1E-013

@.0

9

0.693147180559945E+000

@.0

@.0

Table 5. Applying the CESTAC method for solving Example 4

$J$

$Q_{J+1}$

$\left|Q_{J+1}-Q_{J}\right|$

$\left|Q_{J+1}-Q\right|$

0

0.240909090909090E+001

0.240909090909090E+001

0.1401052072116E-001

1

0.242255892255892E+001

0.1346801346801E-001

0.542507253148E-003

2

0.242322960686387E+001

0.67068430495E-003

0.12817705180E-003

$\vdots$

$\vdots$

$\vdots$

$\vdots$

11

0.242310142981208E+001

0.1E-012

0.1E-013

12

0.242310142981206E+001

0.1E-013

@0

13

0.242310142981206E+001

@0

@0

Example 5. Consider the following integral [53].

$Q=\int_{-1}^{1} \frac{d r}{\log \left(r^{2}\right)}$.

We note that this integral is not convergent on [-1,1]. Table 6, shows the results obtained by the FPA for $\varepsilon=10^{-5}$ and these results are not correct. Also, the number of iteration for different values of ε are presented in Table 7. The results of Algorithm 2 are demonstrated in Table 8. In the last iteration of Table 8 we have $Q_{J+1}=@ .0$. Therefore CADNA library declares that the numerical results are not guaranteed and shows the report of numerical instabilities. The FPA has not these capabilities in comparison with the SA.

Table 6. Numerical results of the FPA for Example 5 with $\varepsilon=10^{-5}$

J

$Q_{J+1}$

$\left|Q_{J+1}-Q_{J}\right|$

0

-1.820478

1.820478

1

-2.175128

0.354650

2

-2.931223

0.756095

$\vdots$

$\vdots$

$\vdots$

28

-6.721456

0.081742

29

-6.771934

0.050477

30

-6.848311

0.076377

$\vdots$

$\vdots$

$\vdots$

66

-8.338342

0.034800

67

-8.362105

0.023763

68

-8.395876

0.033771

$\vdots$

$\vdots$

$\vdots$

Table 7. Number of iteration for different values of ε in Example 5

ε

$10^{-5}$

$10^{-3}$

$10^{-1}$

0.5

1

n

without stopping

without stopping

13

1

1

Table 8. Numerical results for Example 5

J

$Q_{J+1}$

$\left|Q_{J+1}-Q_{J}\right|$

0

-0.182047E+001

0.182047E+001

1

-.2175127E+001

0.354648E+000

2

-0.293122E+001

0.75609E+000

$\vdots$

$\vdots$

$\vdots$

43

-1.000001E+001

0.100000E+001

44

0.999976E+000

0.199996E+001

45

@.0

@.0

Example 6. Let [53]

$Q=\int_{-1}^{1} \tan \left(r^{2}-r\right) d r$.

This integral is not convergent in [-1,1]. In Table 9, we show the numerical results which are obtained by FPA for $\varepsilon=10^{-5}$. Also, in Table 11, number of iteration for different values of ε are presented. The results of CADNA library are reported in Table 10. In this Table 10, the approximate solution for J=47 is equal to the informatical zero @.0. Obtaining numerical noise (@.0) for J=47 means that the result associated to J=47 has no correct digit.

Table 9. Applying the FPA for solving Example 6 with $\varepsilon=10^{-5}$

J

$Q_{J+1}$

$\left|Q_{J+1}-Q_{J}\right|$

0

1.039207

1.039207

1

2.697160

1.657953

2

-10.773429

13.470588

$\vdots$

$\vdots$

$\vdots$

48

-1.398668

37.830112

49

-0.175755

1.222913

50

0.505504

0.681259

$\vdots$

$\vdots$

$\vdots$

98

-9.852456

12.666245

99

-1.113649

8.738807

100

-0.061589

1.052060

$\vdots$

$\vdots$

$\vdots$

Table 10. Numerical results of Example 6

J

$Q_{J+1}$

$\left|Q_{J+1}-Q_{J}\right|$

0

0.103920E+001

0.103920E+001

1

0.269715E+001

0.165795E+001

2

-0.10773E+002

0.13470E+002

$\vdots$

$\vdots$

$\vdots$

46

-0.5E+003

0.5E+003

47

@.0

@.0

Table 11. Number of iteration for different values of ε in Example 6

ε

$10^{-5}$

$10^{-3}$

$10^{-1}$

0.5

1

1.5

n

without stopping

without stopping

without stopping

without stopping

5

0

5. Conclusion

We have considered the CESTAC method based on the SA in discrete case to find the optimal iteration, the optimal approximation of the G-LIR. Also, by using this method the unnecessary iterations can be eliminated. They are the main novelties of the CESTAC method. The advantages of the SA in comparison with the FPA were discussed. The validation of the CESTAC method to implement the G-LIR was guaranteed by proving Theorem 3. Some examples specially the mathematical model of osmosis model were evaluated by the proposed algorithm via the CADNA library. Also, in order to show the flexibility of the CESTAC method in some of examples we compared the numerical results in both cases the SA and the FPA. Moreover, we found the number of iterations for different values of ε. The proposed scheme can be developed to evaluate the optimal value of high dimensional G-LIR, and other kinds of quadrature rules such as Gauss-Hermit, Gauss-Rado and Gauss-Lobatto integration rules.

  References

[1] Noeiaghdam, S. (2019). A novel technique to solve the modified epidemiological model of computer viruses. SeMA Journal, 76(1): 97-108. http://dx.doi.org/10.1007/s40324-018-0163-3 

[2] Noeiaghdam, S., Suleman, M., Budak, H. (2018). Solving a modified nonlinear epidemiological model of computer viruses by homotopy analysis method. Mathematical Sciences, 12(3): 211-222. http://dx.doi.org/10.1007/s40096-018-0261-5

[3] Ozturk, Y., Gulsu, M. (2015). Numerical solution of a modified epidemiological model for computer viruses. Applied Mathematical Modelling, 39(23-24): 7600-7610. https://doi.org/10.1016/j.apm.2015.03.023 

[4] Günerhan, H., Dutta, H., Weaver, M.A., Adel, W., (2020). Analysis of a fractional HIV model with Caputo and constant proportional Caputo operators Chaos. Solitons & Fractals, 139: 110053. https://doi.org/10.1016/j.chaos.2020.110053 

[5] Naik, P.A., Zu, J., Owolabi, K.M. (2020). Global dynamics of a fractional order model for the transmission of HIV epidemic with optimal control. Chaos Solitons Fractals, 138: 109826. https://doi.org/10.1016/j.chaos.2020.109826

[6] Naik, P.A., Yavuz, M., Zu, J. (2020). The Role of Prostitution on HIV Transmission with Memory: A modeling approach. Alexandria Engineering Journal, 59(4): 2513-2531. https://doi.org/10.1016/j.aej.2020.04.016 

[7] Noeiaghdam, S., Ghiasi, E.K. (2019). An efficient method to solve the mathematical model of HIV infection for CD8+ T-cells. International Journal of Mathematical Modelling & Computations, 9(4): 267-281. arXiv:1907.01106.

[8] Guerrero, F., Santonja, F.J., Villanueva, R.J. (2013). Solving a model for the evolution of smoking habit in Spain with homotopy analysis method. Nonlinear Analysis: Real World Applications, 14(1): 549-558. https://doi.org/10.1016/j.nonrwa.2012.07.015 

[9] Sun, M., Tian, L.X., Jia, Q. (2009). Adaptive control and synchronization of a four-dimensional energy resources system with unknown parameters. Chaos Solitons Fract., 39(4): 1943-1949. http://dx.doi.org/10.1016/j.chaos.2007.06.117

[10] Sun, M., Tao, Y., Wang, X., Tian, L. (2011). The model reference control for the four-dimensional energy supply-demand system. Applied Mathematical Modelling, 35(10): 5165-5172. http://dx.doi.org/10.1016/j.apm.2011.04.016

[11] Guerrero Sánchez, Y., Sabir, Z., Günerhan, H., Baskonus, H.M. (2020). Analytical and approximate solutions of a novel nervous stomach mathematical model. Discrete Dynamics in Nature and Society, 2020: 5063271. https://doi.org/10.1155/2020/5063271 

[12] Srivastava, H.M., Günerhan, H. (2019). Analytical and approximate solutions of fractional order susceptible-infected-recovered epidemic model of childhood disease. Mathematical Methods in the Applied Sciences, 42(3): 935-941. https://doi.org/10.1002/mma.5396

[13] Altaf Khan, M., Ullah, S., Farooq, M. (2018). A new fractional model for tuberculosis with relapse via Atangana–Baleanu derivative. Chaos, Solitons and Fractals, 116: 227-238. https://doi.org/10.1016/j.chaos.2018.09.039 

[14] Garud, R.M., Kore, S.V., Kore, V.S., Kulkarni, G.S. (2011). A short review on process and applications of reverse osmosis. Universal Journal of Environmental Research and Technology, 1: 233-238.

[15] Wimalawansa, S.J. (2013). Purification of contaminated water with reverse osmosis: Effective solution of providing clean water for human needs in developing countries. International Journal of Emerging Technology and Advanced Engineering, 3(12): 75-89.

[16] Sanjaya, F., Mungkasi, S. (2017). A simple but accurate explicit finite difference method for the advectiondiffusion equation. Journal of Physics: Conference Series, 909: 012038. https://doi.org/10.1088/1742-6596/909/1/012038 

[17] Fulford, G.R., Broadbridge, P. (2002). Industrial Mathematics: Case Studies in the Diffusion of Heat and Matter. Cambridge University Press, Cambridge.

[18] Maure, O.P., Mungkasi, S. (2019). Application of numerical integration in solving a reverse osmosis model. AIP Conference Proceedings, 2202(1): 020043. https://doi.org/10.1063/1.5141656 

[19] Babolian, E., Masjed-Jamei, M., Eslahchi, M.R. (2005). On numerical improvement of Gauss-Legendre quadrature rules. Appl. Math. Comput., 160(3): 779-789. https://doi.org/10.1016/j.amc.2003.11.031

[20] Lether, F.G. (1978). On the construction of Gauss-Legendre quadrature rules. J. Comput. Appl. Math., 4(1): 47-52. https://doi.org/10.1016/0771-050X(78)90019-0

[21] Petras, K. (1999). On the computation of the Gauss-Legendre quadrature formula with a given precision. J. Comput. Appl. Math., 112(1-2): 253-267. https://doi.org/10.1016/S0377-0427(99)00225-3

[22] Rathod, H.T., Venkatesudu, B., Nagaraja, K.V., Islam, M.I. (2007). Gauss Legendre-Gauss-Jacobi quadrature rules over a tetrahedral region. Appl. Math. Comput., 190(1): 186-194. https://doi.org/10.1016/j.amc.2007.01.014

[23] Deloff, A. (2008). Gauss-Legendre and Chebyshev quadratures for singular integrals. Comput. Phys. Comm., 179(12): 908-914. https://doi.org/10.1016/j.cpc.2008.07.008

[24] Verlinden, P. (1997). Acceleration of Gauss-Legendre quadrature for an integrand with an endpoint singularity. J. Comput. Appl. Math., 77(1-2): 277-287. https://doi.org/10.1016/S0377-0427(96)00131-8

[25] Wang, L.L., Guo, B. (2009). Interpolation approximations based on Gauss-Lobatto-Legendre-Birkhoff quadrature. J. Approx. Theory, 161(1): 142-173. https://doi.org/10.1016/j.jat.2008.08.016

[26] Elliott, D., Johnston, P.R. (2008). Gauss-Legendre quadrature for the evaluation of integrals involving the Hankel function. J. Comput. Appl. Math., 211(1): 23-35. https://doi.org/10.1016/j.cam.2006.11.002

[27] Vignes, J. (1993). A stochastic arithmetic for reliable scientific computation. Math. Comput. Simulation, 35(3): 233-261. https://doi.org/10.1016/0378-4754(93)90003-D

[28] Vignes, J. (2004). Discrete Stochastic Arithmetic for validating results of numerical software. Numerical Algorithms, 37(1-4): 377-390. https://doi.org/10.1023/B:NUMA.0000049483.75679.ce

[29] Chesneaux, J.M. (1994). The equality relations in scientific computing. Numer. Algorithms, 7: 129-143. https://doi.org/10.1007/BF02140678

[30] Chesneaux, J.M. (1992). Stochastic arithmetic properties. IMACS Comput. Appl. Math., 81-91.

[31] Chesneaux, J.M. (1990). CADNA: An ADA tool for round-off errors analysis and for numerical debugging. Proceedings of the Congress ADA in Aerospace, Barcelone.

[32] Lamotte, J.L., Chesneaux, J.M., Jézéquel, F. (2010). CADNA _-C: A version of CADNA for use with C or C++ programs. Computer Physics Communications, 181(11): 1925-1926. https://doi.org/10.1016/j.cpc.2010.07.006

[33] Alt, R., Vignes, J. (1996). Validation of results of collocation methods for ODEs with the CADNA library. Appl. Numer. Math., 21(2): 119-139. https://doi.org/10.1016/0168-9274(96)00004-9

[34] Jézéquel, F., Chesneaux, J.M. (2004). Computation of an infinite integral using Romberg’s method. Num. Algo., 36(3): 265-283. https://doi.org/10.1023/B:NUMA.0000040066.63826.46

[35] Jézéquel, F. (2006). A dynamical strategy for approximation methods. Comptes Rendus Mécanique, 334(6): 362-367. https://doi.org/10.1016/j.crme.2006.04.005

[36] Jézéquel, F., Chesneaux, J.M. (2008). CADNA: a library for estimating round-off error propagation. Computer Physics Communications, 178(12): 933-955.

[37] Scott, N.S., Jézéquel, F., Denis, C., Chesneaux, J.M. (2007). Numerical ’health check’ for scientific codes: the CADNA approach. Computer Physics Communications, 176(8): 507-521. https://doi.org/10.1007/978-3-642-14403-5_35 

[38] Graillat, S., Jézéquel, F., Wang, S., Zhu, Y. (2011). Stochastic arithmetic in multi precision. Math. Comput. Sci., 5: 359-375. https://doi.org/10.1007/s11786-011-0103-4 

[39] Graillat, S., Jézéquel, F., Picot, R. (2015). Numerical validation of compensated summation algorithms with stochastic arithmetic. Electron. Notes Theor. Comput. Sci., 317: 55-69. https://doi.org/10.1016/j.entcs.2015.10.007 

[40] Eberhart, P., Brajard, J., Fortin, P., Jézéquel, F. (2015). High performance numerical validation using stochastic arithmetic. Reliable Computing, 21: 35-52.

[41] Abbasbandy, S., Fariborzi Araghi, M.A. (2004). The use of the stochastic arithmetic to estimate the value of interpolation polynomial with optimal degree. Appl. Numer. Math., 50(3-4): 279-290. https://doi.org/10.1016/j.apnum.2004.01.003

[42] Abbasbandy, S., Fariborzi Araghi, M.A. (2002). Numerical solution of improper integrals with valid implementation. Math. Comput. Appl., 7(1): 83-91. https://doi.org/10.3390/mca7010083

[43] Abbasbandy, S., Fariborzi Araghi, M.A. (2005). A stochastic scheme for solving definite integrals. Appl. Numer. Math., 55(2): 125-136. https://doi.org/10.1016/j.apnum.2004.11.007

[44] Fariborzi Araghi, M.A., Noeighdam, S. (2016). Dynamical control of computations using the Gauss-Laguerre integration rule by applying the CADNA library. Advances and Applications in Mathematical Sciences, 16(1): 1-18.

[45] Fariborzi Araghi, M.A., Noeighdam, S. (2017). A valid scheme to evaluate fuzzy definite integrals by applying the CADNA library. International Journal of Fuzzy System Applications, 6(4): 1-20. https://doi.org/10.4018/IJFSA.2017100101

[46] Fariborzi Araghi, M.A., Noeighdam, S. (2019). Valid implementation of the Sinc-collocation method to solve linear integral equations by the CADNA library. Journal of Mathematical Modeling, 7(1): 63-84. https://doi.org/10.22124/jmm.2018.11608.1200

[47] Noeiaghdam, S., Fariborzi Araghi, M.A. (2017). Finding optimal step of fuzzy Newton-Cotes integration rules by using the CESTAC method. Journal of Fuzzy Set Valued Analysis, 2017(2): 62-85. https://doi.org/10.5899/2017/JFSVA-00383 

[48] Noeiaghdam, S., Fariborzi Araghi, M.A., Abbasbandy, S. (2019). Finding optimal convergence control parameter in the homotopy analysis method to solve integral equations based on the stochastic arithmetic. Numer Algor, 81(1): 237-267. https://doi.org/10.1007/s11075-018-0546-7 

[49] Noeiaghdam, S., Fariborzi Araghi, M.A., Abbasbandy, S. (2020). Valid implementation of Sinc-collocation method to solve the fuzzy Fredholm integral equation. Journal of Computational and Applied Mathematics, 370: 112632. https://doi.org/10.1016/j.cam.2019.112632

[50] Noeiaghdam, S., Sidorov, D., Sizikov, V., Sidorov, N. (2020). Control of accuracy on Taylor-collocation method to solve the weakly regular Volterra integral equations of the first kind by using the CESTAC method. Applied and Computational Mathematics, 19: 87-105. arXiv:1811.09802.

[51] Noeiaghdam, S., Sidorov, D., Muftahov, I., Zhukov, A.V. (2019). Control of accuracy on taylor-collocation method for load leveling problem. The Bulletin of Irkutsk State University. Series Mathematics, 30: 59-72. https://doi.org/10.26516/1997-7670.2019.30.59 

[52] Abramowitz, M., Stegun, I.A. (1972). Handbook of Mathematical Functions, Dover, New York.

[53] Jain, M.K., Iyengar, S.R.K. (2009). Numerical Methods, New Age International (p) Limited, Daryaganj, New Delhi - 110002.

[54] Atkinson, K. (1993). Elementary Numerical Analysis. second ed., John Wiley.

[55] Burden, R.L., Fairs, J.D. (1985). Numerical Analysis. third ed., PWS, Boston.

[56] Davis, P.J., Rabinowitz, P. (1984). Methods of Numerical Integration. Academic Press, New York, 2nd edition.

[57] Chesneaux, J.M., Jézéquel, F. (1998). Dynamical control of computations using the Trapezoidal and Simpson’s rules. Journal of Universal Computer Science, 4(1): 2-10. https://doi.org/10.3217/jucs-004-01-0002