Nonlinear Algebraic Equations Solved by an Optimal Splitting-Linearizing Iterative Method
1 Center of Excellence for Ocean Engineering, Center of Excellence for the Oceans, National Taiwan Ocean University, Keelung, 202301, Taiwan
2 Department of Mathematics, College of Sciences and Humanities in Al-Kharj, Prince Sattam bin Abdulaziz University, Alkharj, 11942, Saudi Arabia
3 Department of Basic Engineering Science, Faculty of Engineering, Menofia University, Shebin El-Kom, 32511, Egypt
4 Department of Marine Engineering, National Taiwan Ocean University, Keelung, 202301, Taiwan
* Corresponding Author: Yung-Wei Chen. Email:
Computer Modeling in Engineering & Sciences 2023, 135(2), 1111-1130. https://doi.org/10.32604/cmes.2022.021655
Received 26 January 2022; Accepted 31 May 2022; Issue published 27 October 2022
AbstractHow to accelerate the convergence speed and avoid computing the inversion of a Jacobian matrix is important in the solution of nonlinear algebraic equations (NAEs). This paper develops an approach with a splitting-linearizing technique based on the nonlinear term to reduce the effect of the nonlinear terms. We decompose the nonlinear terms in the NAEs through a splitting parameter and then linearize the NAEs around the values at the previous step to a linear system. Through the maximal orthogonal projection concept, to minimize a merit function within a selected interval of splitting parameters, the optimal parameters can be quickly determined. In each step, a linear system is solved by the Gaussian elimination method, and the whole iteration procedure is convergent very fast. Several numerical tests show the high performance of the optimal split-linearization iterative method (OSLIM).
Nonlinear partial differential equations (PDEs) frequently appear in many critical natural sciences and engineering technology problems. Mathematical models in physics, mechanics, and other disciplines fields can be applied to describe physical phenomena. Although some phenomena may be modeled by the linear type PDEs, the assumption of more complex conditions on the material property and geometric domain of the problem requires the use of nonlinear type PDEs. In most situations, solving nonlinear PDEs analytically and exactly is impossible. Because of the crucial role of the nonlinear problems, the researchers provided different numerical methods to obtain the solutions of nonlinear PDEs, which are properly transformed into the nonlinear algebraic equations (NAEs) by using numerically discretized methods. As mentioned in [1–4], nonlinear optimization problem, nonlinear programming, nonlinear sloshing problem, and nonlinear complementarity problem (NCP) can be recast to the issues to solve NAEs, with the help of NCP-functions, for example, the Fischer-Burmeister NCP . Indeed, the numerical solution of NAEs is one of the main problems in computational mathematics. Among the many numerical methods, Newton’s method and its improvements have now been widely developed; however, if the initial guess of the solution is incorrect, these algorithms may fail. In general, it is difficult for most NAEs to choose a good initial condition. Therefore, it is necessary to develop a more effective algorithm that is less sensitive to the initial guess of the solution and converges quickly.
For the vector form of the NAEs:
the Newton method is given by
where is an Jacobian matrix with its ij-th component given by . Starting from an initial guess , Eq. (2) with the computed at each iteration can be used to generate a sequence of , . When is convergent under a specified convergence criterion, the solution of Eq. (1) is obtained.
The Newton method possesses a significant advantage in that it is quadratically convergent. However, the main drawbacks are sensitive to the initial point and need to compute at each iteration step. Therefore, some quasi-Newton methods have been developed to overcome these shortcomings of the Newton method. For example, in order to eliminate the requirement of finding and inverting the Jacobian matrix in Newton’s iterative algorithm, Liu et al.  proposed a derivative-free first-order system of nonlinear ordinary differential equations (ODE):
where ν is a nonzero constant and q(t) may generally be a monotonically increasing function of t. The term plays a vital role in stabilizing the controller in their approach. Even if the initial guess of the solution is not proper, it can help obtain a stable solution and accelerate the convergent speed. Liu et al. [2,6,7] showed that the efficiency and performance of the fictitious time integration method (FTIM) combined with the group preserving scheme  to integrate the resultant ODEs in Eq. (3). Atluri et al.  proposed an improvement of Newton’s method, which combined two methods of Newton’s method and FTIM, without the inversion of . Related to the FTIM, some papers are employed to solve nonlinear problems [10,11]. To give a motivation for the present study, we begin with a system of linear algebraic equations (LAEs):
where A and are respectively a non-singular matrix and an input vector.
Let and suppose that D is non-singular. Setting
where denotes the kth step unknown variables . The splitting of A to used in Eq. (6) was proposed by Jacobi to solve the linear algebraic equations (LAEs) iteratively. In general, A can be uniquely decomposed into
where U and L are a strict upper- and lower triangular matrix of A, respectively.
where w is a nonzero relaxation factor. When , it is known as the Gauss-Seidel iterative method. If , , from Eq. (8), can be easily found by the forward substitution method. The main reason is a lower triangular matrix.
Nowadays, there have different splitting techniques to decompose A as
where M is a non-singular matrix. Inserting Eq. (9) into Eq. (4), moving to the right-hand side, taking the value of to be and on the left-hand side taking the value of to be . Then, the resulting iterative scheme expresses as follows:
is called an iteration matrix. The sufficient and necessary condition of the convergence for Eq. (10) is , where ρ is the spectral radius, defined as the absolute value of the largest eigenvalue in the magnitude of . When gets smaller and smaller, the convergence speed is accelerated.
For the LAE, there are many discussions about the split method to determine the best relaxation factor w in SOR [16–20]. In contrast, few results discuss the relationship between the split method and NAEs. In this paper, the split method based on the nonlinear term is applied to reduce the influence of nonlinear terms and further avoid selecting the parameters by a trail and error method. Hence, this article will propose a novel splitting-linearizing technique to iteratively solve the NAE and develop a new algorithm for determining the optimal splitting parameters.
where A and B are both matrices and are constant n-vector. While A is a constant matrix, B is a matrix function of x. When Ax incorporates all linear terms in the NAEs, collects the nonlinear terms.
Then, we can decompose Eq. (11) to
The idea of the splitting technique in [21,22] has been proposed to solve a single nonlinear equation. Future, the nonlinear term in Eq. (11) is first written as . Then, while the term is kept on the left-hand side, the other is moved to the right-hand side, as shown in Eq. (12). Above decomposition is a new splitting application, using a splitting parameter to be determined below.
We give an initial guess to initialize the iteration. At the kth step of the iteration, is known, and Eq. (12) is linearized around to
For , it is a linear equations system, since the coefficients A and are constant matrices. The right-hand side of Eq. (13) is a constant vector because is available from the previous step. A new linear system can be derived and applied to determine the new at the next step by
The above process is a novel splitting-linearizing method for solving the NAEs.
In this study, the Gaussian elimination method solves Eq. (14) to determine at each iterative step. However, before solving Eq. (14), an unknown parameter w is given in Eq. (15). Instead of being given a constant value w in the iterative procedure of the proposed method, which is determined through some trial and error, we can formulate a new method to determine the optimal value of at each step.
For a linear system
to obtain a fast descent direction c for a given residual vector b used in the iterative solution of Eq. (16).
The minimization in Eq. (17) is equivalent to maximize
which is the orthogonal projection of onto y.
and insert them into Eq. (17), we can derive
We can search the minimal value of by
in a given interval by taking , , for example, , and . can be used to adjust the search direction. Inserting into Eqs. (24)–(26), we can compute in Eq. (23), and then we take the optimal value of to obtain the minimization of . The optimal value of in can quickly be obtained. This new algorithm with the optimal value is called the optimal split-linearization iterative method (OSLIM).
In the OSLIM, there is an interesting feature that by using the Cauchy-Schwarz inequality in Eq. (20),
it readily leads to
The procedures of the OSLIM to solve the NAEs are summarized as follows:
(i) Give n and nonlinear equations, , and .
(ii) Give an initial guess .
(iii) For we repeat the follwing compuations:
from Eq. (27).
If converges , this iterative procedure stops; otherwise, go to step (iii) for the next iteration. Here, ε is a given convergence criterion and named the residual norm. Also, the residual norm can be checked by
To further assess the performance of OSLIM, we investigate it using the convergence order. For solving a scalar equation , the numerically computed order of convergence (COC) is approximated by 
where r is a solution of and the sequence is generated from an iterative scheme. In the computation of COC, we store the values of and take , where is the number of iterations for the convergence.
For solving the NAEs, we can generate the sequence , and store them in the program. We define
where is the number of iterations for satisfying the convergence criterion, is exact solution, and is the error in terms of Euclidean norm at the kth step. If the values of are unknown, we take .
We first consider :
Two roots are (1, 1) and (1, −1) obviously; however, there exist other roots.
As noted by Yeih et al. , the Newton method fails because the Jacobi matrix
is singular at . Although the scalar homotopy method  considered the optimal hybrid search directions, the order of the residual norm only satisfies 10−4. When solving this one by the residual-norm based algorithm (RNBA) , the residual norm satisfies the order of 10−7. Above iteration algorithms based on the evolution of the residual vector are not satisfied the given convergence criterion, and the third solution cannot be obtained by these algorithms.
Using the OSLIM to solve this problem, a translation is applied as follows:
and choose sufficiently large to render such that we have
Eq. (37) can be written as
By using the OSLIM, we take , , and , and starting from (1.5, 1.5) and through 34 iterations it converges under as shown in Fig. 1a. The third solution is obtained as = (−0.4776700623, 1.331101541), and the residuals are = (4.44 × 10−16, 6.66 × 10−16). In Fig. 1b, the gradient direction is quickly found at the 5th step, and then the value of quickly tend to the minimum value . The algorithm adjusts the gradient direction to approach the solution to obtain a smaller step size. As shown in Fig. 1c, the values of w change drastically after 15 steps, and the change value of w is between . For this example, COC = 1.2838 is obtained, which indicates that the OSLIM converges faster than the iterative scheme with linear convergence speed.
Considering two-variable nonlinear equations  as follows:
Hirsch et al.  used the continuous Newton algorithm to calculate this problem; however, the numerical results are inaccurate. Under the same situation, the FTIM , modified Newton method , and RNBA  have been applied and obtained three solutions from three different initial guessing values. In Liu et al. , this problem solved by the optimal vector drive algorithm (OVDA) can find five solutions. As stated in the literature , the attracting sets for several fixed points of problem are undefined. Therefore, various methods may give different solutions, even with the same initial value.
or written as
Eqs. (42) and (43) are of the form in Eq. (14), but with different B. For the uniqueness of B, for example, we can put the product term as , such that in the brace is viewed as the coefficient which is inserted into the coefficient matrix B, while x is put into the unknown vector. The principle is that the former variable x is preferred over the latter variable y placed in the unknown vector.
To compare with the result obtained by Atluri et al. , we take , , , , and . In , starting from = (0.1, 0.1) and through 101 steps the solution = (0.134212, 0.811128) is obtained with the residuals being = (1.141 × 10−7, 8.913 × 10−7). By using the OSLIM, we take , and , and it converges with 35 iterations under as shown in Fig. 2a. The solution = (−0.16363472339, 0.23052874358) is obtained, and the residuals = (2.22 × 10−15, 5.33 × 10−15) are much smaller than that in ; moreover, the OSLIM converges faster than that obtained by Atluri et al. . In Fig. 2b, the gradient direction is quickly found at the 6th step, and then the value of quickly tend to the minimum value . The values of w change drastically after 19 steps, and the change value of w is between = [−1, −0.5] as shown in Fig. 2c. For this example, COC = 1.0895 is obtained.
Considering a system of three algebraic equations with three variables:
Obviously is a solution.
By using the OSLIM, we take , and , and starting from = (0.5, 0.5, 0.6), through 33 iterations it converges under as shown in Fig. 3a. The solution (0.9305422841, 1.218366932, 0.8510907842) is obtained, and the residuals are = (0, 3.553 × 10−15, 1.332 × 10−15). The OSLIM converges faster and more accurate than obtained by Atluri et al. , where 140 steps are required, and the residual errors are (6.046 × 10−4, 7.574 × 10−4, 3.005 × 10−6). In Fig. 3b, the gradient direction is quickly found at the 6th step, and then the value of quickly tend to the minimum value . When approaching solutions, the values of w of the algorithm adjusts the iterative step size after 15 steps, and the change value of w is between = [−2, −1.5] as shown in Fig. 3c. For this example, COC = 1.0429 is obtained.
Here, considering the following boundary value problem:
The exact solution is
By using the finite difference discretisation of u at the grid points, we can obtain
where is the element length. In order to conviencely describle the coefficient matrices of B as Eq. (21), , , and can be expressed as follows:
Setting n = 39, , , ε = 10−10 and , the OSLIM converges within 5 iterations as shown in Fig. 4a. In Fig. 4b, the gradient direction is quickly found at the 2nd step when the value of quickly tend to the minimum value . When approaching solutions, the values of w of the algorithm adjusts the iterative step size after the 4th step, and the change value of w is between = [−1, −0.25] as shown in Fig. 4c. The maximum error obtained by the OSLIM is 2.98 × 10−4, which converges faster and is more accurate than Atluri et al. . In , the iteration number is more than 8000 steps under the convergence criterion ε = 10−4, and the errors are in the order of 10−3. COC = 1.9741 is obtained, which reveals that the OSLIM is almost quadratic convergence.
Then, considering a test example given by Krzyworzcka  as follows:
To convince describe the coefficient matrices of B as Eq. (21), , , and can be expressed as follows:
Setting , , ε = 10−10, and initial values = 1, , the OSLIM converges within 10 iterations as shown in Fig. 5a. In Fig. 5b, the gradient direction is quickly found at the 5th step and then the value of quickly tend to the minimum value . When approaching solutions, the values of w of the algorithm adjusts the iterative step size after the 7th steps, and the change value of w is between = [−1, 0.8] as shown in Fig. 5c. The calculation results of , , are shown in Table 1. The OSLIM converges faster than that obtained by Atluri et al. , where 55 steps are required to satisfy a convergence criterion ε = 10−6. COC = 1.9996 is obtained, which shows again that the OSLIM is quadratic convergence.
The following example is given by Roose et al. :
The OSLIM can adopt two different formulations to solve Eq. (53). The two different formulations lead to different coefficient matrices of B. First, we adopt
To convince desdescribee coefficient matrices of B as Eq. (21), can be expressed as follows:
where the number with a superscript k denotes the kth step value, and it is inserted in the coefficient matrix when we derive the linear system. We modify the term
in the second line to
such that we have a second formulation:
To convince describe the coefficient matrices of B as Eq. (21), can be expressed as follows:
Setting n = 10, , , ε = 10−10, and initial values = , , the first OSLIM with the first formulation converges within 14 iterations as shown in Fig. 6a. COC = 1.0651 is obtained. In Fig. 6b, the gradient direction is quickly found at the 4th step and then the value of quickly tend to the minimum value . In Fig. 6c, the value of w of the algorithm adjusts the iterative step size after the 8th step, and the change value of w is between = [−0.1, 0.04] as shown in Fig. 6c. The OSLIM converges faster than that obtained by Atluri et al. , where over 700 steps are required under ε = 10−6. The residual errors of and the calculation results of , , are listed in Tables 2 and 3, respectively.
Setting n = 10, , , ε = 10−10, and initial values = , , the second OSLIM with the second formulation converges within 13 iterations as shown in Fig. 7a. However, the solutions are different from those obtained from the first OSLIM with the first formulation. When n is raised to 50, both OSLIM can find solutions quickly with 12 and 15 iterations, respectively. These two solutions are compared in Fig. 7b.
Considering the following a two-dimensional nonlinear Bratu equation:
and the domain is an ellipse with semi-axes lengths a = 1 and b = 0.5. The Dirichlet boundary condition is imposed. Due to u > 0, we decompose it to
where is the value at the kth step, expanded by
After collocating points inside the domain and on the boundary, simultaneously, to satisfy the Eq. (56) and the boundary condition equation, we can derive NAEs to determine the coefficients iteratively.
By using the OSLIM, we take , and to determine w by using the Newton method, through 24 iterations it converges under ε = 10−5 as shown in Fig. 8a. COC = 1.4926 is obtained. The gradient direction is quickly found at the 7th step, and then the value of quickly tend to the minimum value as shown in Fig. 8b. In order to satisfy , we can find that the algorithm adjusts the search direction after the 3rd step and the change value of w is between = [−0.25, 1.78] as shown in Fig. 8c. The maximum error obtained by the OSLIM is 7.95 × 10−6. Interestingly, the OSLIM can be used to solve a highly nonlinear PDE with high accuracy.
When evaluating the performance of newly developed iterative solutions for solving nonlinear algebraic equations (NAE), there are two critical factors: accuracy and convergence speed. In this article, we first decomposed the nonlinear terms in the NAEs through a splitting parameter and then we linearized the NAEs around the values at the previous step and derived a system of linear algebraic equations (LAEs). Through the maximum orthogonal projection concept at each iterative step, we select the optimal value of the splitting parameter and use the Gaussian elimination method to solve the LAE. The features of the proposed method are summarized as follows: (a) The merit function is in terms of the coefficient matrix, right-side vector and the value of unknown vector at the previous step. (b) It is easy to minimize the search within a preferred interval through some operations. (c) The current OSLIM is insensitive to the initial guess and without needing the inversion of the Jacobian matrix at each iterative step, which can save much of the computational time. We have successfully applied the OSLIM to solve the NAEs resulting from the nonlinear boundary value problem in one- and two-dimensional spatial spaces. Through the numerical test, the convergence order of the OSLIM is between linear and quadratic convergence.
Funding Statement: The corresponding author appreciates the financial support provided by the Ministry of Science and Technology, Taiwan, ROC under Contract No. MOST 110-2221-E-019-044.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
- Liu, C. S., & Atluri, S. N. (2008). A fictitious time integration method (FTIM) for solving mixed complementarity problems with applications to non-linear optimization. Computer Modeling in Engineering & Sciences, 34(2), 155-178. [Google Scholar] [CrossRef]
- Kwak, B. M., & Kwak, J. Y. (2021). Binary NCP: A new approach for solving nonlinear complementarity problems. Journal of Mechanical Science and Technology, 35(3), 1161-1166. [Google Scholar] [CrossRef]
- Shih, C. F., Chen, Y. W., Chang, J. R., & Soon, S. P. (2021). The equal-norm multiple-scale trefftz method for solving the nonlinear sloshing problem with baffles. Computer Modeling in Engineering & Sciences, 127(3), 993-1012. [Google Scholar] [CrossRef]
- Nurhidayat, I., Hao, Z., Hu, C. C., & Chen, J. S. (2020). An ordinary differential equation approach for nonlinear programming and nonlinear complementarity problem. International Journal of Industrial Optimization, 1(1), 1-14. [Google Scholar] [CrossRef]
- Fischer, A. (1992). A special newton-type optimization method. Optimization, 24(3–4), 269-284. [Google Scholar] [CrossRef]
- Liu, C. S., & Atluri, S. N. (2008). A novel time integration method for solving a large system of non-linear algebraic equation. Computer Modeling in Engineering & Sciences, 31(2), 71-83.15. [Google Scholar] [CrossRef]
- Liu, C. S. (2009). A fictitious time integration method for a quasilinear elliptic boundary value problem, defined in an arbitrary plane domain. Computers, Materials & Continua, 11(1), 15-32. [Google Scholar] [CrossRef]
- Liu, C. S. (2001). Cone of non-linear dynamical system and group preserving schemes. International Journal of Non-Linear Mechanics, 36, 1047-1068. [Google Scholar] [CrossRef]
- Atluri, S. N., Liu, C. S., & Kuo, C. L. (2009). A modified newton method for solving non-linear algebraic equations. Journal of Marine Science and Technology, 17(3), 238-247. [Google Scholar] [CrossRef]
- Jones, C., & Tian, H. (2017). A time integration method of approximate fundamental solutions for nonlinear poisson-type boundary value problems. Communications in Mathematical Sciences, 15(3), 693-710. [Google Scholar] [CrossRef]
- Hashemi, M. S., Mustafa, I., Parto-Haghighi, M., & Mustafa, B. (2019). On numerical solution of the time-fractional diffusion-wave equation with the fictitious time integration method. European Physical Journal Plus, 134, 488. [Google Scholar] [CrossRef]
- Demmel, J. W. (1997). Applied numerical linear algebra. Philadelphia: SIAM.
- Saad, Y. (2003). Iterative methods for sparse linear systems. Philadelphia: SIAM.
- Quarteroni, A., Sacco, R., Saleri, F. (2000). Numerical mathematics. New York: Springer Science.
- Hadjidimos, A. (2000). Successive overrelaxation (SOR) and related methods. Journal of Computational and Applied Mathematics, 123(1–2), 177-199. [Google Scholar] [CrossRef]
- Bai, Z. Z., & Chi, X. B. (2003). Asymptotically optimal successive overrelaxation method for systems of linear equations. Journal of Computational Mathematics, 21(5), 603-612. [Google Scholar]
- Wen, R. P., Meng, G. Y., & Wang, C. L. (2013). Quasi-Chebyshev accelerated iteration methods based on optimization for linear systems. Computers & Mathematics with Applications, 66(6), 934-942. [Google Scholar] [CrossRef]
- Meng, G. Y. (2014). A practical asymptotical optimal SOR method. Applied Mathematics and Computation, 242, 707-715. [Google Scholar] [CrossRef]
- Miyatake, Y., Sogabe, T., & Zhang, S. L. (2020). Adaptive SOR methods based on the wolfe conditions. Numerical Algorithms, 84, 117-132. [Google Scholar] [CrossRef]
- Yang S., G., & K, M. (2009). The optimal relaxation parameter for the SOR method applied to the poisson equation in any space dimensions. Applied Mathematics Letters, 22(3), 325-331. [Google Scholar] [CrossRef]
- Liu, C. S. (2020). A new splitting technique for solving nonlinear equations by an iterative scheme. Journal of Mathematics Research, 12(4), 40-48. [Google Scholar] [CrossRef]
- Liu, C. S., Hong, H. K., & Lee, T. L. (2021). A splitting method to solve a single nonlinear equation with derivative-free iterative schemes. Mathematics and Computers in Simulation, 190, 837-847. [Google Scholar] [CrossRef]
- Liu, C. S. (2012). A globally optimal iterative algorithm to solve an ill-posed linear system. Computer Modeling in Engineering & Sciences, 84(4), 383-403. [Google Scholar] [CrossRef]
- Liu, C. S. (2013). An optimal tri-vector iterative algorithm for solving ill-posed linear inverse problems. Inverse Problems in Science and Engineering, 21, 650-681. [Google Scholar] [CrossRef]
- Liu, C. S. (2014). A globally optimal tri-vector method to solve an ill-posed linear system. Journal of Computational and Applied Mathematics, 260, 18-35. [Google Scholar] [CrossRef]
- Yeih, W. C., Ku, C. Y., Liu, C. S., & Chan, I. Y. (2013). A scalar homotopy method with optimal hybrid search directions for solving nonlinear algebraic equations. Computer Modeling in Engineering & Sciences, 90(4), 255-282. [Google Scholar] [CrossRef]
- Chen, Y. W. (2014). A residual-norm based algorithm for solving determinate/indeterminate systems of non-linear algebraic equations. Journal of Marine Science and Technology, 22(5), 583-597. [Google Scholar] [CrossRef]
- Hirsch, M., & Smale, S. (1979). On algorithms for solving f(x) = 0. Communications on Pure & Applied Mathematics, 32(3), 281-312. [Google Scholar] [CrossRef]
- Liu, C. S., & Atluri, S. N. (2011). An iterative algorithm for solving a system of nonlinear algebraic equations,. Computer Modeling in Engineering & Sciences, 73(4), 395-431. [Google Scholar] [CrossRef]
- Krzyworzcka, S. (1996). Extension of the lanczos and CGS methods to systems of nonlinear equations. Journal of Computational and Applied Mathematics, 69(1), 181-190. [Google Scholar] [CrossRef]
- Roose, A., Kulla, V., Lomb, M., Meressoo, T. (1990). Test examples of systems of non-linear equations. Tallin: Estonian Software and Computer Service Company.