## Introduction to Numerical Analysis: Nonlinear Systems of Equations

We limit our discussion here to nonlinear systems of equations composed of equations in unknowns. There are some analytical methods to solve these types of equations; however, in general, a numerical algorithm is the fastest way to find a solution. We will present two algorithms, the fixed-point iteration method and the Newton-Raphson method to solve such a system of equations.

### Fixed-Point Iteration 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.

where is a nonlinear function of the components . 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 after iteration are , and if after iteration the components are: , then, the stopping criterion would be:

where

Note that any other norm function can work as well.

#### Example

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

##### Solution

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

a = Solve[{x1^2 + x1*x2 == 10, x2 + 3 x1*x2^2 == 57}, {x1, x2}, Reals]
N[a]


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

The following is a possible rearrangement:

Using an initial guess of and yields the following:

For the next iteration, we get:

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

Using the same initial guesses, the first iteration produces:

The value of after the first iteration is:

The following Microsoft Excel table shows that convergence to and satisfying the required criterion 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


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 outputs the vector 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]


### Newton-Raphson Method

The Newton-Raphson method is the method of choice for solving nonlinear systems of equations. Many engineering software packages (especially finite element analysis software) that solve nonlinear systems of equations use the Newton-Raphson method. The derivation of the method for nonlinear systems is very similar to the one-dimensional version in the root finding section. Assume a nonlinear system of equations of the form:

If the components of one iteration are known as: , then, the Taylor expansion of the first equation around these components is given by:

Applying the Taylor expansion in the same manner for , we obtained the following system of linear equations with the unknowns being the components of the vector :

By setting the left hand side to zero (which is the desired value for the functions , then, the system can be written as:

Setting , the above equation can be written in matrix form as follows:

where is an matrix, is a vector of components and is an -dimensional vector with the components . If is invertible, then, the above system can be solved as follows:

#### Example

Use the Newton-Raphson method with to find the solution to the following nonlinear system of equations:

##### Solution

In addition to requiring an initial guess, the Newton-Raphson method requires evaluating the derivatives of the functions and . If , then it has the following form:

Assuming an initial guess of and , then the vector and the matrix have components:

The components of the vector can be computed as follows:

Therefore, the new estimates for and are:

The approximate relative error is given by:

For the second iteration the vector and the matrix have components:

The components of the vector can be computed as follows:

Therefore, the new estimates for and are:

The approximate relative error is given by:

For the third iteration the vector and the matrix have components:

The components of the vector can be computed as follows:

Therefore, the new estimates for and are:

The approximate relative error is given by:

Finally, for the fourth iteration the vector and the matrix have components:

The components of the vector can be computed as follows:

Therefore, the new estimates for and are:

The approximate relative error is given by:

Therefore, convergence is achieved after 4 iterations which is much faster than the 9 iterations in the fixed-point iteration method. The following is the Microsoft Excel table showing the values generated in every iteration:

The Newton-Raphson method can be implemented directly in Mathematica using the “FindRoot” built-in function:
View Mathematica Code

FindRoot[{x1^2 + x1*x2 - 10 == 0, x2 + 3*x1*x2^2 == 57}, {{x1, 1.5}, {x2, 3.5}}]


The following “While” loop in Mathematica provides an iterative procedure that implements the Newton-Raphson Method in Mathematica and produces the table shown:
View Mathematica Code

x = {x1, x2};
f1 = x1^2 + x1*x2 - 10;
f2 = x2 + 3 x1*x2^2 - 57;
f = {f1, f2};
K = Table[D[f[[i]], x[[j]]], {i, 1, 2}, {j, 1, 2}];
K // MatrixForm
xtable = {{1.5, 3.5}};
xrule = {{x1 -> xtable[[1, 1]], x2 -> xtable[[1, 2]]}};
Ertable = {1}
NMax = 100;
i = 1;
eps = 0.000001

While[And[Ertable[[i]] > eps, i <= NMax],
Delta = -(Inverse[K].f) /. xrule[[i]];
xtable = Append[xtable, Delta + xtable[[Length[xtable]]]];
xrule = Append[
xrule, {x1 -> xtable[[Length[xtable], 1]],
x2 -> xtable[[Length[xtable], 2]]}];
er = Norm[Delta]/Norm[xtable[[i+1]]]; Ertable = Append[Ertable, er]; i++]
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]], Ertable[[i + 1]]}, {i, 1, Length[xtable] - 1}];
T2 = Prepend[T2, Title];
T2 // MatrixForm


### Linear vs. Nonlinear Analysis of Structures: An Application

Most structural engineering problems in which the displacements and forces inside a structure are sought can be formulated using linear systems of equations. This is done by invoking certain assumptions that simplify the problem. However, when a structure undergoes large deformations, the problem has to be formulated using nonlinear systems of equations. The following is an example of a truss structure that aims at clarifying the difference.

#### Example

Two steel truss members with areas and are connected to a steel cable with an area as shown in the figure. All the connections are hinged connections allowing free relative rotations. Points , , and are fixed by means of hinges preventing displacements but allowing free rotations. Assuming that the horizontal displacement of point is to the right and the vertical displacement of point is downwards. Assume that Young’s modulus , , and . The force in each member is equal to where where is the member number, is the undeformed length, is the deformed length, is the longitudinal stress, is the strain, and is Young’s modulus in member . Find the values of and at equilibrium using: 1) Linear analysis where the displacements are assumed to be small, and 2) Nonlinear analysis. Solve twice, once with , and another time with . Do the two solution strategies give the same answer?

##### Solution

