Linear and Nonlinear Algebraic Equation Systems

Linear Systems
In root of equations section, we determined the value x that satisfied a single equation, f(x) = 0. Now, we
deal with the case of determining the values 𝑥1 , 𝑥2 , . . . 𝑥𝑛 that simultaneously satisfy a set of equations:
𝑓1 (𝑥1 , 𝑥2 , … , 𝑥𝑛 ) = 0
𝑓2 (𝑥1 , 𝑥2 , … , 𝑥𝑛 ) = 0
⋮
𝑓𝑛 (𝑥1 , 𝑥2 , … , 𝑥𝑛 ) = 0
Such systems are either linear or nonlinear.
1 Linear Algebraic Equations and Matrices
An important problem in structural engineering is that of
finding the forces and reactions associated with a statically
determinate truss. Following figure shows an example of
such a truss. The forces (F) represent either tension or
compression on the members of the truss. External
reactions (𝐻2 , 𝑉2 , 𝑎𝑛𝑑 𝑉3) are forces that characterize how
the truss interacts with the supporting surface. The hinge at
node 2 can transmit both horizontal and vertical forces to
the surface, whereas the roller at node 3 transmits only
vertical forces. It is observed that the effect of the external
loading
of 1000 lb is distributed among the various members of
the truss.
Free-body force diagrams are shown for each node in the
figure. The sum of the forces in both horizontal and
vertical directions must be zero at each node, because the
system is at rest. Therefore, for node 1,
∑𝐹𝐻 = 0 = −𝐹1 𝑐𝑜𝑠30𝑜 + 𝐹3 𝑐𝑜𝑠60𝑜 + 𝐹1,ℎ
∑𝐹𝑉 = 0 = −𝐹1 sin 30𝑜 − 𝐹3 sin 60𝑜 + 𝐹1,𝑣
For node 2,
∑𝐹𝐻 = 0 = 𝐹2 + 𝐹1 𝑐𝑜𝑠30𝑜 + 𝐹2,ℎ + 𝐻2
∑𝐹𝑉 = 0 = 𝐹1 sin 30𝑜 + 𝐹2,𝑣 + 𝑉2
For node 3,
∑𝐹𝐻 = 0 = −𝐹2 − 𝐹3 𝑐𝑜𝑠60𝑜 + 𝐹3,ℎ
∑𝐹𝑉 = 0 = 𝐹3 sin 60𝑜 + 𝐹3,𝑣 + 3
where 𝐹𝑖,ℎ is the external horizontal force applied to node i (where a positive force is from left to right) and
𝐹1,𝑦 is the external vertical force applied to node i (where a positive force is upward). Thus, in this problem,

Page 1 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
the 1000-lb downward force on node 1 corresponds 𝑡𝑜 𝐹1,𝑦 = −1000. For this case all other
𝐹𝑖,𝑦 ’𝑠 𝑎𝑛𝑑 𝐹𝑖,ℎ ’𝑠 are zero. Note that the directions of the internal forces and reactions are unknown. Proper
application of Newton’s laws requires only consistent assumptions regarding direction. Solutions are
negative if the directions are assumed incorrectly. Also note that in this problem, the forces in all members
are assumed to be in tension and act to pull adjoining nodes together.
A negative solution therefore corresponds to compression. This problem can be written as the following
system of six equations and six unknowns:
0.8660 0 −0.5 0 0 0 𝐹1 0
0.5 0 0.866 0 0 0 𝐹2 −1000
−0.8660 −1 0 −1 0 0 𝐹3 0
=
−0.5 0 0 0 −1 0 𝐻2 0
0 1 0.5 0 0 0 𝑉2 0
[ 0 0 −0.8660 0 0 −1] { 𝑉3 } { 0 }

2 Matrix Algebra Overview

A matrix consists of a rectangular array of elements represented by a
single symbol. As depicted in the following figure. The matrix has m
rows and n columns and is said to have a dimension of m by n (or m
× n). It is referred to as an m by n matrix.
Matrices where m = n are called square matrices. For example, a 3 × 3
matrix is:
𝑎11 𝑎12 𝑎13
[𝐴] = [𝑎21 𝑎22 𝑎23 ]
𝑎31 𝑎32 𝑎33
The diagonal consisting of the elements 𝑎11 , 𝑎22 , and 𝑎33 is termed the principal or main diagonal of the
matrix. A symmetric matrix is one where the rows equal the columns—that is, 𝑎𝑖𝑗 = 𝑎𝑗𝑖 for all i’s and j’s.
For example,
5 1 2
[𝐴] = [1 3 7]
2 7 8
is 3 × 3 symmetric matrix. A diagonal matrix is a square matrix where all elements off the main diagonal
are equal to zero, as in:
𝑎11 0 0
[𝐴] = [ 0 𝑎22 0 ]
0 0 𝑎33
An identity matrix is a diagonal matrix where all elements on the main diagonal are equal to 1, as in
1 0 0
[𝐼] = [0 1 0]. An upper triangular matrix is one where all the elements below the main diagonal are
0 0 1
zero, as in:

Page 2 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
𝑎11 𝑎12 𝑎13
[𝐴] = [ 0 𝑎22 𝑎23 ]
0 0 𝑎33
A lower triangular matrix is one where all elements above the main diagonal are zero, as in
𝑎11 0 0
[𝐴] = [𝑎21 𝑎22 0 ]
𝑎31 𝑎32 𝑎33
The product of two matrices is represented as [C] = [A][B], where the elements of [C] are defined as
𝑛

