Abstract
Implementation of interval arithmetic in complex problems has been hampered by the tedious programming exercise needed to develop a particular implementation. In order to improve productivity, the use of interval mathematics is demonstrated using the computing platform INTLAB that allows for the development of interval-arithmetic-based programs more efficiently than with previous interval-arithmetic libraries. An interval-Newton Generalized-Bisection (IN/GB) method is developed in this platform and applied to determine the solutions of selected nonlinear problems. Cases 1 and 2 demonstrate the effectiveness of the implementation applied to traditional polynomial problems. Case 3 demonstrates the robustness of the implementation in the case of multiple specific volume solutions. Case 4 exemplifies the robustness and effectiveness of the implementation in the determination of multiple critical points for a mixture of methane and hydrogen sulfide. The examples demonstrate the effectiveness of the method by finding all existing roots with mathematical certainty.
1. Introduction
There are a large number of problems that require the computation of stationary points and roots of equations. These are found in the fields of optimization [1], economics and finance [2, 3], thermodynamics [4, 5], applied mathematics [6], and similar contributions over the last thirty years. The application of interval arithmetic involves from algebraic equations with known solutions to more complicated systems representing physical phenomena. Among the earlier problems, second- and third-order polynomial problems serve to illustrate the effectiveness of the implementation. Similarly, more complicated multidimensional problems serve to illustrate stability and robustness of the implementation. In particular, the determination of critical points is of interest. Critical point computations are a well-known highly nonlinear problem that has been studied for a long time. The literature is endowed with some worthy analyses of the problem. Michelsen and Heidemann [7] discuss the calculation of critical points from cubic two-constant equations of state. They numerically determine the critical points of mixtures solving the highly nonlinear critical point equations. Hoteit et al. [8] claim an efficient algorithm based on bisection, secant, and inverse quadratic programming methods where their method is apparently faster but incapable of handling the presence of multiple critical points without restarting the program after each critical point determination. Nichita and Gomez [9] discuss the use of the tunneling global optimization method to determine all the of global minima where they concentrate on the determination of the critical points of mixtures. They use temperature and molar volume as primary variables. Their method lacks the strength of being able to stop when all critical points have been found, and this follows from the fact that at the beginning of the process the number of critical points to be determined needs to be specified. Michelsen [10] develops some earlier ideas that deal with the computation of phase envelopes with use of finite differences. Michelsen [11] discusses the isothermal flash separation and critical point computations with the use of tangent plane analysis. Michelsen and Kistenmacher [12] describe the disadvantage of using composition-dependent binary interaction coefficients in computations. Michelsen and Heidemann [13] present the calculation of tricritical points using tangent-plane distance analysis with some extra algebra in the mathematical development. The results are generated using the Peng-Robinson and Soave-Redlich-Kwong (SRK) equations of state. Michelsen [14] illustrates the determination of the critical points and phase boundaries using the real-arithmetic Newton-Raphson method with temperature and pressure as variables; in his analysis rather close initial estimates to the solution are needed for convergence.
There are also other more general contributions such as that of Cai et al. [15] who developed expressions for the spinodal criterion, critical criterion, and various stability tests for systems containing one discrete component and one polydisperse polymer. Lindvig et al. [16] propose the EFV (Entropic-FV) model for predicting the miscibility behavior of paints. Cismondi and Michelsen [17] use a Newton procedure for the calculation of critical lines, critical endpoints, and three-phase lines for binary mixtures. Carstensen and Petković [18] propose the use of a hybrid method combining normal (Nourein's method) and interval arithmetic to arrive to the solution of polynomial equations efficiently. Maranas and Floudas [19] propose an approach based on convex lower bounding and domain partitioning to achieve convergence to all solutions within a domain. Stradi et al. [4] discuss the computation of critical points using the INTBIS library [20].
This paper presents the solution of nonlinear problems using interval arithmetic in INTLAB (INTerval LABoratory) [21] by using the interval-Newton Generalized-Bisection method (IN/GB). The advantage of this approach is that INTLAB is developed over MATLAB [22], a well-known and powerful computational package widely used in mathematics and engineering applications [23–25]. Although, with the cost of interpretation overhead [26], INTLAB facilitates interval computations by using the MATLAB user interface and debugging tools, this makes production times smaller and results generation faster [27]. Section 2 of the paper provides a description of the problems to be solved in the paper, Section 3 proposes the methodology to solve the problems using the IN/GB method, Section 4 presents the results and discussion comparing the IN/GB method with the simple bisection method, and Section 5 highlights the conclusions derived from this study.
2. Problem Description
This section starts with two examples to demonstrate the use and performance of the IN/GB method. The first deals with two second-order equations where the solution is to be found over a relatively large domain. The second illustrates a particularly important case of a repeated root that causes the first derivative to be zero, and consequently it would traditionally create problems during the process of solution for the real-arithmetic Newton-Raphson method [23]. The third and fourth examples are related to the behavior of fluids in the two phase and critical regions.
Consider the general case where all of the roots to the nonlinear problem have to be found within a given domain in the absence of an initial guess.
Case 1. The first problem is a set of two second-order equations as follows: The domain for the search is the same for both variables: . Three solutions are found for these two equations with mathematical certainty in a single run of the IN/GB method implemented over the INTLAB program.
Case 2. The second problem is a third-order polynomial with three real roots, where one of the roots is repeated:
This is an important example because the derivative of the function is zero at the repeated root. The traditional real-arithmetic Newton-Raphson method [23] would fail if we needed to apply the method near the root where the derivative is close to zero. The domain of the search is .
The interval Newton is expected to have problems because of its reliance on the interval derivative of the function to find a solution, but the combination with the bisection method should solve the problem.
Case 3. The third problem deals with the determination of volume solutions for the Peng-Robinson equation of state [28] for carbon dioxide at a given temperature and pressure, which is given by the following expression that related pressure, temperature, and volume:
This can also be written in the form
The necessary conditions for this problem are the following: ,
,
,
,
,
,
and the interval for the volume search is given by .
This is a classic chemical engineering problem where, for a given equation of state, there may be more than one molar volume that satisfies the equation of state at a given temperature and pressure. In particular, if the conditions are such that two phases coexist, like in this case, then there will be three solutions. The traditional interpretation is that the solutions represent the liquid molar volume (smallest volume solution) and the vapor molar volume (largest volume solution). The middle volume solution is an unstable solution, and thus not observable.
Case 4. This is a highly nonlinear problem subject to analysis by different research groups for a number of years [4, 7, 14, 17]. In this case, there is a mixture of methane and hydrogen sulfide where the problem is to determine all of the critical points for a given composition and domain ranges for temperature and molar volume using the Redlich-Kwong equation of state [28]. The problem was solved in a previous development with the use of the INTBIS library and some extensive programming.
The domains of the variables under study are indicated as follows: ,
,
,
.
The criticality conditions are written as follows:
The mass balance for a mixture of constant composition is given by
Meanings are given in the appendix.
3. Methodology
We employed the interval mathematics platform INTLAB to implement the interval-Newton Generalized-Bisection (IN/GB) method for the solution of nonlinear problems. The advantage of this procedure is that INTLAB is user friendly and runs over MATLAB; this makes INTLAB a very promising platform by putting the tools of MATLAB at the disposal of the INTLAB programmer. In the past, one of the major problems with the use of interval mathematics was the use of libraries that required extensive programming and lacked utilities for fast debugging. As a result, large time investments needed to be made in programming and debugging prior to the implementation and generation of first results, a situation that obviously derives in making interval arithmetic less attractive to a larger scientific audience. However, with INTLAB, the situation is different: program sizes are much smaller, debugging is more efficient, and results are generated more rapidly.
3.1. The Interval-Newton Generalized-Bisection Method
An interval-Newton Generalized-Bisection method (IN/GB) is presented to demonstrate the use of the INTLAB platform for nonlinear problems of different levels of complexity.
The interval IN/GB method is a powerful computational approach that allows for the computation of all of the roots of a nonlinear problem without the need for initial guesses for the variables [1, 29, 30].
The IN/GB method requires only that we specify the domain of the variables of interest. The program tests a sequence of enclosures determined through the IN/GB method. If there is a root, then the IN/GB method determines a very thin slice in which the root lies. On the other hand, if there is no root, then the algorithm also determines with mathematical certainty the absence of solutions.
The most important component part of the IN/GB method is the Newton method in interval arithmetic. This is written as follows:
This equation is of the form , and consequently suitable for solution with simultaneous linear equations solvers.
In this equation, capital letters are intervals and small case letters are real numbers. is the interval extension of the Jacobian matrix, which is determined by evaluating the Jacobian matrix, , using the interval domain , rather than the real variable , where is generally the midpoint of . This system of equations is solved using the Gauss-Seidel method [1], or a similar algorithm to determine , where is the image generated by the application of the IN/GB method. If and intersect, then the process is repeated, noting that in successive applications it is the intersection of and () that serves as initial interval. The process continues until a sufficiently thin interval around the solution is derived. Other implementations, when sufficiently closed to the solution, switch to the real-arithmetic Newton-Raphson method to accelerate convergence [20].
If is a subset of , then there is a single root in the domain of interest [29]. Similarly, if does not intersect , then there is no root in the domain, and a new subdomain is tested for the presence of roots.
It is important to indicate at this point that there are several problems not generally mentioned in the literature that may occur when using the IN/GB method. If the interval is too large, then the image will contain and be larger than the original interval , generating the problem where the image is larger than the original interval and no convergence is achieved with the interval-Newton portion of the method alone. The other problem is where the image is equal to ; in this case the search will not advance, and no solution will be achieved within the permissible execution time. In both cases, the action taken is to proceed with the partition of the box, as shown in Table 1, and restart the search for each subdomain generated. The partition is done with the Bisection method [23] over each dimension of the domain followed by a range test applied to all subdomains generated.
Application of the Bisection method generates two intervals from each interval that is bisected. Consequently, in order to explore all of the possibilities, the combinations of these intervals are needed to determine the existence and values of all possible solutions.
The existence of a root is ascertained by evaluating the original function over the interval, or subdomain, of interest, using interval arithmetic [30]. If the resulting interval contains zero, , then there may exist at least one root in that domain. This process of testing for the existence of a root in the domain of interest is called range test. If the evaluation of the function over the interval of interest, , results in an interval image that does not contain zero, then there is a mathematical guarantee that there is no root within the interval .
For example, if three intervals are bisected, then six subintervals are generated, and eight combinations of them would need to be tested first in the range test to ascertain whether there may be a root in the subdomain, and second in the IN/GB method search for the solution.
Generally, multiple applications of the Bisection method, part of the IN/GB method, are needed when large variable domains are used. This situation is particularly true when zeros are contained in diagonal elements of the interval Jacobian expression, . The presence of these zeros generates interval images, , in the IN/GB method with infinite spans for some or all variables. Consequently, bisection is needed to divide and eliminate those subdomains where infinite spans are generated. These computational problems are managed efficiently in INTLAB. In INTLAB, the programming environment allows for processing warnings identified by the flags infinite (inf) and not-a-number (NaN) without an error that would stop the program execution. The flag inf would identify quantities that exceed real number representation such as and also results generated by the division of a number by zero: . The flag NaN would identify quantities that are mathematically undefined such as or intervals that do not intersect. Consequently, if these flags occur, then bisection is applied because they are indicative of problems with the interval-Newton portion of the algorithm. This capacity to handle numerically undefined quantities saves extensive programming time and debugging effort invested on the solution of a particular problem.
The search for a solution can only end with one of two outcomes, either a root is found within an interval or no root is found with mathematical certainty. The main steps of the algorithm are represented in Figure 8.
4. Results and Discussion
4.1. Case 1
Table 2 presents the results obtained using the IN/GB and the Bisection methods applied to Case 1. The IN/GB method determined all of the roots and searched the entire domain with 685 subdivisions. The Bisection method with interval arithmetic was applied, and convergence to all roots was achieved with 820 subdivisions of the original interval. It is important to mention that, before applying either the IN/GB method or the Bisection method, a range test is applied to the subdomain under consideration to determine whether a root may be found within the subdomain. Three solutions are found for these two equations. The algorithm developed found three solutions without the need for initial guesses over a fairly large domain. This case is a baseline case where the solutions can be calculated also with the real-arithmetic Newton method, but then at least three different initial guesses would need to be provided.
Figure 1 shows the evolution of the error with the subinterval number for the Bisection method; Figure 2 shows the evolution of the error with the subinterval number for the IN/GB method. The error limit is represented with a horizontal line at 1 * 10−6, and the trends are approximated with a third-order polynomial in a semilog plot for the Bisection and the IN/GB methods. The interval-Newton part of the IN/GB method may generate some very large errors particularly when the derivative, or the Jacobian matrix, contains a zero in the diagonal-element intervals. Consequently, in order to keep the scale to a legible size, large errors described as inf and nondefined quantities described as NaN appear with an error value of 10 on the graph’s -axis. For both the bisection method and the IN/GB method, the solutions are found toward the end of the search of subdomains.
4.2. Case 2
Table 3 presents the results for Case 2, where the IN/GB method determines all of the roots to this third-order equation using 17723 subdivisions to complete the search over the entire domain. The bisection method requires 17981 subdivisions to search the entire domain and find all of the roots. It is important to notice that the first root is found much faster by the IN/GB method by a factor of more than a hundred times over the bisection method. However, the convergence speed is smaller at the second root where both algorithms need closely the same number of subdivisions. For this simple example, the roots are evident and can be determined by inspection. However, if using the real-arithmetic Newton-Raphson method, there would be no immediate way to determine that one of the roots is repeated, and this situation would require additional time to test additional initial guesses or execute a polynomial deflation procedure with the roots already found.
Figure 3 shows the evolution of the error with the subinterval number for the Bisection method. Figure 4 shows the evolution of the error with the subinterval number for the IN/GB method.
In the case of the Bisection method, it is particularly different the convergence profile, compared to Case 1; both cases have been approximated with a third-order polynomial determined with a least-squares fit. The IN/GB method presents some very large errors at the beginning of the search, and a log-log plot is used to better depict this behavior. Large error values are constrained to a value of 10 for graphical clarity, and they are indicative of inf and NaN results; in these cases, the algorithm switches from the interval-Newton to the Bisection method. In this case, the roots are found by IN/GB method both at the beginning and at the end of the process of testing all generated subintervals as indicated by the values below 1 * 10−6.
4.3. Case 3
Table 4 presents the results for Case 3. This is a particularly interesting example because it exposes some of the difficulties in interval computations. In this example, the interval-Newton method does not result in an interval, , with the desired characteristics: that be of smaller size than the original interval, , and that and intersect. During the search process with the interval-Newton, this portion of the IN/GB method generates NaN or inf messages, and the program switches to the Bisection method. This point is very important because it clearly demonstrates that if for some reason the interval-Newton portion of the algorithm is incapable of dealing with the problem then the other half of the algorithm (i.e., the Bisection method) takes over and provides a solution, making the IN/GB a very efficient method.
Once an interval is bisected, a range test on each subdomain created is performed, and if a root is detected (i.e., ), then that interval is stacked for further processing. This process will continue until either each subdomain is eliminated for not containing zero (i.e., ) or a root is found within the subdomain with a prescribed error tolerance. In this example, the IN/GB method finds the solutions and scans the entire domain by dividing the domain in 4332 subsets. Consequently, the graph in Figure 5 for the Bisection and the IN/GB methods is the same.
4.4. Case 4
Table 5 presents the results for Case 4. This is a more elaborate example because for a multicomponent mixture there is no procedure to determine a priori whether there will be none, one, or more than one critical points. The IN/GB method is much faster (Figure 6) than the Bisection method (Figure 7), and the order for finding the solutions differs between the methods. The IN/GB finds the solutions and scans the entire domain by dividing the domain in 345733 subsets, while the Bisection method finds the solutions and scans the entire domain by dividing the domain in 865970 subsets.
It is illustrative to plot the error in the IN/GB method versus the number of subdivisions of the initial domain. There are some very large errors, particularly at the beginning, which are cut to a value of 10 units, in order to preserve a reasonable scaling in the figures. As it was mentioned previously, for those cases with very large errors (inf) or not-a-number warnings (NaN), the algorithm proceeded to bisecting those domains in order to continue with the search. The relative error used to determine solutions was of 1 * 10−6, which is indicated in the graphs with a flat line.
The Bisection method converges linearly in a semilog plot toward the solutions which are found at the end of the process, as it is expected. However, the IN/GB finds solutions prior to the end of the process and with fewer subdivisions of the original interval.
5. Conclusions
This paper has introduced the use of INTLAB in computations for the determination of several solutions without the need for initial guesses and with a single run of the algorithm. The IN/GB method was implemented in INTLAB with four examples. The first two examples were algebraic equations that serve as a reference to verify the effectiveness of the algorithm and describe some important findings of the IN/GB method. The third and fourth cases are related to the computation of multiple solutions that represent properties of substances and fluids. Convergence to a solution is faster with the IN/GB method in general, but there are exceptions to this rule particularly when a singular Jacobian matrix may occur.
Appendix
where the following are considered. is the total number of components. is the second-order derivative of the Helmholtz free energy with respect to the number of moles of species and . is the change in total number of moles of species . This is used in the context of a Taylor series expansion of the Helmholtz free energy in terms of composition. is the ideal gas constant. is the absolute temperature. is the total number of moles. is the number of moles of species , analogous meaning for other subindices. is the parameter used in the equations for the computation of critical points; it is defined as is the mole fraction of component , analogous meaning for other subindices. is the average energy parameter of the equation of state for the mixture of components. This is computed using standard van der Waals mixing rules: is the energy of interaction parameter between species and . is the energy parameter of species . The meaning is the same for other subindices. It is computed using the following correlation: is the accentric factor of species . is the critical temperature of species . is the critical pressure of species . is the accentric factor of species . is the binary interaction parameter between species and . is the parameter used in the equations for the computation of critical points; it is defined as is the parameter used in the equations for the computation of critical points; it is defined as is the parameter used in the equations for the computation of critical points; it is defined as is the average van der Waals mole volume of the equation of state. This is computed using standard van der Waals mixing rules: is the van der Waals mole volume of species . This is computed as follows: for the Redlich-Kwong-Soave equation of state, and for the Peng-Robinson equation of state. is the parameter used in the equations for the computation of critical points; it is defined as is the parameter used in the equations for the computation of critical points; it is defined as are the auxiliary functions used in the computation of critical points; they are defined as is the dimensionless volume; it is defined as are the parameters of the equation of state. Their values are calculated with the following formulas:
where for the Redlich-Kwong-Soave equation of state and for the Peng-Robinson equation of state.
Acknowledgment
The author wishes to acknowledge the helpful comments and review provided by Professor Emmanuel Haven at the University of Leicester, UK.