Open Educational Resources

Nonlinear Systems of Equations: Fixed-Point Iteration Method

The Method

Similar to the fixed-point iteration method for finding roots of a single equation, the fixed-point iteration method can be extended to nonlinear systems. This is in fact a simple extension to the iterative methods used for solving systems of linear equations. The fixed-point iteration method proceeds by rearranging the nonlinear system such that the equations have the form.

    \[\begin{split} x_1&=f_1(x_1,x_2,\cdots,x_n)\\ x_2&=f_2(x_1,x_2,\cdots,x_n)\\ &\vdots\\ x_n&=f_n(x_1,x_2,\cdots,x_n) \end{split} \]

where \forall i\leq n:f_i is a nonlinear function of the components x_1, x_2,\cdots, x_n. By assuming an initial guess, the new estimates can be obtained in a manner similar to either the Jacobi method or the Gauss-Seidel method described previously for linear systems of equations. Similar to linear systems of equations, the Euclidean norm can be used to check convergence. So, if the components of the vector x after iteration i are x_1^{(i)}, x_2^{(i)},\cdots,x_n^{(i)}, and if after iteration i+1 the components are: x_1^{(i+1)}, x_2^{(i+1)},\cdots,x_n^{(i+1)}, then, the stopping criterion would be:

    \[ \frac{\|x^{(i+1)}-x^{(i)}\|}{\|x^{(i+1)}\|}\leq \varepsilon_s \]


    \[\begin{split} \|x^{(i+1)}-x^{(i)}\|&=\sqrt{(x_1^{(i+1)}-x_1^{(i)})^2+(x_2^{(i+1)}-x_2^{(i)})^2+\cdots+(x_n^{(i+1)}-x_n^{(i)})^2}\\ \|x^{(i+1)}\|&=\sqrt{(x_1^{(i+1)})^2+(x_2^{(i+1)})^2+\cdots+(x_n^{(i+1)})^2}\\ \end{split} \]

Note that any other norm function can work as well.


Use the fixed-point iteration method with \varepsilon_s=0.0001 to find the solution to the following nonlinear system of equations:

    \[ x_1^2+x_1x_2=10\qquad x_2+3x_1x_2^2=57 \]


The exact solution in the field of real numbers for this system can actually be obtained using Mathematica as shown in the code below.

View Mathematica Code
Clear[x1, x2];
a = Solve[{x1^2 + x1*x2 == 10, x2 + 3 x1*x2^2 == 57}, {x1, x2}, Reals]
View Python Code
import sympy as sp
x1, x2 = sp.symbols('x1 x2')
a = sp.nonlinsolve([x1**2 + x1*x2 - 10, x2 + 3*x1*x2**2 - 57], [x1, x2])

The fixed-point iteration numerical method requires rearranging the equations first to the form:

    \[ x=\left(\begin{array}{c}x_1\\x_2\end{array}\right)=f(x)=\left(\begin{array}{c}f_1(x_1,x_2)\\f_2(x_1,x_2)\end{array}\right) \]

The following is a possible rearrangement:

    \[ \left(\begin{array}{c}x_1\\x_2\end{array}\right)=\left(\begin{array}{c}\frac{10-x_1^2}{x_2}\\57-3x_1x_2^2\end{array}\right) \]

Using an initial guess of x_1=1.5 and x_2=3.5 yields the following:

    \[ x_1=\frac{10-(1.5)^2}{3.5}=2.21429\Rightarrow x_2=57-3(2.21429)(3.5)^2=-24.37516 \]

For the next iteration, we get:

    \[ x_1=\frac{10-(2.21429)^2}{-24.37516}=-0.20910\Rightarrow x_2=57-3(-0.20910)(-24.37516)^2=429.709 \]

Continuing the procedure shows that it is diverging. A different rearrangement for the equations has the form:

    \[ \left(\begin{array}{c}x_1\\x_2\end{array}\right)=\left(\begin{array}{c}\sqrt{10-x_1x_2}\\\sqrt{\frac{57-x_2}{3x_1}}\end{array}\right) \]