𝐶𝑖𝑗 = ∑ 𝑎𝑖𝑘 𝑏𝑘𝑗

𝑘=1

where n = the column dimension of [A] and the row dimension of [B]. That is, the 𝑐𝑖𝑗 element is obtained by
adding the product of individual elements from the 𝑖 𝑡ℎ row of the first matrix, in this case [A], by the jth
column of the second matrix [B]. The following figure depicts how the rows and columns line up in matrix
multiplication.

The transpose of a matrix involves transforming its rows into columns and its columns into rows. For
example, for the matrix [𝐴]3 × 3 matrix, the transpose, designated [A]T, is defined as
𝑎11 𝑎21 𝑎31
[𝐴}𝑇 == [𝑎12 𝑎22 𝑎32 ]
𝑎13 𝑎23 𝑎33
2.1 Representing Linear Algebraic Equations in Matrix Form
It should be clear that matrices provide a concise notation for representing simultaneous linear equations.
For example, a 3 × 3 set of linear equations,
𝑎11 𝑥1 + 𝑎12 𝑥2 + 𝑎13 𝑥3 = 𝑏1
𝑎21 𝑥1 + 𝑎22 𝑥2 + 𝑎23 𝑥3 = 𝑏2
𝑎31 𝑥1 + 𝑎32 𝑥2 + 𝑎33 𝑥3 = 𝑏3
can be expressed as [𝐴]{𝑥} = {𝑏} where [A] is the matrix of coefficients:
𝑎11 𝑎12 𝑎13
[𝐴] = [𝑎21 𝑎22 𝑎23 ]
𝑎31 𝑎32 𝑎33
{b} is the column vector of constants: {b}T = [b1 b2 b3 ] and {x} is the column vector of unknowns:
{𝑥}𝑇 = [𝑥1 𝑥2 𝑥3 ]

Page 3 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
3 Solving Small Numbers Of Equations
Before proceeding to the computer methods, we will describe several methods that are appropriate for
solving small (𝑛 ≤ 3 ) sets of simultaneous equations and that do not require a computer. These are the
graphical method, Cramer's rule and elimination of unknowns.
3.1 Graphical Method:
A graphical solution is obtainable for two equations by plotting them on two axis plot with on axis
corresponding to x1 and the other to x2. Because we are dealing with linear systems, each equation is a
straight line. x2

EX: use the graphical method to solve 10

3𝑥1 + 2 𝑥2 = 18
8
−𝑥1 + 2𝑥2 = 2 3x
1+
2x
Solution: 6 2 =
18 Solution: x1=4; x2=3
3 1
𝑥2 = − 𝑥1 + 9 𝑥2 = 𝑥1 + 1 4
2 2
The solution is the intersection of the two lines at 𝑥1 = 4 2 2 x2 =2
-x 1+
and 𝑥2 = 3. This result can be checked by substituting
these values into the original equations 0
0 x1
1 2 3 4 5 6
3.2 Cramer's Rule:
Cramer's rule is another solution technique that is best suited to small numbers of equations. Cramer's rule;
this rule states that each unknown in a system of linear algebraic equations may be expressed as of two
determinants with denominator D and with the numerator obtained from D by replacing the column of
coefficients of the unknown in equation by the constants b1, b2, . . . ,bn. for example, in a system of
equations:
𝑎11 𝑥1 + 𝑎12 𝑥2 + 𝑎13 𝑥3 = 𝑏1
𝑎21 𝑥1 + 𝑎22 𝑥2 + 𝑎23 𝑥3 = 𝑏2
𝑎31 𝑥1 + 𝑎32 𝑥2 + 𝑎33 𝑥3 = 𝑏3
𝑥1 , 𝑥2 and 𝑥3 would be computed as
𝑏1 𝑎12 𝑎13 𝑎11 𝑏1 𝑎13 𝑎11 𝑎12 𝑏1
𝑏
| 2 𝑎22 𝑎23 | 𝑎
| 21 𝑏2 𝑎23 | 𝑎
| 21 𝑎22 𝑏2 |
𝑎11 𝑎12 𝑎13
𝑏3 𝑎32 𝑎33 𝑎31 𝑏3 𝑎33 𝑎31 𝑎32 𝑏3
𝑥1 = , 𝑥2 = , 𝑥1 = where 𝐷 = |𝑎21 𝑎22 𝑎23 |
𝐷 𝐷 𝐷 𝑎31 𝑎32 𝑎33
3.3 Elimination of Unknowns
The elimination of unknowns by combining equations is an algebraic approach that can be illustrated for a
set of two equations:
𝑎11 𝑥1 + 𝑎12 𝑥2 = 𝑏1
𝑎21 𝑥1 + 𝑎22 𝑥2 = 𝑏2
The basic strategy is to multiply the equations by constants so that one of the unknowns will be eliminated
when the two equations are combined. The result is a single equation that can be solved for the remaining
Page 4 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
unknown. This value can then be substituted into either of the original equations to compute the other
variable. For example, the first equations might be multiplied by 𝑎21 and second by 𝑎11 to give:
𝑎21 𝑎11 𝑥1 + 𝑎21 𝑎12 𝑥2 = 𝑎21 𝑏1
𝑎11 𝑎21 𝑥1 + 𝑎11 𝑎22 𝑥2 = 𝑎11 𝑏2
which can be solved for
𝑏 𝑎12 𝑎 𝑏1
| 1 | 𝑎11 𝑏2 − 𝑎21 𝑏1 | 11 | 𝑎22 𝑏1 − 𝑎21 𝑏2
𝑏2 𝑎22 𝑎21 𝑏2
𝑥2 = 𝑎 = , 𝑥1 = 𝑎 =
11 𝑎12 𝑎11 𝑎22 − 𝑎21 𝑎12 11 𝑎12
|𝑎 | |𝑎 | 𝑎11 𝑎22 − 𝑎21 𝑎12
21 𝑎22 21 𝑎22

