Solving interval systems of equations obtained during the numerical solution of boundary value problems

Numerical solution of boundary value problems modelled by differential or integral equations is reduced to solving linear system of equations. Even in the systems of equations without intervals, the solution does not have to be unique, e.g., in the case of a singular (or at least ill-conditioned) matrix. Such problems are even more inconvenient in the interval linear systems of equations. In this paper, the various methods of solving such systems are tested. These systems have been generated during the numerical solution of boundary value problems modelled by Navier–Lamé equations.

of algebraic equations is not yet the solution of the problem. To obtain unique solution of linear system of equations, the matrix determinant should be different from zero, and the matrix should be well conditioned. It is possible, that linear system of equations obtained as the result of approximation of differential (integral) equations may not meet conditions determining its unique solution.
In recent years, boundary problems are modelled by integral equations using interval numbers [for example (Piasecka Belkhayat and Korczak 2014;Zalewski et al. 2009)]. Such modelling of problems is more reliable, because it gives opportunity to include measurement errors in modelling the shape of boundary and boundary conditions. However, problem of interval modelling the shape of boundary is still open. Many authors [for example (Dessombz et al. 2001;Muhanna et al. 2005;Piasecka Belkhayat and Korczak 2014;Zalewski et al. 2009)] defined only boundary conditions and input parameters as intervals. In process of modelling boundary problems, intervals appears as measurement errors in data definition. Numerical solution of such problems is reduced to solving interval systems of equations. Solving these systems is more complicated than solving the ones with precise elements. It requires application of interval arithmetic (Hayes 2003;Krämer 2013;Moore et al. 2009).
Methods of solving systems of equations were examined in case of rounding errors (for example, generated by hardware decoding of floating-points operations) (Krämer et al. 2006;Rump 2014). The authors applied interval arithmetic to solve such systems. Rounding errors are described by very narrow intervals. They are significantly smaller then measurement errors and may be ignored in modelling of boundary problems. In this case, direct application of methods described in Krämer et al. (2006) and Rump (2014) gives wide width of solution intervals. It is useless in practical point of view.
In this paper, we apply three selected numerical methods for solving interval linear systems of equations, i.e. the Gauss elimination method (Neumaier 1990), the LU decomposition method (Goldsztejn and Chambert 2007) and the Jacobi iterative method (Markov 1999). We also apply directed interval arithmetic (Markov 1995) to minimize width of interval solutions. Because of lack of implementation of directed interval arithmetic in well-known libraries, such as C-XSC or Intlab, we implement own C++ class based on double-precision floating-point type.
In this paper, the basic concepts of interval numbers are presented. Selected numerical methods for solving interval linear systems of equations are described. The aim of this study is to test the effectiveness of mentioned methods for solving interval systems of equations obtained during solving boundary value problems. To solve such problems, classical boundary element method (BEM) (Brebbia et al. 1984) or its modification based on parametric integral systems of equations (PIES) (Zieniuk 2013) can be used. Analysis of the results is performed in terms of their accuracy, depending on the size and conditioning of the matrix of the system of equations.

Basic concepts of interval numbers
Interval numbers were determined by the interval, whose left endpoint had to be smaller than the right one. These numbers are called proper numbers, classical intervals or just intervals and they are presented as follows Moore et al. 2009): To operate on these numbers special interval arithmetic was created and it is generally defined by the formula (where • ∈ {+, −, ·, /}) (Hayes 2003): (2) Our objective is to estimate the range of a function. Formula (2) provides the exact description of the range of each arithmetic operation over given intervals. During the development of various methods using interval numbers (Dawood 2011; Hansen and Walster 2004;Neumaier 1990), it turns out that it is impossible to obtain opposite and inverse element of such numbers. Thus problem occurs in solving simple interval equation as well as in interval systems of equations. For this purpose, an extension of classical intervals by improper numbers has been proposed (Markov 1995). These are numbers which the upper endpoint is smaller than the lower one. Such an extension was called directed or extended intervals (Popova 2001). In addition, two operators were added to get inverse and opposite element of the given number. They allow to obtain additional arithmetic operators of subtraction and division: (3) x/y, x/y for x < 0, y < 0 x/y, x/y for x > 0, y < 0 x/y, x/y for x < 0, y > 0 x/y, x/y for x ∈ 0, y < 0 x/y, x/y for x ∈ 0, y > 0 (4)

Methods of solving interval systems of equations
General matrix form of interval system of equations can be written as follows (Neumaier 1990;Shary 2012): where A is a matrix of interval coefficients, x is a vector of interval solutions, and b is a vector of intervals. There are many methods of solving such system of equations. However, we should apply directed intervals, due to the way of modelling described in Chapter 4. Therefore in this paper, we use the LU method presented in Goldsztejn and Chambert (2007), the Jacobi method discussed by Markov (1999) and the Gauss elimination method (Neumaier 1990) as one of the most commonly used methods. Selected methods, tested and fully analysed, were presented below.