Using the same initial guesses, the first iteration produces:

    \[ x_1=\sqrt{10-1.5\times 3.5}=2.1794 \Rightarrow x_2=\sqrt{\frac{57-3.5}{3(2.1794)}}=2.86051 \]

The value of \varepsilon_r after the first iteration is:

    \[ \varepsilon_r=\frac{\sqrt{(1.5-2.1794)^2+(2.8605-3.5)^2}}{\sqrt{(2.1794)^2+(2.86051)^2}}=0.2595 \]

The following Microsoft Excel table shows that convergence to x_1=2.0 and x_2=3.0 satisfying the required criterion \varepsilon_r<\varepsilon_s=0.0001 is achieved after 9 iterations.


The following code does the same thing in Mathematica to produce the table below:

View Mathematica Code
xtable = {{1.5, 3.5}};
ErrorTable = {1};
Nmax = 100;
eps = 0.0001;
i = 2;
While[And[ErrorTable[[i - 1]] > eps, i <= Nmax], 
  x1new = Sqrt[10 - xtable[[i - 1, 1]]*xtable[[i - 1, 2]]]; 
  x2new = Sqrt[(57 - xtable[[i - 1, 2]])/3/x1new]; 
  xnew = {x1new, x2new};
  xtable = Append[xtable, xnew]; 
  er = Norm[xnew - xtable[[i - 1]]]/Norm[xnew]; 
  ErrorTable = Append[ErrorTable, er]; i = i + 1];
Title = {"Iteration", "x1_in", "x2_in", "x1_out", "x2_out", "Er"};
T2 = Table[{i, xtable[[i, 1]], xtable[[i, 2]], xtable[[i + 1, 1]], xtable[[i + 1, 2]], ErrorTable[[i + 1]]}, {i, 1, Length[xtable] - 1}];
T2 = Prepend[T2, Title];
T2 // MatrixForm
View Python Code
import numpy as np
import pandas as pd
xtable = np.array([[1.5, 3.5]])
ErrorTable = [1]
Nmax = 100
eps = 0.0001
i = 1
while ErrorTable[i - 1] > eps and i <= Nmax: 
  x1new = np.sqrt(10 - xtable[i - 1, 0]*xtable[i - 1, 1]) 
  x2new = np.sqrt((57 - xtable[i - 1, 1])/3/x1new)
  xnew = [x1new, x2new]
  xtable = np.concatenate((xtable,[xnew])) 
  er = np.linalg.norm(xnew - xtable[i - 1])/np.linalg.norm(xnew) 
  i += 1
T2 = [[i+1, xtable[i][0], xtable[i][1], xtable[i + 1][0], xtable[i + 1][1], ErrorTable[i + 1]] for i in range(len(xtable)-1)]
pd.DataFrame(T2, columns=["Iteration", "x1_in", "x2_in", "x1_out", "x2_out", "Er"])

MATLAB files for the fixed-point iteration example:
Download MATLAB file 1 (fpisystem.m)
Download MATLAB file 2 (g1.m)
Download MATLAB file 3 (g2.m)

The example here shows that the fixed-point iteration method is not guaranteed to give a possible solution. In fact, the initial guess and the form chosen affect whether a solution can be obtained or not. Note that the “FixedPointList” built-in function in Mathematica can be used to implement the method with an initial guess. Notice in the code below how the function g outputs the vector x as a list and that the second component uses the output of the first component:

View Mathematica Code
g[x_] := (x1 = Sqrt[10 - x[[1]]*x[[2]]]; {x1, Sqrt[(57 - x[[2]])/3/x1]})
FixedPointList[g[#] &, {1.5, 3.5}, 20]
View Python Code
import numpy as np
def g(x):
  x1 = np.sqrt(10 - x[0]*x[1])
  return [x1, np.sqrt((57 - x[1])/3/x1)]
x = [1.5, 3.5]
for i in range(20):
  x = g(x)

Lecture Video

Leave a Reply

Your email address will not be published. Required fields are marked *