4 Gauss Elimination
The elimination of unknown was used to solve simultaneous equations. The procedure consisted of two
steps:
1- The equations were manipulated to eliminate one of the unknowns from the equations. The result of
this elimination step was that we had one equation with one unknown.
2- This equation could be solved directly and the result back-substituted into one of the original
equations to solve for the remaining unknown.
A general set of n equations:
𝑎11 𝑥1 + 𝑎12 𝑥2 + ⋯ + 𝑎1𝑛 𝑥𝑛 = 𝑏1 (1)
𝑎21 𝑥1 + 𝑎22 𝑥2 + ⋯ + 𝑎2𝑛 𝑥𝑛 = 𝑏2 (2)
⋮
𝑎𝑛1 𝑥1 + 𝑎𝑛2 𝑥2 + ⋯ + 𝑎𝑛𝑛 𝑥𝑛 = 𝑏𝑛 (𝑛)
Forward Elimination of Unknown:
The first phase is designed to reduce the set of equations
to an upper triangular system. The initial step will be to
eliminate the first unknown 𝑥1 , from the second through
the nth equations. To do this, multiply Eq.(1) by a21 / a11
to give:
a 21 a a
a 21 x1  a 21 x 2  ...........  21 a1n x n  21 b1
a11 a11 a11
Now, this equation can be subtracted from Eq.(2) to
give:

 a   a  a
 a 22  21 a12  x 2  ............   a 2 n  21 a1n  x n  b2  21 b1 '
or a 22 x 2  ...........  a 2' n x n  b2'
 a11   a11  a11

The same procedure is then repeated for the remaining equations until we get the upper triangle matrix.

Back Substitution:
Page 5 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems

bnn 1
𝑥𝑛 can be determined as follows: x n   n 1
the procedure, which repeated to evaluate the remaining x's
a nn
can be represented by the following formula:
𝑛
𝑏𝑖𝑖−1 − ∑ 𝑖−1
𝑎𝑖𝑗 𝑥𝑗
𝑗=𝑖+1
𝑥𝑖 = for 𝑖 = 𝑛 − 1, 𝑛 − 2, ⋯ ,1
𝑎𝑖𝑖𝑖−1
The logarithm for the elimination and back substitution is as follow:
a- forward elimination: b- Back substitution:
DO k =1,n-1 xn = bn /an,n
DO i = k+1,n DO i = n-1, 1, -1
factor = ai,k / ak,k Sum = 0
DO j = k+1, n DO j = i+1 , n
ai,j = ai,j – factor * ak,j Sum = sum + ai,j . xj
ENDDO ENDDO
bi = bi – factor * bk Xi = (bi – sum)/ai,i
ENDDO ENDDO
ENDDO

EX: use Gauss Elimination to solve:

3x1  0.1x 2  0.2 x3  7.85
0.1x1  7 x 2  0.3x3  19.3 carry six significant figures during the computation.
0.3x1  0.2 x 2  10 x3  71.4
Solution:
“Kept as exercise in class”

4.1 Pitfalls of Elimination Method:

Whereas, there are many systems of equations that can be solved with Gauss elimination, there are some
pitfalls that must be explored before writing a general computer program to implement the method.

4.1.1 Division by Zero:

During both the elimination and the back substitution phases, it is possible that a division by zero can occur.
For example, solve these equations:
2𝑥2 + 3𝑥3 = 8
4𝑥1 + 6 𝑥2 + 7𝑥3 = −3
2𝑥1 + 𝑥2 + 6𝑥3 = 5
The first row would involve division by 𝑎11 = 0. Problems also can arise when a coefficient is very close
to zero. The technique of pivoting has been developed to partially avoid these problems.

