Method of change of variable decomposition to solve and invert a system of linear equations

Method of change of variable decomposition to solve and invert a system of linear equations

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...

320KB Sizes 0 Downloads 28 Views

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).