Interval Arithmetic for Nonlinear Problem Solving

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 ofmultiple critical points for amixture ofmethane and hydrogen sulfide.The examples demonstrate the effectiveness of the method by finding all existing roots with mathematical certainty.


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] International Journal of Engineering Mathematics present the calculation of tricritical points using tangentplane 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][24][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.

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 () = 0 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:  ∈ [−100, 100],  ∈ [−100, 100].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  ∈ [−100, 100].
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: 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.

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.

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, 0 ⊂ (), 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  1000 and also results generated by the division of a number by zero: 1.0/0.0.The flag NaN would identify quantities that are mathematically undefined such as 0.0/0.0 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. 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 realarithmetic Newton method, but then at least three different initial guesses would need to be provided.

Results and Discussion
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.3 presents the results for Case 2, where the IN/GB method determines all of the roots to this thirdorder 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.

Case 2. Table
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.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.

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.

4. 1 .Figure 1 :
Figure 1: Evolution of the relative error (-axis) with the subinterval number for the Bisection method for Case 1.

Figure 2 :Figure 3 :
Figure 2: Evolution of the relative error (-axis) with the subinterval number for the IN/GB method for Case 1.

Table 1 :
Interval-Newton Generalized-Bisection method results and actions.
N: IN/GB method, B: Bisection method.