Page 6 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
4.1.2 Ill-Condition Systems:
The adequacy of the solution depends on the condition of the system. The following graph depicts three
cases that can pose problems when solving sets of linear equations.
x2 1 x2
x2=
1 x2 = 1.1
+ x 2= + x 1+
x 2
- 0.5x 1 -0.5x 1 - 0.46 1
2 x2 =
0.5 2x 2 = +
x + x2= -x 1+ - 0.5x 1
-0.5 1
x1 x1 x1
No solution Infinite solutions Ill conditioned system
First one shows the case when two equations represent parallel lines. For such situations, there is no solution
because the lines never cross. Second depicts the case where the two lines are coincident. For such situation,
there are an infinite number of solutions. Both types of systems are said to be singular. In addition, systems
that are very close to being singular can also cause problems. These systems are said to be ill-conditioned.
Well-conditioned systems are those where a small change in one or more of the coefficients results in a
similar small change in the solution. Ill-conditioned systems are those where small change in coefficients
result in large change in solution.
EX: solve the following system:
𝑥1 + 2𝑥2 = 10
1.1𝑥1 + 2𝑥2 = 10.4
Then solve it again, but with the coefficient of 𝑥1 in the second equation modified to 1.05.
210  210.4 110.4  1.110
Solution: x1  4 x2  3
12  21.1 12  21.1
With the slight change of the coefficient 𝑎12 from 1.1 to 1.05the result is changed dramatically:
210  210.4 110.4  1.110
x1  8 x2  1
12  21.05 12  21.05
4.2 Techniques For Improving Solutions:
The following techniques can be incorporated into Gauss elimination algorithm to overcome some of the
pitfalls discussed in the previous section:
4.2.1 Use of More Significant Figures:
The simplest remedy for ill-conditioning is to use more significant figures in the computation.
4.2.2 Pivoting
Obvious problems occur when a pivot element is zero because the normalization step leads to division by
zero. Problems may also arise when the pivot element is close to zero because if the magnitude of the pivot
element is small compares to the other element, then round-off error can be introduced. Therefore, before
each row is normalized, it is advantageous to determine the largest available coefficient in the column below
the pivot element. The rows can then be switched so that the largest element is the pivot element. The simple
algorithm of the pivoting is as follows:
Page 7 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
p = k
big =│ak,k│
DO ii = k+1 , n
dummy = │aii,k│
IF( dummy > big)
big = dummy
p = ii
END IF
END DO
IF( p≠ k)
DO jj = k,n
dummy = ap,jj
ap,jj = ak,jj
ak,jj = dummy
END DO
dummy = bp
bp = bk
bk = dummy
END IF
5 LU Decomposition:
One motive for introducing LU decomposition is that it provides an efficient means to compute the matrix
inverse. The inverse has a number of valuable applications in engineering practice.
As described in the previous section, Gauss elimination is deigned to solve systems of linear algebraic
equations, [𝐴]{𝑋} = {𝐵}. LU decomposition methods separate the time-consuming elimination of the matrix
[A] from manipulations of the right-hand side {B}. Thus, once [A] has been decomposed multiple right-
hand-side vectors can be evaluated in an efficient manner.
5.1 Overview of LU Decomposition:
Just as was the case with Gauss elimination, LU decomposition requires pivoting to avoid division by zero.
[𝐴] {𝑋} − {𝐵} = 0
Suppose that this equation could be expressed as an upper triangular system:
u11 u12 u13   x1   d1 
0 u    
 22 u 23   x 2   d 2 
 0 0 u 33     
 x 3  d 3 
Recognize that this is similar to the manipulation that occurs in the first step of Gauss elimination. The
previous equation can be expressed as:
[𝑈]{𝑋} − {𝐷} = 0 − − − − − − − −(1)
Now, assume that there is a lower diagonal matrix with 1's on the diagonal,
1 0 0
L  l12 1 0 .
l 31 l 32 1

This has property :[𝐿]{[𝑈]{𝑋} − {𝐷}} = [𝐴] {𝑋} − {𝐵} . If this equation holds, it follows from the rules
for matrix multiplication that [𝐿][𝑈] = [𝐴] and [𝐿]{𝐷} = [𝐵] -----(2)
A two-step strategy for obtaining solution can be based on:
Page 8 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
1- LU decomposition step:
[A] is factored into lower [L] and upper [U] triangular matrices.
2- Substitution step:
[L] and [U] are used to determine a solution {X} for right hand side
{B}. This step itself consists of two steps. First Eq.(2) is used to
generate an intermediate vector {D} by forward substitution. Then
the result is substituted into Eq.(1) which can be solved by back
substitution for {X}.
now, let's show how Gauss elimination can be implemented in this
way.
a11 a12 a13 
U    0 '
a 22 ' 
a 23 
which is the desired upper triangular format.
 0 0 a 33 
''

1 0 0 '
L   f 21 1 0 where the f 21  21
a
,
a
f 31  31 and f 32
a32
 '
a11 a11 a 22
 f 31 f 32 1

We could check out that the L U   A

EX: Derive an LU decomposition based on the Gauss elimination for:
 3  0.1  0.2  7.85 
A  0.1 7  0.3 and B   19.3
 
0.3  0.2 10   71.4 

Solution: after forward elimination, the following upper triangular matrix was obtained:
3  0.1  0.2 
First step of forward elimination the matrix A will be : 0 7.00333  0.293333

0  0.19 10.02 

3  0.1  0.2 
Second step : U   0 7.00333  0.293333

0 0 10.012 

Therefore,
 1 0 0
  0.19
L  0.033333 0 where f 21 
0.1 0.3
1  0.033333 , f 31   0.1 and f 32   0.02713
3 3 7.00333
 0.1  0.02713 1

Is L U   A ?

3  0.1  0.2   x1   7.85 

   
The forward elimination leads to 0 7.00333  0.293333  x 2    19.5617

0 0 10.012    
 x3   70.0843 


Page 9 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems

 1 0 0  d 1   7.85 
0.0333333    
 1 0 d 2    19.3 This will lead to
 0.1  0.02713 1   
d 3   71.4 