Gauss elimination method
The best-known and widely used method in the real numbers domain is the Gauss elimination method. A major step in solving interval system of equations using this method (Neumaier 1990) is to transform the matrix coefficients from the formula (5) to the following form: ⎡ To achieve the matrix A and the vector b, the following formulas are used: where i = 1, 2, . . . , n and j, k = i + 1, i + 2, . . . , n.
With an upper triangular matrix (6) it is easy to determine the value of the solution vector. Starting from the last row, elements of the matrix A are equal to the value in corresponding row of the vector b. The only difference between the presented method and the method with the real numbers is to use the interval arithmetic.

The LU decomposition method
An extension of the Gauss elimination method is the LU decomposition method. In this method, the matrix A of the system of equation (5) can be written as the product of two triangular matrices: lower L and upper U (Goldsztejn and Chambert 2007): The solution of interval system of equations is reduced to solving two interval systems of equations with the triangular matrices: The elements of the matrix L and U are calculated alternately. First, the values in the row of the U matrix are calculated, then the column of the L matrix. The corresponding elements of the matrix for all i ∈ {1, 2, . . . , n} can be calculated: Analysing the Eq. (11), it can be seen that the method will not work if the value of any element on the diagonal of the matrix U will be equal to zero (Goldsztejn and Chambert 2007).

The Jacobi method
Another method of solving linear systems of equations is the iterative Jacobi method (Markov 1999). In this method, the interval system of Eq. (5) is solved by the following formula: where x i is the next approximation of the solution vector, A is the matrix of interval coefficients, b is a vector of intervals, T and T −1 are matrices defined as follows: Substituting (13) and (14) into (12) we can find solution of the system (5): Calculation are performed until the difference between calculated value and the value from previous iteration is small enough.