The units adopted will be for forces, for lengths, and for . The following is the displaced shape of the structure. The middle hinge move horizontally to the right a distance and vertically downards a distance . There are two unknowns in the problem and . These can be found using two equations of equilibrium at the middle hinge. The two equations are the sum of the vertical and horizontal forces, each equal to zero.

##### Linear Analysis

In linear analysis, the truss member length is assumed to only change due to displacements along its longitudinal axis. In other words, the forces in the truss member do not change if it simply rotates. In addition, the forces are assumed to be aligned with the original geometry of the structure before deformations. The strains in each member are given by: , , and . Therefore, the forces in each member can be calculated as follows:

where , , and . Assuming small deformations, equilibrium at point can be analyzed using the following forces diagram:

The first equation is the sum of horizontal forces equals to zero:

The second equation is the sum of the vertical forces equal to zero:

These equations are linear and can be written in matrix form as follows:

Setting yields a solution and , while setting yields and . Notice that in a linear analysis, the displacements are directly proportional to the applied load, that is, the ratio between at and at is equal to . The same applies to .

The following is the Mathematica code to formulate and solve the above system of equations using the “Solve” function.
View Mathematica Code

Clear[F, Ee];
Ee = 200 (10^9);
eps1 = -u1/3;
A1 = A2 = 200/1000000;
A3 = 50/1000000;
eps2 = u1/3;
Cth = 3/5;
Sth = 4/5;
eps3 = (-u1*Cth + u2*Sth)/5;
F1 = Ee*A1*eps1;
F2 = Ee*A2*eps2;
F3 = Ee*A3*eps3;
Eq1 = Expand[(F1 + F3*Cth - F2)]
Eq2 = Expand[F - F3*Sth]
sol1 = Solve[{Eq1 == 0.0, Eq2 == 0.0 /. F -> 10000}]
sol2 = Solve[{Eq1 == 0.0, Eq2 == 0.0 /. F -> 500000}]

##### Nonlinear Analysis

In a nonlinear analysis, the analysis takes into consideration the displaced position of the structure without any simplified assumptions. Looking at the first truss member, the displaced structure is shown in the following figure.

For this member, the original length , the final length , the strain and the angle of inclination is such that and .

For the second truss member, the displaced structure is shown in the following figure.

For this member, the original length , the final length , the strain and the angle of inclination is such that and .

For the third truss member, the displaced structure is shown in the following figure.

For this member, the original length , the final length , the strain and the angle of inclination is such that and .

Finally, the equations of equilibrium can be written taking into consideration the inclination of each member:

The first equation is the sum of horizontal forces equals to zero:

The second equation is the sum of the vertical forces equal to zero:

The following is a screenshot of the Mathematica output of the two equations:

Setting yields a solution and , while setting yields and . For small deformations under small applied loads, the linear analysis and nonlinear analysis produce exactly the same results. In this situation, the rotations of the members are very small such that the original shape and the deformed shape can be assumed to be the same. However, under large loads, the linear analysis is no longer valid. The two solutions at are far from each other with the nonlinear analysis solution giving less displacement. In this particular situation, as the structure deforms, the nonlinear effects cause the structure to be more rigid and therefore exhibit less deformations. Cable structures with large loads and large deformations in the cables need to be analyzed using a nonlinear analysis rather than a linear analysis.

The following is the Mathematica code to formulate and solve the above nonlinear system of equations using the “FindRoot” function with an initial guess of , and .

View Mathematica Code
Ee = 200 (10^9);
A1 = A2 = 200/1000000;
A3 = 50/1000000;
L1 = Sqrt[(3 - u1)^2 + u2^2];
Sth1 = u2/L1;
Cth1 = (3 - u1)/L1;
L2 = Sqrt[(3 + u1)^2 + u2^2];
Sth2 = u2/L2;
Cth2 = (3 + u1)/L2;
L3 = Sqrt[(3 - u1)^2 + (4 + u2)^2];
Cth3 = (3 - u1)/L3;
Sth3 = (4 + u2)/L3;
eps1 = (L1 - 3)/3;
eps2 = (L2 - 3)/3;
eps3 = (L3 - 5)/5;
F1 = Ee*A1*eps1;
F2 = Ee*A2*eps2;
F3 = Ee*A3*eps3;
Eq1 = Simplify[F1*Cth1 - F2*Cth2 + F3*Cth3]
Eq2 = Simplify[F - F1*Sth1 - F2*Sth2 - F3*Sth3]

sol1 = FindRoot[{Eq1 == 0, Eq2 == 0 /. F -> 10000}, {{u1, 0}, {u2, 0}}]
sol2 = FindRoot[{Eq1 == 0, Eq2 == 0 /. F -> 500000}, {{u1, 0}, {u2, 0}}]



### Problems

1. Use your implementation of the Newton-Raphson method to find 4 different sets of roots for the following nonlinear system of equations:

Use . Compare with the “FindRoot” built-in function in Mathematica.

2. Use your implementation of the Newton-Raphson method to find 1 set of roots for the following nonlinear system of equations:

Use . Compare with the “FindRoot” built-in function in Mathematica.

3. Compare the convergence of the Newton-Raphson method and the fixed-point iteration method with initial guesses of and to find a solution to the following nonlinear system of equations:

4. Use the Newton-Raphson method to find a solution to the following nonlinear set of equations:

Use

5. The following structural system is composed of two truss elements with unit cross sectional area. The left and right supports are hinged supports. The middle connection where the load is applied is also a hinged connection. Find the displacement at equilibrium of the following structural system. Assume a linear material with where is the stress in the bar, is the corresponding strain, and is the initial length. Use the Newton-Raphson method.

6. Solve the following nonlinear system of equations using the Newton-Raphson method and

Hint, if :