#12 (07/16/2025)

Gauss-Jordan elimination method

The linearity principle:

a11 x1 + a12x2 + …+a1nxn
=
c1
a21 x1 + a22x2 + …+a2nxn
=
c2,
(1)

(λ1 a11 + λ2 a12)x1+(λ1 a21 + λ2 a22)x2+…+(λ1 a1n + λ2 a2n)xn = λ1 c1 + λ2 c2.
(2)
Example

2x+3y+4z = 6
3x+5y+2z = 5
4x+3y+30z=32.
⎫
⎪
⎬
⎪
⎭
(3)

⇒

1x+0y+0z = ?
0x+1y+0z = ?
0x+0y+1z = ?.
⎫
⎪
⎬
⎪
⎭
(4)
Ref. line number x y z = Comment
(1) 2 3 4 6
(2) 3 5 2 5
(3) 4 3 30 32
(4) 1 1.5 2 3 (1)÷2
(5) 3 5 2 5
(6) 4 3 30 32
(7) 1 1.5 2 3
(8) 0 0.5 -4 -4 (5)−(4)×3
(9) 0 -3 22 20 (6)−(4)×4
(10) 1 1.5 2 3
(11) 0 1 -8 -8 (8)÷0.5
(12) 0 -3 22 20
(13) 1 0 14 15 (10)−(11)×1.5
(14) 0 1 -8 -8
(15) 0 0 -2 -4 (12)−(11)×(−3)
(16) 1 0 14 15
(17) 0 1 -8 -8
(18) 0 0 1 2 (15)÷(−2)
(19) 1 0 0 -13 (16)−(18)×14
(20) 0 1 0 8 (17)− (18)×(−8)
(21) 0 0 1 2
Finally we obtain
x = −13,     y=8,     z=2.
The number of multiplications required for the Gauss-Jordan elimination method is about n3. For n=10, this amounts to 1,000. Compare this number with the one required for Cramer's rule.

Inverse matrix by the Gauss-Jordan elimination method (optional)

Apply the Gauss-Jordan elimination method for
( A | I)
(5)
i.e.

⎛
⎜
⎜
⎜
⎝
2
3
4
1
0
0
3
5
2
0
1
0
4
3
32
0
0
1
⎞
⎟
⎟
⎟
⎠ 		→⎛
⎜
⎜
⎜
⎝
1
2/3
2
1/2
0
0
0
1/2
−4
−3/2
1
0
0
−3
22
−2
0
1
⎞
⎟
⎟
⎟
⎠ 		→⎛
⎜
⎜
⎜
⎝
1
0
14
5
−3
0
0
1
−8
−3
2
0
0
0
−2
−11
6
1
⎞
⎟
⎟
⎟
⎠ 		→⎛
⎜
⎜
⎜
⎝
1
0
0
−72
39
7
0
1
0
41
−22
−4
0
0
1
11/2
−3
−1/2
⎞
⎟
⎟
⎟
⎠
(6)
We thus obtain
A−1 = ⎛
⎜
⎜
⎜
⎝
−72
39
7
41
−22
−4
11/2
−3
−1/2
⎞
⎟
⎟
⎟
⎠
(7)

Implementation

#include <stdio.h>
#define N 3
int main()
{
 double a[N][N+1]={{2, 3, 4, 6},{3, 5, 2, 5},{4, 3, 30, 32}};
  double pivot,d;
  int i,j,k;

for(k=0; k<N; k++)
{
 pivot=a[k][k];

 for(j=k; j<N+1; j++) a[k][j]=a[k][j]/pivot;
 for(i=0; i<N;  i++)
  {
    if(i != k)
     {
       d=a[i][k]; 
        for(j=k; j<N+1; j++) a[i][j]=a[i][j]-d*a[k][j];
     }
   }
}

for(i=0; i<N; i++) printf("x[%d]=%lf\n", i+1, a[i][N]);
return 0;
}

LU decomposition (optional)

Any matrix, A, can be uniquely factorized as

A = L U,
(8)
where

L = ⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
1
0
0
…
0
l21
1
0
…
0
l31
l32
1
…
0
:
:
:
···
:
ln1
ln2
ln3
…
1
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠ 		,   U = ⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
u11
u12
u13
…
u1n
0
u22
u23
…
u2n
0
0
u33
…
u3n
:
:
:
···
:
0
0
0
…
unn
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠ 		.
(9)
The matrix, L, is a lower triangular matrix and the matrix, U, is an upper triangular matrix. Note that the diagonal elements of L are set to be 1. This decomposition is called LU decomposition (or LU factorization) and provides an effective way of solving simultaneous equations which is more efficient than the Gauss-Jordan elimination method.
The decomposition above is unique as it is possible to find L and U directly as

⎛
⎜
⎜
⎜
⎜
⎜
⎝
1
0
0
0
l21
1
0
0
l31
l32
1
0
l41
l42
l43
1
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		⎛
⎜
⎜
⎜
⎜
⎜
⎝
u11
u12
u13
u14
0
u22
u23
u24
0
0
u33
u34
0
0
0
u44
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		= ⎛
⎜
⎜
⎜
⎜
⎜
⎝
a11
a12
a13
a14
a21
a22
a23
a24
a31
a32
a33
a34
a41
a42
a43
a44
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		.
(10)
Writing each element explicitly gives

⎧
⎪
⎪
⎨
⎪
⎪
⎩
u11 = a11,
u12=a12,
u13=a13,
u14=a14
l21 u11 = a21,
l21 u12+u22 = a22,
l21 u13+u23 = a23,
l21 u14+u24 = a24
l31u11 = a31,
l31u12+l32u22 = a32,
l31u13+l32u23+u33 = a33,
l31u14+l32u24+u34 = a34
l41u11 = a41,
l41u12+l42u22 = a42,
l41u13+l42u23+l43u33 = a43,
l41u14+l42u24+l43u34+u44 = a44,
(11)
from which one can solve the unknowns, lij and uij, as

