A fast numerical algorithm for finding all real solutions to a system of N nonlinear equations in a finite domain

Numerical Algorithms https://doi.org/10.1007/s11075-024-01908-7 RESEARCH A fast numerical algorithm for finding all real solutions to a system of N nonlinear equations in a finite domain Fernando Chueca-Díez1·Alfonso M. Gañán-Calvo2,3 Received: 4 January 2024 / Accepted: 30 July 2024 © The Author(s) 2024 Abstract A highly recurrent traditional bottleneck in applied mathematics, for which the most popular codes (Mathematica, Matlab, and Python as examples) do not offer a solution, is to find all the real solutions of a system of nnonlinear equations in a certain finite domain of the n-dimensional space of variables. We present two similar algorithms of minimum length and computational weight to solve this problem, in which one resembles a graphical tool of edge detection in an image extended to ndimensions. To do this, we discretize the n-dimensional space sector in which the solutions are sought. Once the discretized hypersurfaces (edges) defined by each nonlinear equation of the n-dimensional system have been identified in a single, simultaneous step, the coincidence of the hypersurfaces in each n-dimensional tile or cell containing at least one solution marks the approximate locations of all the hyperpoints that constitute the solutions. This makes the final Newton-Raphson step rapidly convergent to all the existent solutions in the predefined space sector with the desired degree of accuracy. Keywords Nonlinear systems ·Algebraic equations ·All-solutions ·Domain discretization 1 Introduction Numerical root finding methods such as Newton-Raphson, Broyden’s [2], or Levenberg-Marquardt [9] have proven countless times an invaluable tool in the STEM BAlfonso M. Gañán-Calvo [email protected] Fernando Chueca-Díez [email protected] 1University of Bristol, Queen’s Building, University Walk, Bristol BS8 1TR, UK 2Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos, ETSI, Universidad de Sevilla, Camino de los Descubrimientos, Sevilla 41092, Spain 3Laboratory of Engineering for Energy and Environmental Sustainability, Universidad de Sevilla, Camino de los Descubrimientos, Sevilla 41092, Spain 123 Numerical Algorithms fields. It has permitted the user arrive at a viable solution of a system of equations with relative little computational complexity where finding their analytical solutions becomes, in most cases, an impossible task. A common feature of these numerical methods is that they require an initial set of values which the user must choose to begin the iterative process. The choice of these initial values can, not only greatly effect the number of iterations required to arrive at an approximate solution (thereby incrementing the computational effort), but will dictate towards which solution the process will converge or if it will converge at all. Sometimes an estimate of the solution is known, either from past experience or from trial and error, but in general this is not the case. Moreover, the system of equations often have a large number of equally important solutions, requiring efficient methods to ensure that all solutions are found in the absence of analytic formulae (e.g. the solution of general polynomial equations of order higher than four, as Galois and Abel independently showed). Without a failproof system that guarantees the collection of all solutions in the domain, a search process can become whimsical and unreliable, making the use of an algorithm that can reliably detect all solutions a necessity [15]. An example of such a necessity is to find all the equilibrium solutions of non-linear systems of ODEs. This is critical to assess the global phase map of the general solutions of the ODEs. Another typical example is the solution of the complex frequencies and wavelengths of normal modes of dynamical systems defined by dispersion relations, for example in fluid problems with free surfaces. Yet another fundamental example where our method could be particularly useful may be the efficient solution of a system of partial differential equations, such as the steady incompressible Navier- Stokes equations in a given domain, defining the original variables (i.e., velocity and pressure) and their derivatives as unknowns, and relating the latter to the original variables by additional algebraic equations (e.g., using appropriate finite differences or pseudo-spectral approximation expressions). This latter use of our method will be explored in a future work. Therefore, the goal is to obtain initial approximate values for all solutions within the selected domain while minimizing the iterative process of the root finding method, with minimal computational and implementation effort. Other methods exist to attain this objective [1,6–8,11,13–15]. However, our proposed algorithm has proven to have a substantially lower complexity level in our own experience, resulting in easier implementation and much lower computational steps. 1.1 Basic definitions A nonlinear system of equations F(x)may be defined as follows without any loss of generality, F(x)= ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ f1(x) f2(x) . . . fn(x) ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ (1) 123 Numerical Algorithms with x∈S=[a1,b1]∗[a2,b2]...[an,bn]⊂Rnand f1,f2,..., fnbeing nonlinear scalar functions such that f:S→R.Ifx∗is a solution to the system of equations we can say, F(x∗)= ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ f1(x∗)=0 f2(x∗)=0 . . . fn(x∗)=0 ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ (2) where x={x1, ..., xn}. Thus, the problem entails discovering the roots of every function in the equation system within a desired sector of the ndimensional space of variables. Each function defines a (n−1)-dimensional manifold within the n-dimensional domain. Every intersection of all these manifolds (hyper-surfaces) identifies a solution to the equation system. We take advantage of edge detection ideas in graphics algorithms by extending the concept of edge in a planar image (i.e. two dimensions) to a hyper-surface in ndimensions. We will see that the hyper-surface corresponding to the roots of each equation of system (2) lies in the subset of cells where the function fichanges its sign. 2 The proposed algorithm and initial identification schemes We begin by assuming that the n-dimensional space sector S, of arbitrary shape, has been reduced to an n-hypercube by the appropriate mapping of the original sector (see Fig. 1). It is assumed that the mapping into a hypercube conserves the original space connectivity of Sand involves a much simpler operation than the root finding. In most cases, the mapping simply entails a linear re-scaling of an original hyperprism of arbitrary volume corresponding to the original variables. We then divide the hypercube into (N−1)ncells, where Nrepresents the level of discretization of each dimension, and generate the discretized set Acomprising Nnpoints (i.e., Acontains n×Nnscalar values). The optimal memory allocation of the set Aallows for the subsequent one-step evaluation of the nequations in all Nnpoints to be completed in a highly efficient manner. This can be achieved by utilizing current powerful algebraic manipulation libraries such as Cupy, which are available in Python and take advantage of GPU computing. By computing fi(A)at each vertex of the Nndiscretised hyper-cells, we obtain nvalues at each vertex. These nvalues can be considered the new coordinates in a mapped space. In a strict sense, the function Fperforms a mapping of the set Ainto a distorted set Bgenerated by the transformation B=F(A), which can be as complex as the complexity of F implies. From a computational perspective, the discrete mapping B=F(A)represents the sole function call made prior to the identification of the points in Awhich are closer to those x∗such that F(x∗)=0. In this paper, we present two different schemes, designated as A and B, for the identification of those hyper-points in Awhich are closer to solutions x∗. Both methods utilize the idea of computing the local changes of sign of fi(x)around each point in A, and identifying those points showing changes in all i-directions. By leveraging the highly efficient matrix operations offered by current codes as CUDA, along with 123 Numerical Algorithms mapping discretisation Fig. 1 Initial operation scheme some ingenuity presented here as schemes A and B, this operation can be completed at an extremely rapid pace. In this initial identification, all solutions are contained into a reduced number of cells. Both schemes, as stated before, involve minimal computational steps, making them extremely efficient and easy to implement to find an initial solution approach. Subsequently, the vertex with lower coordinate values of each solution-containing cell is used as the seed in a classical Newton-Raphson scheme, which can rapidly converge to the final desired precision. Even though several solutions may appear repeated (depending on the level of discretization and the local gradient of each function fi), this scheme is advantageous as it requires significantly fewer computations compared to similar methods which use iterative processes to hone in on the solution. Obviously, our methodology can be nested within the solution-containing cells if Fpresented fractal-like distributions of roots or involved singularities. 2.1 Method A Method Ais a product-based method. We first obtain nsets Bi=F(Ai=1,...,n), where Aiis Adisplaced (or perturbed) one distance between successive points in the positive i-direction (see Fig. 2). Then, we generate the set Ci=sign1(BBi), where function sign1() converts positive and negative values to 0 and 1 respectively, and is the inner product. In simple terms, Ciindicates the occurrence of a change in sign around each point in each positive i-direction of every scalar function fj(for j=1, ..., n). Consequently, the vector Ciis composed of n×Nnscalar values. Next, we define the 123 Numerical Algorithms Fig. 2 Graphical representation of scheme A. The arrows indicate the directions of the ncoordinates (here, n=3) along which the discretised domain A is shifted by one distance between successive points set G(A)=n i=1Ci, which has n×Nnscalar components as well, corresponding to the changes of sign in any direction of each scalar function fjat each of the the Nn points. In formal terms, the compound mapping G(A)=n i=1sign1(F(A)F(Ai))indicates whether a point in Ais closest to any root of each scalar fjfunctions. In other words, the values which are non-zero of the mapping G(A)identify the points of A which are closest to the hyper-surface defined by fj(x)=0. Among those points, those whose transformation G(A)have all scalar components equal to non-zero are closest to a solution of F(x)=0. Computationally, the automatic identification of those points that lay more closely to the solutions is performed by defining the scalar function (A)(A)as the inner product of the nscalar components of G(A)at each point of A. Thus, those points in Asuch that (A)= 0 are the points used as seeds in the final Newton-Raphson step. A detailed example of the steps taken in this scheme is given in Appendix 7. 123 Numerical Algorithms A subsidiary alternative formulation of this scheme, not developed here, is to generate the set Ci=sign2(BBi), where alternatively the function sign2() converts positive and negative values to 1 and 0 respectively, and is the inner product. Then we define the alternative set G(A2)=n i=1Ci, which also has n×Nnscalar components, similarly corresponding to the sign changes in any direction of each scalar function fjat each of the Nnpoints. By defining (A2)(A2)as the summation of the n scalar components of G(A2)at each point, the closest points to the solutions are such that (A2)=0. 2.2 Method B Another approach, designated as Method B, is a summation-based method. This method begins with the application of the function sign to all n×Nnscalar values in matrix B, and the latter result is divided by 2n(2nis the number of vertices in each hypercell). Subsequently, the resulting values of each i-component at the 2n vertices of each cell are summed and the absolute value of the result is taken. Since the number of hyper-cells is (N−1)n, the aforementioned operation yields a set of n×(N−1)nvalues that may be referred to as G(B). Similarly to Method A,ifthe i-component of G(B)is equal to 0, this indicates that the corresponding cell contains solutions of the equation fi(x)=0. In other words, the values 0 of each i−component of G(B)identify the cells which lay at the hyper-surface defined by fi(x)=0. Similarly, those cells where all i−components of G(B)are 0 are those containing solutions of F(x)=0. The way to perform computationally the last identification is by the summation of the nscalar values of G(B)at each cell of A. This operation defines a final scalar function (B)for each cell of A. Consequently, the points at the center of cells such that (B)=0 can be used as the seeds for the final Newton-Raphson step. It can now be seen why Method Awas described as a multiplication based method, and Method Bas a summation based method as each make extensive use of these two operators, respectively. It is noteworthy that both methods Aand Bcan be reduced to a single or very few code lines in the most efficient languages, such as Python. 3 Results The GPU’s vector architecture simplifies our initial root location method to a single step following either Method Aor B, rendering comparisons with other multi-step methods moot. Instead, we will exclusively evaluate its performance using one of the most commonly used languages, Python. The implementation of the scheme in the GPU is done by using two different libraries for each method, where Cupy was utilized in Method Aand Numba in Method B. Both of these libraries make use of CUDA (Compute Unified Device Architecture). Consequently, it is necessary to install this software and libraries in the device where the algorithm is to be used (if no changes are to be made). Furthermore, the scheme was also developed in Mathematica and 123 Numerical Algorithms MATLAB, but was found to be around a factor of 100 and 10 slower, respectively. Therefore, these much slower implementations are not included here. 3.1 Performance comparison For this first subsection, a set of four nonlinear systems of equations where selected to test the effectiveness of the algorithm. All the experiments were preformed using a AMD Ryzen 7 5800HS with Radeon Graphics 3.20 GHz and 16.0 GB of installed RAM. For each experiment, the domain of each variable was either linearly discretised into up to 500 points or the domain size was multiplied by 5. Furthermore, the experiments were repeated 50 times, with the average recorded. The median was utilised to minimise extreme values. Experiment 1: Effati-Nazemi: The nonlinear system of equations has been taken from the papers of Effati [4] where f1(x1,x2)=cos(2x1)−cos(2x2)−0.4=0 f2(x1,x2)=2(x2−x1)+sin(2x2)−sin(2x1)−1.2=0(3) The domain used was −2≤xi≤2, −10 ≤xi≤10, −100 ≤xi≤100 and the number of roots inside this domain was unspecified by the literature (Table 1). Experiment 2: Solution of geometry size of thin wall rectangle girder section From Luo et al. [10], this is a non-smooth system of equations which gives a singular physical solution to the geometry size of a thin wall rectangle. If we take the non-physical approach as seen in Sidarto [12] there is a section of the domain which is undefined in the real set. Therefore, as seen here, the algorithm would fail as it only works for real functions (Table 2). f1(x1,x2,x3)=x1x2−(x1−2x3)(x2−2x3)−165 =0 f2(x1,x2,x3)=x1x3 2 12 −(x1−2x3)(x2−2x3)3 12 −9369 =0 f3(x1,x2,x3)=2(x2−x3)2(x1−x3)2x3 x2+x1−2x3 −6835 =0 (4) where D={(x1,x2,x3):−40 ≤xi≤40 i=1,2,3}. Their are 6 solutions according to the literature. The three variable system of equations can be simplified Table 1 Results for Experiment 1 Domain N◦solutions Time Method A(s) Time Method B(s) −2≤xi≤2 1 0.0030 0.0240 −10 ≤xi≤10 13 0.0125 0.0276 −100 ≤xi≤100 127 0.1000 0.1190 123 Numerical Algorithms into a two variable system of equations rearranging and eliminating x2[12] which is smooth and contained in the real set. Therefore by utilizing this form of the system we are able to arrive at the solutions (Table 2). g1(x1,x3)=x1x3 2 12 −(x1−2x3)(x2−2x3)3 12 −9369 =0 g2(x1,x3)=2(x2−x3)2(x1−x3)2x3 x2+x1−2x3 −6835 =0 x2=2x3−x1+165 2x3 (5) Experiment 3: The reactor problem The system of equations below is found in Tsoulos et al. [13] “models two continuous nonadiabatic stirred tank reactors", as described by Floudas [5]: f1(x1,x2)=(1−R)( D 10(1+β1)−x1)e 10x1 1+10x1 γ−x1 f2(x1,x2)=x1−(1+β2)x2+(1−R)( D 10 −β1x1−(1+β2)x2)e 10x2 1+10x2 γ (6) where the parameters γ=1000, D=22, βi=2, and R varies from 0.935 to 0.995 with increments of 0.005. 0 ≤xi≤1 and the number of solutions depend on the value of R varying from 1 to 7 (Table 3). Experiment 4: Chen et al. Taking the model of Chen et al. [3] found in Sidarto and Kania [12], the system of equations seen below, f1(x1,x2)=ex1−x2−sin(x1+x2) f2(x1,x2)=x2 1x2 2−cos(x1+x2)(7) has a domain of D=−10 ≤xi≤10, where the number of solutions is unspecified by the literature. Figure 3shows the contour of F(x1,x2)=0 where the black circles are the solutions found by the algorithm. We can quickly draw some conclusions from the results of these experiments. At first glance we can see Method Ais around a factor of 10 times faster than Method B (Table 4). We can also see that both methods are extremely consistent as are affected very little by the form of the function being computed. In this section a very fine mesh was utilized to ensure that all the solutions were found, but in reality such a fine mesh Table 2 Results for Experiment 2 Domain N◦solutions Time Method A(S) Time Method B(s) −40 ≤xi≤40 6 0.01600 0.1208 123 Numerical Algorithms Table 3 Results form Experiment 3 RN ◦solutions Time Method A(s) Time Method B(s) 0.935 1 0.0030 0.0260 0.940 1 0.0030 0.0260 0.945 3 0.0040 0.0270 0.950 5 0.0040 0.0260 0.955 5 0.0040 0.0260 0.960 7 0.0050 0.0270 0.965 5 0.0040 0.0259 0.970 5 0.0040 0.0260 0.975 5 0.0040 0.0260 0.980 5 0.0040 0.0260 0.985 5 0.0050 0.0260 0.990 1 0.0030 0.0260 0.995 1 0.0030 0.0270 Table 4 Results for Experiment 4 Domain N◦solutions Time Method A(s) Time Method B(s) −10 ≤xi≤10 6 0.0060 0.0210 Fig. 3 Contour plot of fi(x1,x2)=0forExperiment 4 where the black circles are the solutions arrived at utilizing the algorithm 123