Computers & Structures Vol. 42, No. 5, pp. 84e852, Printed in Gnat Britain.
1992 0
0045-7949/92 $5.00 + 0.00 1992 Pergamon Press plc
TECHNICAL NOTE METHOD OF CHANGE OF VARIABLE DECOMPOSITION TO SOLVE AND INVERT A SYSTEM OF LINEAR EQUATIONS N. BENCHARIF and S. F. NC Department of Civil Engineering, University of Ottawa, Ottawa, Ontario, Canada KIN 9B4 (Received 9 January 1991) Abstract-Solving a system of simultaneous linear equations is perhaps one of the most common mathematical operations for all engineers and scientists. In this paper a new method to carry out this mathematical operation is presented using a change of variable technique. The technique has been well tested and proven to be extremely efficient both in terms of computer time and storage. Essentially, the method allows the user either to transform the system half-way and then solve it by back substitution or continue the operation and calculate the inverse of the system by inverting a diagonal matrix [D] and premultiplying it by a transformed matrix [Cl. The new method has been proven to be more efficient than the Gauss-Jordan or LU decomposition methods commonly employed to solve simultaneous linear equations.
1. INTRODUCTION
unknown vector {Y} using back substitution, i.e.
The solution of simultaneous equations can basically be accomplished either by direct or iterative techniques. In this paper we will investigate a new direct method based on change of variable as a solution to a system of simultaneous linear equations. The best known direct methods yield their solutions by forward elimination and backsubstitution, (e.g., Gauss, Gauss-Jordan, and Crout methods) involving the right-hand side of the system in the solution. The Cholesky decomposition method is more appropriate in solving a system of simultaneous equations when inversion is not required. If the inversion of the system is needed, both the Gauss-Jordan and Cholesky methods require the same number of operations. In this paper a new technique is presented to solve systems of simultaneous equations using change of variables. Consider a system of equations in its explicit form
1X]=
{X1 =
[aw’P).
(5)
Because [D] is a diagonal matrix, its inverse is nothing but the inverse of its elements, and the premultiplication of this matrix by the [C] matrix will yield [A]-‘, the inverse of the original matrix. The technique presented here results in one square matrix and one vector, and the matrix is inverted using the locations of the original matrix. The efficiency of the method is demonstrated with an approximate gain of 50% in computer time than the Gauss-Jordan method and 50% gain in computer storage than the LU method. The development of the method and the alternative ways of its application to solve systems of simultaneous equations will now be demonstrated.
2.
(2)
where [D] is a diagonal matrix. The second one will be a combination of the first unknown vector {X} and the new unknown vector {Y} where {X} will represent the upper triangular elements excluding the diagonal elements and {Y} will represent the lower triangular elements including the diagonals
r
(4)
By simple substitution of eqn (2) in (4) we get the equation
The technique will change the original system into two systems, the first one will be in terms of a new variable vector {Y}
PI{ YI = {BI,
[Cl{Y).
THE CHANGE
OF VARIABLE
METHOD
This”method consists of solving a system of equations using a new technique. The procedure is better described as the ‘technique of change of variables’. By changing the variables from x to y a diagonal system in terms of y would result. The solution of this system entails back substitution in the transformed system to solve for x. However, if back substitution is employed, the final matrix will be the inverted matrix of the system. All the steps in the mathematical operation will now be shown.
N
i=N
C %Yj j=1
xi=
I
i
1 ci,j.Yjf
j=t
1 j-i+
2.1. Solution of the system (3)
N
Ci,jXj,
For example, a system with three equations and three unknowns is chosen in the general form
i c N.
t
From this stage, solution can be performed without any further transformation. In the case where inversion is required, eqn (3) will be expressed only in terms of the new 849
Technical Note
850
Starting with the first equation, since the purpose is to have all the elements equal to zero except the diagonal term, we Put x, =y,
aI2
aI3
alI
all
--x2--x3
Similarly the second system may be formed from equations (6) (10) and (13) as
[Z]=[G;
;;I
E;Z].
and by making use of the first equation allYI = b,.
(7)
Considering the second equation and substituting the value of x, gives aI2 azl Y, - -x2 alI
aI3 - -x3 ai1
*22x2+ *23x3
+
(8)
>
[a{‘v
after collecting terms it becomes *,,Y,+(
where the c terms in the second system correspond to the constants of the equations defining the three equations for XI, ~2,~3. From the first system we solve for y and then we back substitute in the second system to determine the x vector. If we separate the variables we get the following expressions
= [Cl,{ YI
(19)
and
-$+*22)x2+(
--~+4r)xr.
To change the variable and have all the terms zero except the diagonal term *zl
[ - (*,3/*,,)
[Cl”{X)= P4Pl-‘P~.
(9)
+
*,,I
x2=-[-(*~2/*,,)+*221Y’+y2-[-(a,2/a,,)+a22]X3~
(10)
Replacing the second equation results in
(11) Now, consider the third equation, and replace the value of x, first then the value of x2. Once all the calculations are completed we can then obtain GYI + 532Y2+ 533x3.
(12)
(20)
From the above system the left-hand side matrix coefficients [Cl, corresponds to an upper triangular matrix with diagonal elements all equal to one, and the right-hand side coefficients [C],[D]-’ corresponds to the inverse of a lower triangular matrix. If the present technique is compared with the LU decomposition method presented by Press et al. [I], or the methods presented by Westlake [2], it can be observed that the two decompositions are totally different. 2.2. Inversion of the system
It should now be noted that the second system resulting from eqns (6), (10) and (13) is in terms of x and y. The final operation is to express the whole system in terms of y only. The first equation defining x, is in terms of x2 and x3. First by substituting the expression of x2 in the expression of x,, it becomes Y, + c,2(c2,y, +y2 + c23x3) + c,3x3
(21)
(c,2c21 + l)y, + c,2y2 + (c,2c23 + c,3)x3
(22)
To have all the terms zero except the diagonal term we put
XI [I[
41 x3= -T+YI-5Y2+Y3.
and the general system will reduce to Substituting x3 into the third equation gives C,,Y, = b3 1
x2
(14)
(c12c21+ l)Yl
‘12y2
(c,2c23 + c,3)x3
C2IYI
Y2
c23x3
C3IYI
c32y2
Y3
=
x3
where (*3i*l2
4, = a3, +
-
(15)
- aI2 + *22*1l
*3
a32 -
q2 =
*32)*21
a3, =
-
aI1
a32 *23+
-
(16) +(c12c23
43
+
*3,
aI3
all*u - aI2
+
*33*11
(17)
r*d9
1O 0
0
1
iiy31 =Ei].
(-%+,2)2
+
c!3)(c31~l
+
c32~2+~3)
(23)
and, after simplification it becomes
alI
Now we have the first system formed from equations (7) (11) and (14) 0
Substituting the x3 expression in x, and x2, row 1 becomes (c12c2, + l)Yl + Cl2Y2
62
aI1 a3142
1.
0
[l +
cl2c21 +
+
c3,(c,2c23
1%
+
+
c,3)b,
c32(cIZc23
+
c13)b2
+
(c12c23 +
c13)Y3’
(24)
Similarly, row 2 becomes c,,Y,
+Y2+c23(c31YI
(25)
+c32Y2+Y3).
After collecting terms we get
Writing it in condensed form we get cc21 +
[(c,2c21
+
I) +
c23c31h
+
c1 +
(18)
Pl{YI = PI.
The
c31 (c,2c23
(c2, + c23c31) c31
+
c,3)1
system now becomes
k,2 + c32(c12c23+ c,3)1
(c,2c23 + c,3)
(l + c23c32)
c23
c32
1
C23c32)Y2 +
c23Y3.
(26)
851
Technical Note 8y, - 25.6x2 + 19.2x, = 1
In explicit form this can be. written as (27)
{X] = [Cl{ I?.
(40)
then
Once all the terms are expressed in terms of y, we can calculate the y vector from the fint system
- 16
8y, - 25.
- 6.4
-7j-yYl+Y2-7773
.
>
l
+19.2-x,=
(41) - 28.571433, - 25.6~~+ 4.57143~~= 1.
(Y2 /=( 0 (-:+a+’ b’1 1
From this we put the value of xr equal to
0
0
(42)
- 25.6
- 28.57143
where in the condensed form we put
x3=
{Y} = [D]_‘(B).
-
(43)
Yl--Y2+Y’.
4.57143
(28) After we substitute it in the previous equation we get
Thus, the x vector will be equal to 4.57143yr = 1.
[CIP-‘{~h
WI =
(44)
(29)
From this relation we conclude that the product of [C][D]-’ is nothing but the inverse of [A] and its determinant is equal to the determinant of the diagonal matrix [D] [A]_’ = [C][D]_‘.
Then we have formed two systems the first system from eqns (33), (3% and (44)
(30)
2.3. Numerical application
And the second system from eqns (31), (36) and (43)
Consider a 3 x 3 system of equation -16 20
-16
-16
24 8
2 -xl
‘1
x,
= 1
20 _xr
-1
-8
-32
2
x=Yl-~x2-5jxl
- 16 x2= - xy,
- 6.4 +Y,---3
11.2
Starting with the first unknown, and by putting -
-16
2
x3 = -
xI=Yl-~x2-5jx3.
By substitution, three equations equation becomes
would result: the first
2
- 16x, + 2x, = 1
From above we determine the values of y then by applying eqn (3) we get the x values. To invert the system we have to replace the x, expression in the x, expression
(32) XI =Yl-T
2oy, = 1.
+24x2-8x,=1
(34)
>
2 --x3
>
20 (45)
(35) x2 =
- 16
- 6.4 11.2
+yz--xx,.
(36)
- 6.4 11.2
.
1.42857y, + y2 + 0.57143(6.25y, + 5.6~~+ yr)
(47)
x, = 5.00001y, + 4.20001y2+ 0.57143~~.
(48)
Then x, xl = 2.14286~’ + 0.8~~+ 0.35714(6.25y, + 5.6~~+ y,)
We substitute it in the previous equation giving -syl+yz--xj
(46)
Now we have to replace the xr expression in x2 and the new
From this expression we put the value of x2 equal to
(
11.2
Xl
-16y,+11.2x2-6.4x,=1.
-16y,+11.2
- 6.4
x, = 2.14286~~+ 0.8~~+0.35714x,.
y,-=$x2-ix,
EYI
7-5Yl+Y2--x’
(
(33)
The second equation becomes
x2= -
- 16
- 16
>
-1
28.57143 - 25.6 4.57143 Yl--Y2+Y3.
(49) -6.4x, = 1
x, = 4.37499~’ + 2.79998~~+ 0.35714~~.
(50)
(37) Then the system will look like
11.2y,= 1.
(38)
The third equation results after we substitute the value of xl and then the value of x, in the third equation -16
2 -
Yl--57x2-205 >
32x* + 20x, = 1 (39)
I[1 y,
5.00001 4.20001 0.57143 6.25 1.0 4.37499 5.6 2.79998 0.35714
y,
.
Y3
Then if we replace the first system in the second system we
852
Technical Note Table 1. Comnarison of the methods
get XI X2
[
x3
System of equations
Methods
60 x 60
Gauss-Jordan LU decomposition Present technique
5 4 4
60 x 60 2(60 x 60) 60 x 60
135 x 135
Gauss-Jordan LU decomposition Present technique
12 8 9
135 x 135 2(135 x 135) 135 x 135
240 x 240
Gauss-Jordan LU decomposition Present technique
52 29 32
240 x 240 2(240 x 240) 240 x 240
375 x 375
Gauss-Jordan LU decomposition Present technique
216 109 121
375 x 375 2(375 x 375) 375 x 375
= 1
Computer time (CPU set)
5.00001 4.20001 0.57143
6.25 .4.37499
5.6 1.0 2.79998 0.35714 I 0
0
11.2
0
‘20 0
X
-’
1
0 4.57143 I[
0
1 1
or we can write
5.00001 4.20001 0.57143 6.25 5.6 1.0 4.37499 2.79998 0.35714 1 r I *
20
X
0
0
i r
Oii3 O 0
0
1 ~ 4.57143
Thus the final system is
-1[ x, xl
.x3
=
The methods are based on triangularization [2] which transforms the original system A in general to LDU or the condensed form LU where the diagonal system belongs to L, CI,or both, but the transformation is represented only on the left-hand side. The technique presented in this paper transforms the original system in both the left- and righthand system [eqns (20)]. The solution of simultaneous equations without inversion using the new technique considers the same capacity in space and operation count with the former LU decomposition method. However, when the inversion of the system is required the two methods agree only on the number of operation count. The former LU decomposition considers an additional N x N matrix to perform and store the inverse whereas the new technique computes the inverse in place. A comparison between Gauss-Jordan, LU decomposition and the present technique to solve a system of simultaneous equations in terms of time computation and storage capacity is shown in Table 1. Both the LU decomposition technique and the technique presented in this paper share the advantage that the right hand side of the system is not considered in the transformation, and the solution can be performed for N right-hand side vectors. The advantage of the new technique on the former one in space storage is significant. The efficiency of the new technique to solve large systems has been tested up to 375 simultaneous equations. Comparing the new technique with the previous ones gives a gain of 50% in computer time with Gauss-Jordan method and about 50% in saving of computer storage with the regular LU decomposition method.
0.31250 0.25000 0.25000 0.21875 0.50000 0.07812 0.37500 0.21876 0.12500
Then if this system gives the exact solution of the problem, the 3 x 3 matrix is the inverse of the A matrix. The solution from this system is x, = 0.54687, x2 = 0.75000, x3 = 1.03125. 3. The
Storage capacity
COMPARISON
AND CONCLUSION
advantage of the LU decomposition method or the Gauss and Gauss-Jordan methods is well established [I].
REFERENCES
1. W. H. Press, B. P. Flannery, S. W. T. Vetterling, Numerical Recipes, Computing. Cambridge University (1986). 2. J. R. Westlake, A Handbook of
A. Teukolsky and The Art of Science
Press, New York
Numerical Matrix Inversion and Solution of Linear Equations. John Wiley,
New York (1968).