Analysis of solutions of interval linear system of equations obtained during numerical solving of boundary value problems
To solve boundary value problems, we can apply boundary element method (Brebbia et al. 1984) based on boundary integral equations (BIE) and similarly parametric integral equations system (Zieniuk 2013;Zieniuk et al. 2012Zieniuk et al. , 2014Kużelewski and Zieniuk 2014) which is a modification of classical BIE. Numerical solution of the problems using mentioned methods is reduced to solve linear system of equations. The system is defined on boundary of solving problem. However, first of all we should define the problem. It means that the shape of boundary should be presented as an interval. Such intervals include measurement errors ε which are associated with defining the boundary (Fig. 1). They are not included in classical modelling (use of real numbers). In practice, to define this problem the boundary should be presented as some band with a width ε, which includes measurement error. It is assumed that the precise boundary is located within the band (bold line in the Fig. 1 labelled as b). The bounds of interval band are the inner edge a and the outer edge c. Solution of a boundary problem with a shape of boundary defined as interval is reduced to solving interval system of equations. In this paper, modelling of the problem described by Navier-Lamé equations was considered. The aim of the research is to examine the effectiveness of methods of solving interval system of equations. Such system can be obtained by MEB (or PIES). Process of modelling of boundary problem is carried out two times, separately for the inner and outer edges. As result of this strategy, two linear systems of equations are obtained (16). The first of these relates to the inner edge (a on the Fig. 1) and the second to the outer edge (c on the Fig. 1    Interval system of equations is obtained combining elements matrices and vectors in (16). We decided to examine the efficiency of methods described in Sect. 3. To solve such a system, the classical and directed interval arithmetic were used. Finally, obtained solutions were compared.

Differences between traditional and directed intervals arithmetic
The main objective of traditional interval computations is to find enclosure [all possible values of x from Eq. (5)]. Directed interval arithmetic is focused on minimizing of solution error and it gives narrower intervals than classic one. Solution of the system of two equations using classical interval arithmetic (Hansen and Walster 2004) gives significantly wider intervals than using directed intervals (Markov 1999). Therefore, it was decided to check the case of linear system of 12 equations. For this purpose, the boundary value problem defined in the triangular area was considered. The system of equations generated during the numerical solution of the problem has been solved by the LU decomposition method. The results of application classical interval arithmetic and directed one are presented in Table 1.
It turns out that for system of 12 equations interval width of results for classical intervals is much greater than for directed ones, as well as in case of 2 equations. It is also noted that the middle value of classical interval solution is more shifted from the real number (b in Fig. 1) than directed one. It is connected with the fact that traditional intervals gives all possible solutions, which is usually useless in practice. Additionally, it is more accurate to use directed intervals because modelling problem in a way presented in Sect. 4 may produce improper intervals. Therefore, we decided to focus only on application of directed interval arithmetic.

Tests for systems with a small number of equations
To test the effectiveness of three selected methods of solving systems of equations described in Sect. 3, the problem defined in the triangular area with the generated system of 24 equations was examined. First to solve the system, the Jacobi iterative method was applied. Unfortu- nately, the method does not always converge to the correct result. It also happened in this case. Because it often happens that the Jacobi method fails, we decided to examine properties of generated system. For this purpose, the condition number has been calculated. It turned out that this is quite high value, much different from the 1. Therefore, to improve the condition of the matrix, we decided to make the following preconditioning of the system (5) (Neumaier 1990): where mid([x, x]) = (x + x)/2 is the middle value of interval. Therefore, the equation is multiplied by the inverse matrix of middle values. After the implementation of such transformation, the system has been solved by the Jacobi method. Because it is the iterative method which results does not always converge to the correct solution, it was decided to test the LU decomposition method and the Gauss elimination method. In case of solving the preconditioned system of equations, the results obtained by these methods are very close to the ones from the Jacobi method. However, the direct application (without preconditioning) of the LU decomposition or the Gauss elimination method slightly decreases the solution interval. This phenomenon is observed only at systems with a small number of equations. To verify the correctness and accuracy of the results obtained by different methods, relative error of the solutions is presented in (Fig. 2).
In this case, relative error is specified by absolute value of interval solution to precise one with application of degenerate intervals (Markov 1995). The exact value became a value of a precise solution (b in the Fig. 1), therefore, the formula of relative error was described in the following form: where: x interval solution of system equations, x 0 precise solution. Figure 2 shows, that after the preconditioning, the relative error of all considered methods is very similar. However, direct application of the LU decomposition or the Gauss elimination method to the generated system reduces slightly the error value. The correctness of these solutions was checked additionally by multiplying them by a matrix of coefficients A and then the comparison with

Tests for linear systems with a large number of equations
In the case of system of 24 equations, the differences between the methods are quite hardly visible. Therefore, it was decided to examine boundary value problem which required solution of the system of 96 equations. It turned out that in case of such system, conditioning of the matrix is deteriorating. Application of the Gauss elimination method gave slightly worse results than the LU decomposition method. In the Jacobi method, the convergence problem occurred. Therefore, the further analysis focuses on the LU decomposition method. Figure 3 shows relative error of selected solutions calculated by (18). Solutions were obtained by the LU method and preconditioned LU method (17). It turned out that preconditioning influences the improvement of results. Significant difference between selected maximum values of the relative errors is seen in Fig. 3. In addition, in case of the LU* method, improvements of the middle value of error can be noted. This is closer to the expected zero than in the case of the LU method. Another disadvantage of direct application of the LU method is that interval solution not always contains the precise solution. This problem does not appear in preconditioned LU method. Therefore, despite slightly worse value in the case of systems of 24 equations, we should apply preconditioned methods.

Influence of measurement error size ε into the solution
Until now, systems of equations have been created by taking relatively small measurement error (the value ε in Fig. 1). The value of ε = 0.001, and maximum length of the boundary of triangular area was 1. Therefore, it was decided to examine the change of solutions when the value of ε is changed. Figure 4 shows the relationship between the lower and upper endpoint of the interval solution, and the precise solution for different values of measurement error ε.
It can be seen, that for each element of the solution vector, the precise solution (grey bar) is always between the ends of the interval (black and white bars). In addition, it can be seen that the change of the width of selected errors affects the intervals width of the obtained solutions. Respectively for the value ε = 0.1, a quite large differences of the interval solution were observed. Sometimes the values were change from negative to positive. In case of the ε = 0.01 bars are closer, but still difference between some values can be clearly seen. Only in the case of ε = 0.001 difference is almost invisible.

Conclusions
This paper presents the effectiveness of the methods used to solve interval systems of equations generated during the numerical solution of boundary value problems. Comparison of three selected methods of solving such systems was performed using directed interval arithmetic. The Gauss elimination and the LU decomposition methods work well with systems of 24 equations. However, in case of 96 equations it may happen that the interval does not contain a precise solution. The reason for this behaviour of methods is ill-conditioned matrix obtained from PIES. Therefore, it seems reasonable to use preconditioned methods. After transformation, the effectiveness of both methods improved, whereas the middle of the interval relative error was closer to zero. In case of solving ill-conditioned linear systems of 96 equations, the LU method allows to get a little bit narrower intervals than the Gauss elimination method. In case of the Jacobi method, the main problem turns out to be the convergence of the method. Preconditioning increases the likelihood of convergence significantly, but does not guarantee that each problem can be solved. As a result, it was found, that the best method for solving interval systems of equations is preconditioned LU decomposition method. The main objective of this study was try to get a correct solution from practical point of view.
The intervals are generated by measurement errors which appear during the modelling of boundary value problem. The intervals appeared in the process of data definition. This paper is a necessary step to our further work on uncertainty in solving boundary value problems. For further work, the application of fuzzy sets (Reiser et al. 2013) can be considered. In this case, interval solutions should be analysed as suitable alpha-cuts.