⎧
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎩
u11=a11,
u12=a12,
u13=a13,
u14=a14
l21= a21
u11
 ,
u22 = a22− l21u12,
u23=a23− l21 u13,
u24=a24−l21u14
l31= a31
u11
 ,
l32= a32− l31u12
u22
 ,
u33=a33− l31 u13− l32 u23,
u34=a34− l31 u14− l32 u24
l41= a41
u11
 ,
l42 = a42−l41u12
u22
 ,
l43= a43−l41u13−l42u23
u33
 ,
u44=a44−l41u14−l42u24−l43u34.
(12)
Note that the immediate results for uij and lij must be used to solve for the next result.
To solve
A x = L U x = c,
(13)
  1. Solve L y = c first.
  2. Then, solve U x = y.
Finding L and U from A requires less effort than finding A−1 and can accomplish the same goal as finding A−1.

⎛
⎜
⎜
⎜
⎜
⎜
⎝
1
0
0
0
l21
1
0
0
l31
l32
1
0
l41
l42
l43
1
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		⎛
⎜
⎜
⎜
⎜
⎜
⎝
y1
y2
y3
y4
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		= ⎛
⎜
⎜
⎜
⎜
⎜
⎝
c1
c2
c3
c4
⎞
⎟
⎟
⎟
⎟
⎟
⎠
(14)
This can be solved easily as

⎧
⎪
⎪
⎨
⎪
⎪
⎩
y1 = c1
l21 y1 + y2 = c2
l31 y1 + l32y2 + y3 = c3
l41 y1 +l42 y2 + l43y3 + y4 = c4
(15)
from which one obtains

⎧
⎪
⎪
⎨
⎪
⎪
⎩
y1 = c1
y2 = c2 −l21 y1
y3 = c3− l31 y1 − l32y2
y4 = c4 − l41 y1−l42 y2 − l43y3
(16)
For U x=y,

⎛
⎜
⎜
⎜
⎜
⎜
⎝
u11
u12
u13
u14
0
u22
u23
u24
0
0
u33
u34
0
0
0
u44
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		⎛
⎜
⎜
⎜
⎜
⎜
⎝
x1
x2
x3
x4
⎞
⎟
⎟
⎟
⎟
⎟
⎠ 		= ⎛
⎜
⎜
⎜
⎜
⎜
⎝
y1
y2
y3
y4
⎞
⎟
⎟
⎟
⎟
⎟
⎠
(17)
First

⎧
⎪
⎪
⎨
⎪
⎪
⎩
u11 x1 + u12 x2 + u13 x3 + u14 x4 = y1
u22 x2+u23x3 + u24x4 = y2
u33 x3 + u34x4 = y3
u44 x4=y4
(18)
from which one obtains

⎧
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎩
x4 = y4
u44
x3 = y3− u34x4
u33
x2 = y2−u23x3 − u24x4
u22
x1 = y1−u11 x1 − u12 x2 − u13 x3
u11
(19)
This procedure is called backward substitution.
It can be shown that the number of operations (multiplications, divisions) to decompose A into A = LU is about n3/3 and to obtain x from Ly=c and Ux=y, n2 is required so the total number of operations is about n3/3 which is 1/3 of the number required for the Gauss-Jordan elimination method.

Notes

  1. If A is a symmetric matrix (aij=aji), LU-decomposition is called the Cholesky decomposition. The number of operations required for the Cholesky decomposition is n3/6.
  2. If A is a sparse matrix, so are L and U.
  3. If A is a triangular diagonal matrix (typical of the finite element method), so are L and U.

Gauss-Seidel method (Important)

If the diagonal elements of the matrix, A, are generally larger than non-diagonal elements, the linear equation, A x = c, can be solved by an iterative method called the Gauss-Seidel method.
Example:

⎧
⎪
⎨
⎪
⎩
7x+y+2z
=
10
x+8y+3z
=
8
2x+3y+9z
=
6.
(20)

⎧
⎪
⎨
⎪
⎩
x
=
(10−y−2z)/7
y
=
(8−x−3z)/8
z
=
(6−2x−3y)/9.
(21)

⎧
⎪
⎨
⎪
⎩
xn+1
=
(10−yn−2zn)/7
yn+1
=
(8−xn−3zn)/8
zn+1
=
(6−2xn−3yn)/9,
(22)
with x0 = y0 = z0 = 0 as the initial approximation. This scheme is also called the Jacobi method. A slightly better iteration scheme is to use the most immediate values for the next approximation, i.e.

⎧
⎪
⎨
⎪
⎩
xn+1
=
(10−yn−2zn)/7
yn+1
=
(8−xn+1−3zn)/8
zn+1
=
(6−2xn+1−3yn+1)/9,
(23)
This scheme is called the Gauss-Seidel iterative method. It is easy to implement this method: 
#include <stdio.h>
int main()
{
double x, y, z;
int i,n;

x=y=z=0.0;

printf("Enter # of iteration = ");
scanf("%d", &n);

for (i=0;i<n;i++)
{
x = (10-y-2*z)/7;
y = (8-x-3*z)/8.0;
z = (6-2*x-3*y)/9.0;
}

printf("x = %lf, y= %lf, z=%lf\n", x,y,z);
return 0;
}


File translated from TEX by TTH, version 4.03.
On 15 Jul 2025, 22:55.