d1  7.85
0.0333333d1  d2   19.3
0.1d 1`  0.02713d 2  d3  71.4

 7.85 
 
𝑑1 = 7.85 , 𝑑2 = −19.5617 and 𝑑3 = 70.0843  D   19.5617
 70.843 
 

3  0.1  0.2   x1   7.85 

U X   D  0  0 7.00333  0.293333  x2    19.5617
 which can be solved by back
0 0 10.012    
 x3   70.0843 


substitution
 3 
X     2.5 
7.0003
 
6 The Matrix Inverse:
If matrix [A] is square, there is another matrix [A]-1 called the inverse of [A], for which [A]-1 [A]= [A] [A]-
1
=[I].
6.1 Calculating The Inverse:
The inverse can be computed in a column-by-column fashion by generating solutions with unit vectors as
the right- constants hand-side. For example, if the right-hand-side constant has a 1 in the first position and
zeros elsewhere,
The inverse can be computed in a column-by-column fashion by generating solutions with unit vectors as
the right- constants hand-side. For example, if the right-hand-side constant has a 1 in the first position and
1
 
zeros elsewhere, b  0 , the resulting solution will be the first column of the matrix inverse and so on.
0
 

 3  0.1  0.2
EX: employ LU decomposition to determine the matrix inverse for the system A  0.1 7  0.3
0.3  0.2 10 

Recall that the decomposition resulted in the following lower and upper triangular matrices:
3  0.1  0.2   1 0 0
U   0 7.00333  0.293333
 
L  0.033333 1 0
0 0 10.012   0.1  0.02713 1

Solution: the first column of the matrix inverse can be determined by performing the forward–substitution
solution procedure with unit vector as the right hand side vector.
Page 10 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems

 1 0 0  d1  1
0.0333333    
 1 0 d 2   0 and solved with forward substitution for [D]T = [1 -0.0333 -
 0.1  0.02713 1    
d 3  0
0.1009 ]. This vector can then be used as the right hand of the equation:
3  0.1  0.2   x1   1 
0 7.00333  0.293333  x    0.0333
 2   

which can be solved by back substitution
0 0 10.012   x3   0.1009
   


 0.33249   0.33249 0 0
X    0.00518 A
1
  0.00518 0 0
 0.01008  0.01008 0 0

To determine the second column:

 1 0 0  d1  0
0.0333333    
 1 0 d 2   1 This can be solved for [D] and the result can be used to determine
 0.1  0.02713 1    
d 3  0
{X}T = [0.004944 0.142903 0.00271] and the same procedure will be followed to determine he third
 0.33249 0.004944 0.006798
column of the inverse matrix. The resulted matrix is: A   0.00518 00.142903 0.004183

1

 0.01008 00.00271 0.09988 

7 Iterative Methods for Systems of Equations

7.1 Cholesky Decomposition:
Symmetric matrix is one where 𝑎𝑖𝑗 = 𝑎𝑗𝑖 for all 𝑖 and 𝑗. In other words, [𝐴] = [𝐴]𝑇 . Such systems occur
commonly in both mathematical and engineering problem contexts. They offer computational advantages
because only half the storage is needed and in most cases only half the computation time is required for their
solution. One of the most popular approaches involves Colesky decomposition. This algorithm is based on
the fact that a symmetric matrix can be decomposed as in [𝐴] = [𝑈]𝑇 [𝑈]. That is the resulting triangular
factors are the transpose of each other. The result can be simply expressed by the iterative relation. For
the 𝑖 𝑡ℎ row:

𝑖−1
2
𝑎𝑖𝑗 − ∑𝑖−1
𝑘=1 𝑢𝑘𝑖 𝑢𝑘𝑗
𝑢𝑖𝑖 = √𝑎𝑖𝑖 − ∑ 𝑢𝑘𝑖 ____________ (1) , 𝑢𝑖𝑗 = for 𝑗 = 𝑖 + 1, … … . . , 𝑛 __________(2)
𝑢𝑖𝑖
𝑘=1

 6 15 55 
EX: apply Cholesky decomposition to the symmetric matrix A  15 55 225 .
55 225 979

Solution:
for the first row (𝑖 = 1) Eq.(1) is employed to compute:
𝑢11 = √𝑎11 = √6 = 2.44949

Page 11 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
Eq.(2) is employed to compute:
𝑎12 15
𝑢12 = = = 6.1237
𝑢11 2.4494
𝑎13 55
𝑢13 = = = 22.454
𝑢11 2.4494
For the second row (𝑖 = 2):

2
𝑢22 = √𝑎22 − 𝑢12 = √55 − (6.1237)2 = 4.1833

𝑎23 − 𝑢12 𝑢13 255 − 6.123724 ∗ 22.454

𝑢23 = = = 20.9165
𝑢22 2.4494
For the second row (𝑖 = 3):

2 2
𝑢33 = √𝑎33 − 𝑢13 − 𝑢23 = 0√979 − (22.454)2 − (20.9165)2 = 6.110101

2.44949 6.123724 22.45366

Thus the Cholesky decomposition yields [𝑈} = [ 0 4.1833 20.9165 ] .Note [𝐴] = [𝑈]𝑇 [𝑈]
0 0 6.110101
.
The algorithm of the Cholesky decomposition:
DO k = 1,n
DO i = 1,k-1
sum = 0
DO j = 1,i-1
Sum = sum + aij . akj
END DO
aki = ( aki – sum ) / aii
END DO
sum=0
DO j = 1, k-1
2
sum = sum + akj
END DO
akk  akk  sum
END DO
7.3 Gauss-Seidel:
Iterative or approximate methods provide an alternative to the elimination methods described to this point.
Such approaches are similar to the techniques we developed to obtain the roots of a single equation. Those
approaches consisted of guessing a value and then using symmetric method to obtain a refined estimate of
the root. Because the present part deals with a similar problem- obtaining the values that simultaneously
satisfy a set of equations- we might suspect that such approximate methods could be useful in this context.
The Gauss-Seidel method is the most commonly used iterative method. Assume that we are given a set of n
equations: [𝐴] {𝑋} = {𝐵}.
Suppose that we are limit out selves to a 3 × 3 set of equations. If the diagonal elements are all nonzer, the
first equation can be solved for 𝑥1 the second for 𝑥2 and the third for 𝑥3 to yield :

Page 12 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
𝑗−1 𝑗−1 𝑗 𝑗−1 𝑗 𝑗
𝑗 𝑏1 − 𝑎12 𝑥2 − 𝑎13 𝑥3 𝑗 𝑏2 − 𝑎21 𝑥1 − 𝑎23 𝑥3 𝑗 𝑏3 − 𝑎31 𝑥1 − 𝑎32 𝑥2
𝑥1 = ; 𝑥2 = ; 𝑥3 =
𝑎11 𝑎22 𝑎33
where j and j − 1 are the present and previous iterations.
Now, we can start the solution process by choosing guesses for the x's. A simply way to obtain initial
guesses is to assume that they are all zero. These zeros can be substituted into the previous equations, which
can be used to calculate a new value for 𝑥1 = 𝑏1 /𝑎11 . Then we substitute this new value of 𝑥1 along with
the previous guess of zero for 𝑥3 to compute a new value for 𝑥2 . The process is repeated for third equation to
calculate a new estimate for 𝑥3 . Then we return to the first equation and repeat the entire procedure until our
solution converges closely enough to the true value. Convergence can be checked using the criterion
𝑗 𝑗−1
𝑥 −𝑥
𝜀𝑎,𝑖 =| 𝑖 𝑗𝑖 | × 100 ≤ 𝜀𝑠
𝑥𝑖
For all 𝑖, where 𝑗 and 𝑗 − 1 are the present and previous iterations.

EX: use the Gauss-Seidel method to obtain the solution of the same system used before:
3x1  0.1x 2  0.2 x3  7.85
0.1x1  7 x 2  0.3x3  19.3
0.3x1  0.2 x 2  10 x3  71.4
recall that the true solution is x1 = 3, x2 = -2.5 and x3 = 7
Solution: first solve each of the equation for its unknown on the diagonal:
7.85  0.1x 2  0.2 x3  19.3  0.1x1  0.3x3 71.4  0.3x1  0.2 x 2
x1  , x2  and x3 
3 7 10
7.85  0  0
By assuming that 𝑥2 𝑎𝑛𝑑 𝑥3 are zero, they can be used to compute x1   2.616667
3
This value with the assumed value of 𝑥3 = 0 can be substituted into second equation to calculate
 19.3  0.1(2.616667)  0
x2   -2.794524
7
The first iteration is completed by substituting the calculated values for 𝑥1 𝑎𝑛𝑑 𝑥2 into third equation to
yield
71.4  0.3(2.616667)  0.2(2.794524)
x3   7.005610
10
For second iteration, the same process is repeated to compute:
7.85  0.1(2.794524)  0.2(7.00561)
x1   2.990557  t  0.31%
3
 19.3  0.12.990557   0.37.005610 
x2   2.499625  t  0.015%
7
71.4  0.32.990557   0.2 2.499625
x3   7.000291  t  0.0042%
10

Page 13 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
The method is therefore converging on the true solution. Additional iteration could be applied to improve
the answers. In an actual problem we would not know the true answer. Consequently, we could estimate the
2.990557  2.616667
absolute error for example  a ,1  100%  12.5% . For 𝑥2 𝑎𝑛𝑑 𝑥3 , the error estimates
2.990557

are  a,2  11.8% and  a,3  0.076% .

7. Jacobi Method:
The difference between the Jacobi and Gauss-Seidel method is that in Jacobi method, we compute a set of
new x's on the basis of a set of old x's. Thus as new values are generated they are not immediately used but
they are retained for next iteration.
Algorithm for Gauss-Seidel:
DO i = 1,n
dummy = ai,i
DO j = 1 , n
ai,j = ai,j/dummy
END DO
bi = bi/dummy
END DO
DO i = 1 , n
sum = bi
DO j = 1 , n
IF I ≠ j THEN
sum = sum – ai,j* xj
END DO
xi = sum
END DO

Page 14 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
8. NONLINEAR SYSTEMS
The following is a set of two simultaneous nonlinear
equations with two unknowns:
𝑥12 + 𝑥1 𝑥2 = 10
𝑥2 + 3𝑥1 𝑥22 = 57
In contrast to linear systems which plot as straight lines,
these equations plot as curves on an 𝑥2 versus 𝑥1 graph.
As in the following figure, the solution is the intersection
of the curves.

8.1. Successive Substitution

A simple approach for solving a system of nonlinear equations is to use the same strategy that was employed
for fixed-point iteration and the Gauss-Seidel method. That is, each one of the nonlinear equations can be
solved for one of the unknowns. These equations can then be implemented iteratively to compute new
values which (hopefully) will converge on the solutions. This approach, which is called successive
substitution, is illustrated in the following example.
Example:
Use successive substitution to determine the roots of:
𝑥12 + 𝑥1 𝑥2 − 10 = 0
𝑥2 + 3𝑥1 𝑥22 − 57 = 0
Note that a correct pair of roots is 𝑥1 = 2 𝑎𝑛𝑑 𝑥2 = 3. Initiate the computation with guesses of 𝑥1 =
1.5 and 𝑥2 = 3.5.
Solution: 10 November 2022
The previous equations can be solved :
10 − 𝑥12
𝑥1 = ; 𝑥2 = 57 − 3𝑥1 𝑥22
𝑥2
The new value of the 𝑥1 and𝑥2 :
10 − 1.52
𝑥1 = = 2.21429 and 𝑥2 = 57 − 3 ∗ 2.1429 ∗ 3.52 = −24.37516
3.5
Thus, the approach seems to be diverging. This behavior is even more pronounced on the second iteration:
10 − 2.214292
𝑥1 = = −0.20910 and 𝑥2 = 57 − 3 ∗ (−0.20910) ∗ (−24.37516)2 = 429.709
−24.37516
Now we will repeat the computation but with the original equations set up in a different format. For
example, an alternative solution of two equations is:

57 − 𝑥2
𝑥1 = √10 − 𝑥1 𝑥2 ; 𝑥2 = √
3𝑥1

The new value of the 𝑥1 and𝑥2 :

Page 15 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems

57 − 3.5
𝑥1 = √10 − 1.5 ∗ 3.5 = 2.17945 ; 𝑥2 = √ = 2.86051
3 ∗ 2.17945

57 − 2.86051
𝑥1 = √10 − 2.17945 ∗ 2.86051 = 1.94053 ; 𝑥2 = √ = 3.04955
3 ∗ 1.94053

Thus, the approach is converging on the true values of 𝑥1 = 2 𝑎𝑛𝑑 𝑥2 = 3.

The previous example illustrates the most serious shortcoming of successive substitution that is,
convergence often depends on the manner in which the equations are formulated. Additionally, even in those
instances where convergence is possible, divergence can occur if the initial guesses are insufficiently close
to the true solution. These criteria are so restrictive that fixed-point iteration has limited utility for solving
nonlinear systems.

8.2. Newton-Raphson
Just as fixed-point iteration can be used to solve systems of nonlinear equations, other open root location
methods such as the Newton-Raphson method can be used for the same purpose. Recall that the Newton-
Raphson method was predicated on employing the derivative (i.e., the slope) of a function to estimate its
intercept with the axis of the independent variable that is, the root. We used a graphical derivation to
compute this estimate. An alternative is to derive it from:
𝑓(𝑥𝑖+1 ) = 𝑓(𝑥𝑖 ) + (𝑥𝑖+1 − 𝑥𝑖 )𝑓 ′ (𝑥𝑖 )
where 𝑥𝑖 is the initial guess at the root and 𝑥𝑖+1 is the point at which the slope intercepts the 𝑥 axis. At this
intercept, 𝑓 (𝑥𝑖+1 ) by definition equals zero and the previous equation can be rearranged to yield
𝑓(𝑥𝑖 )
𝑥𝑖+1 = 𝑥𝑖 −
𝑓 ′ (𝑥𝑖 )
The multiequation form is derived in an identical fashion. However, a multivariable Taylor series must be
used to account for the fact that more than one independent variable contributes to the determination of the
root. For the two-variable case, a first-order Taylor series can be written for each nonlinear equation as:
𝜕𝑓1,𝑖 𝜕𝑓1,𝑖
𝑓1,𝑖+1 = 𝑓1,𝑖 + (𝑥1,𝑖+1 − 𝑥1,𝑖 ) + (𝑥2,𝑖+1 − 𝑥2,𝑖 )
𝜕𝑥1 𝜕𝑥2
𝜕𝑓2,𝑖 𝜕𝑓2,𝑖
𝑓2,𝑖+1 = 𝑓2,𝑖 + (𝑥1,𝑖+1 − 𝑥1,𝑖 ) + (𝑥2,𝑖+1 − 𝑥2,𝑖 )
𝜕𝑥1 𝜕𝑥2
𝜕𝑓2,𝑖 𝜕𝑓
𝑓1,𝑖 − 𝑓2,𝑖 1,𝑖
𝜕𝑥2 𝜕𝑥2
𝑥1,𝑖+1 = 𝑥1,𝑖 − − − − − − − − − − −(8.1)
𝜕𝑓1,𝑖 𝜕𝑓2,𝑖 𝜕𝑓1,𝑖 𝜕𝑓2,𝑖
× − ×
𝜕𝑥1 𝜕𝑥2 𝜕𝑥2 𝜕𝑥2
𝜕𝑓1,𝑖 𝜕𝑓
𝑓2,𝑖 − 𝑓1,𝑖 2,𝑖
𝜕𝑥1 𝜕𝑥1
𝑥2,𝑖+1 = 𝑥2,𝑖 − − − − − − − − − − −(8.2)
𝜕𝑓1,𝑖 𝜕𝑓2,𝑖 𝜕𝑓1,𝑖 𝜕𝑓2,𝑖
× − ×
𝜕𝑥1 𝜕𝑥2 𝜕𝑥2 𝜕𝑥2

Page 16 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
The denominator of each of last two equations is formally referred to as the determinant of the Jacobian of
the system.
Example: Use the multiple-equation Newton-Raphson method to determine roots of:
𝑥12 + 𝑥1 𝑥2 − 10 = 0
𝑥2 + 3𝑥1 𝑥22 − 57 = 0
Initiate the computation with guesses of 𝑥1 = 1.5 and 𝑥2 = 3.5.
Solution: First compute the partial derivatives and evaluate them at the initial guesses of 𝑥1 𝑎𝑛𝑑 𝑥2 :
𝜕𝑓1,0 𝜕𝑓1,0
= 2𝑥1 + 𝑥2 = 2 ∗ 1.5 + 3.5 = 6.5 ; = 𝑥1 = 1.5
𝜕𝑥1 𝜕𝑥2
𝜕𝑓2,0 𝜕𝑓2,0
= 3𝑥22 = 3 ∗ 3.52 = 36.75 ; = 1 + 6𝑥1 𝑥2 = 1 + 6 ∗ 1.5 ∗ 3.5 = 3.25
𝜕𝑥1 𝜕𝑥2
Thus, the determinant of the Jacobian for the first iteration is 6.5 ∗ 3.25 + 1.5 ∗ 36.75 = 156.1250.
The values of the functions can be evaluated at the initial guesses as:
𝑓1,0 = 1.52 + 1.5 ∗ 3.5 − 10 = −2.5
𝑓2,0 = 3.5 + 3 ∗ 1.5 ∗ 3.52 − 57 = 1.625
These values can be substituted into Eq. (8.1&8.2) to give:

−2.5 ∗ 3.25 − 1.625 ∗ 1.5

𝑥1 = 1.5 − = 2.03603
156.125
1.625 ∗ 6.5 − (−2.5) ∗ 36.75
𝑥2 = 3.5 − = 2.84388
156.125

Thus, the results are converging to the true values of 𝑥1 = 2 and 𝑥2 = 3. The computation can be repeated
until an acceptable accuracy is obtained.

The two-equation Newton-Raphson approach can be generalized to solve n simultaneous equations. To do

this, the equations can be written for the kth equation as:

𝜕𝑓𝑘,𝑖 𝜕𝑓𝑘,𝑖 𝜕𝑓𝑘,𝑖 𝜕𝑓𝑘,𝑖 𝜕𝑓𝑘,𝑖 𝜕𝑓𝑘,𝑖

𝑥1,𝑖+1 + 𝑥2,𝑖+1 + ⋯ + 𝑥𝑛,𝑖+1 = −𝑓𝑘,𝑖 + 𝑥1,𝑖 + 𝑥2,𝑖 + ⋯ + 𝑥𝑛,𝑖
𝜕𝑥1 𝜕𝑥2 𝜕𝑥𝑛 𝜕𝑥1 𝜕𝑥2 𝜕𝑥𝑛
where the first subscript k represents the equation or unknown and the second subscript denotes whether the
value or function in question is at the present value (i) or at the next value (𝑖 + 1). Notice that the only
unknowns in the previous equation are the 𝑥𝑘,𝑖+1 terms on the left-hand side. All other quantities are
located at the present value (i) and, thus, are known at any iteration. Consequently, the set of equations
generally represented by the previous equation (i.e., with k = 1, 2, . . . , n) represents a set of linear
simultaneous equations that can be solved numerically by the elimination methods elaborated in previous
sections.
Matrix notation can be employed to express the previous equation briefly as:
[𝐽]{𝑥𝑖+1 } = −{𝑓} + [𝐽]{𝑥𝑖 } − − − − − −(8.3)

Page 17 of 18
2 May 2023
 Linear and Nonlinear Algebraic Equation Systems
where the partial derivatives evaluated at i are written as the Jacobian matrix consisting of the partial
derivatives:
𝜕𝑓1,𝑖 𝜕𝑓1,𝑖 𝜕𝑓1,𝑖
⋯
𝜕𝑥1 𝜕𝑥2 𝜕𝑥𝑛
𝜕𝑓2,𝑖 𝜕𝑓2,𝑖 𝜕𝑓2,𝑖
[𝐽] = 𝜕𝑥1 ⋯
𝜕𝑥2 𝜕𝑥𝑛
⋮ ⋮ ⋮ ⋮
𝜕𝑓𝑛,𝑖 𝜕𝑓𝑛,𝑖 𝜕𝑓𝑛,𝑖
⋯
[ 𝜕𝑥1 𝜕𝑥2 𝜕𝑥𝑛 ]
The initial and final values are expressed in vector form as:
{𝑥𝑖 }𝑇 = [𝑥1,𝑖 𝑥2,𝑖 ⋯ 𝑥𝑛,𝑖 ] and {𝑥𝑖+1 }𝑇 = [𝑥1,𝑖+1 𝑥2,𝑖+1 ⋯ 𝑥𝑛,𝑖+1 ]
Finally, the function values at i can be expressed as {𝑓}𝑇 = [𝑓1,𝑖 𝑓2,𝑖 ⋯ 𝑓𝑛,𝑖 ].
Recall that the single-equation version of the Newton-Raphson method is:
𝑓(𝑥𝑖 )
𝑥𝑖+1 = 𝑥𝑖 −
𝑓 ′ (𝑥𝑖 )
If equation (8.3) is solved by multiplying it by inverse of the Jacobian, the result is:
{𝑥𝑖+1 } = {𝑥𝑖 } − [𝐽]−1 {𝑓}
Comparison of the last two equations clearly illustrates the parallels between the two equations. In essence,
the Jacobian is analogous to the derivative of a multivariate function.
Such matrix calculations can be implemented very efficiently in MATLAB. We can illustrate this by using
MATLAB to duplicate the calculations from the previous Example.” This is kept as exercise in class”
Problem: Apply Newton’s Method to find both solutions of the system of three equations:
2𝑥12 − 4𝑥1 𝑥22 + 3𝑥32 + 6𝑥3 + 2 = 0
𝑥12 + 𝑥22 − 2𝑥2 + 2𝑥32 − 5 = 0
3𝑥12 − 12𝑥1 + 𝑥22 + 3𝑥32 + 8 = 0

Page 18 of 18
2 May 2023