MATRICES AND MATRIX OPERATIONS: Unit 19

Dr. Wlodzislaw Kostecki

The Papua New Guinea University of Technology (PNGUT)

Department of Electrical and Communication Engineering

Lae, Morobe Province

Papua New Guinea

Copyright © 2000 by Wlodzislaw Kostecki

All rights reserved

-------------------------------------------------------------------

(19) Solution of systems of linear algebraic inhomogeneous equations using Cramer ’ s rule

OBJECTIVES :

• To define a system of linear algebraic inhomogeneous equations.

• To show how a system of such equations may be written in the matrix form.

• To introduce the concepts of the coefficient matrix , solution vector , and inhomogeneous vector .

• To introduce Cramer’ s rule for solving systems of linear algebraic inhomogeneous equations.

• To provide a numerical example for the application of Cramer’ s rule.

• To show how systems of linear algebraic inhomogeneous equations may be solved in a simple way using Maple ’ s functions linsolve or solve .

• To demonstrate that so-called ill-conditioned systems of linear inhomogeneous equations do not pose problems if tackled with Maple .

> restart : with(linalg, copyinto, det, genmatrix, linsolve, multiply, stackmatrix) :

Consider the following system of linear algebraic inhomogeneous equations in variables:

where the term inhomogeneous denotes that at least one of the numbers ... is not .

Since the number of variables equals the number of equations of this system, it is said to be a square system of linear equations.

This system of equations may be written in the matrix form. To this end, do the following.

• Define the square matrix containing coefficients as the coefficient matrix [ A ]:

> A := matrix(4, 4, [ a[11], a[12], a[13], a[14], a[21], a[22], a[23], a[24], a[31], a[32], a[33], a[34], a[41], a[42], a[43], a[44] ]) : A = matrix(A) ;

• Define the column matrix containing variables as the solution vector [ X ]:

> X := matrix(4, 1, [ x[1], x[2], x[3], x[4] ]) : X = matrix(X) ;

• Define the column matrix containing numbers as the inhomogeneous vector [ K ]:

> K := matrix(4, 1, [ k[1], k[2], k[3], k[4] ]) : K = matrix(K) ;

* * *

N.B. Both the "solution vector" and "inhomogeneous vector" are column structures and as such, cannot be named "vectors" in Maple – refer to Unit (5). However, because of a customary practice used in textbooks, both names are maintained in this Unit.

* * *

With the matrices defined as above, the system of the equations may now be written in the matrix form

[ A ] [ X ] = [ K ]

or,

> matrix(A) * matrix(X) = matrix(K) ;

Provided that the matrix [ A ] is non-singular , i.e. its determinant is not , it follows that the inverse to the matrix [ A ] exists. Multiplying each side of the matrix equation [ A ] [ X ] = [ K ] by Inv [ A ] yields

(Inv [ A ] ) [ A ] [ X ] = (Inv [ A ] ) [ K ]

Since (Inv [ A ] ) [ A ] = [ U ] [see Unit (14)] and [ U ] [ X ] = [ X ] [see Unit (9)], the above simplifies to yield the solution vector in the form

[ X ] = (Inv [ A ] ) [ K ]

or, in "like-in-a-book" form,

> matrix(X) = Inv(matrix(A)) * matrix(K) ;

* * *

N.B. When a solution vector [ X ] exists whose elements simultaneously satisfy all the equations in the system, the equations are said to be consistent . If no solution vector exists having this property, then the equations are said to be inconsistent .

* * *

A useful method of solution for small systems of such equations (maximum equations where computations are to be performed manually) is known as Cramer’ s rule. It states that any unknown is found as

= Det([ A , j <– K ]) / Det([ A ])

The determinant in the numerator is defined as the determinant of the modified coefficient matrix [ A ], in which the th column is replaced by the inhomogeneous vector [ K ], as indicated by the notation [ A , j <– K ] .

Exemplarily, the numerator for the unknown has the form

> i := 1 : j := 2 :

> `det(A, j <– K)` := copyinto(K, A, i, j) : Det(A, ` j <– K`) = Det(matrix(`det(A, j <– K)`)) ;

In the denominator, the determinant of the coefficient matrix [ A ] is called the characteristic determinant . For this particular case, it is given by the following expression:

> `det(A)` := det(A) : Det(A) = `det(A)` ;

Numerical example

Solve the following system of equations using Cramer’ s rule:

> Eq[1] := x[1] - 3*x[2] + x[3] - x[4] = -3 :

> Eq[2] := 2*x[1] - x[2] + 3*x[3] + x[4] = 2 :

> Eq[3] := x[1] + x[2] + 2*x[3] + 2*x[4] = 3 :

> Eq[4] := x[1] - 2*x[2] + 4*x[3] + x[4] = 0 :

> Eq[1] ; Eq[2] ; Eq[3] ; Eq[4] ;

Step 1 . Assign name and value to the number of equations and variables:

> No_Eq := 4 : No_x := 4 : 'No_Eq' = No_Eq ; 'No_x' = No_x ;

Step 2 . Construct the coefficient matrix [ A ]:

Method 1 . Using the function genmatrix and either of the two alternative methods of specifying the equations and variables:

• enter manually the names of all equations enclosed in brackets, followed by the names of all variables enclosed in brackets, viz.

> A := genmatrix([Eq[1], Eq[2], Eq[3], Eq[4]], [x[1], x[2], x[3], x[4]]) : A = matrix(A) ;

• create a bracketed sequence (list) of equations, followed by bracketed sequence (list) of variables, viz.

> `seq(Eq[i])` := [seq(Eq[i], i=1..No_Eq)] : `seq(x[i])` := [seq(x[i], i=1..No_x)] :

> A := genmatrix(`seq(Eq[i])`, `seq(x[i])`) : A = matrix(A) ;

Method 2 . Requiring several steps and the for -loop construct:

(a) Count the number of operands of the left-hand side of each equation:

> for j to No_Eq do if op(select(type, convert(lhs(Eq[j]), list), numeric)) <> NULL then ops(lhs_Eq[j]) := 1 else ops(lhs_Eq[j]) := nops(lhs(Eq[j])) fi : print(No_operands(lhs_Eq[j]) = ops(lhs_Eq[j])) : od :

(b) Create lists of operands of the left-hand side of each equation:

> for j to No_Eq do if ops(lhs_Eq[j]) = 1 then L[j] := [lhs(Eq[j])] else L[j] := convert(lhs(Eq[j]), list) fi : print(evaln(L[j]) = L[j]) : od :

(c) Construct row matrices containing numeric operands of the elements of each list:

> for i to No_x do n := 0 : for j to No_Eq do if has(L[i], x[j]) = false then a[i,j] := 0 else n := n+1 : a[i,j] := op(n, L[i])/x[j] fi : od : RM[i] := matrix(1, No_x, [[seq(a[i,m], m=1..No_x)]]) : print(RM[i] = matrix(RM[i])) : od :

(d) Create a sequence of the row matrices:

> `seq(RM[i])` := seq(matrix(RM[k]), k=1..No_Eq) : sequence(RM['i']) = `seq(RM[i])` ;

(e) Construct matrix [ A ] by stacking the row matrices:

> A := stackmatrix(`seq(RM[i])`) : A = matrix(A) ;

N.B. Method 2 is universal and works even if some coefficients in the simultaneous equations are . It also works even if one coefficient is non- in any equation (the case purely abstractive from mathematical point of view).

If none of the coefficients is , the following short variant of this method may be used:

> a := matrix(No_x, No_Eq) :

> for i to No_Eq do for j to No_x do a[i,j] := op(j, lhs(Eq[i]))/x[j] od : od :

> A := matrix(a) : A = matrix(A) ;

Step 3 . Create a copy of matrix [ A ] for the future use (in Step 7):

> C_A := copy(A) : C_A = matrix(C_A) ;

Step 4 . Define and input the solution vector [ X ]:

Method 1 . Using the for -loop construct:

> x := array(1..No_Eq) : for i to No_x do x[i] := x[i] od : X := convert(x, matrix) : X = matrix(X) ;

Method 2 . Using manual inputting of each vector element:

> X := matrix(4, 1, [x[1], x[2], x[3], x[4]]) : X = matrix(X) ;

Step 5 . Define and input the inhomogeneous vector [ K ]:

Method 1 . Using the for -loop construct:

> k := array(1..No_Eq) : for i to No_x do k[i] := rhs(Eq[i]) od : K := convert(k, matrix) : K = matrix(K) ;

Method 2 . Using manual inputting of each vector element:

> K := matrix(4, 1, [rhs(Eq[1]), rhs(Eq[2]), rhs(Eq[3]), rhs(Eq[4])]) : K = matrix(K) ;

The system of the equations may now be written in the matrix form, i.e.

> matrix(A) * matrix(X) = matrix(K) ;

Step 6 . Compute the characteristic determinant of the coefficient matrix [ A ]:

> `det(A)` := det(A) : Det(A) = `det(A)` ;

Step 7 . Find the unknowns ... :

(a) Finding

> `det(A, 1 <– K)` := copyinto(K, A, 1, 1) : Det(A, `1 <– K`) = Det(matrix(`det(A, 1 <– K)`)) ;

> `det(A, 1 <– K)` := det(`det(A, 1 <– K)`) : Det(A, `1 <– K`) = `det(A, 1 <– K)` ;

> x[1] := `det(A, 1 <– K)`/`det(A)` : 'x[1]' = x[1] ; x[1] := 'x[1]' :

(b) Finding

> A := matrix(C_A) :

> `det(A, 2 <– K)` := copyinto(K, A, 1, 2) : Det(A, `2 <– K`) = Det(matrix(`det(A, 2 <– K)`)) ;

> `det(A, 2 <– K)` := det(`det(A, 2 <– K)`) : Det(A, `2 <– K`) = `det(A, 2 <– K)` ;

> x[2] := `det(A, 2 <– K)`/`det(A)` : 'x[2]' = x[2] ; x[2] := 'x[2]' :

(c) Finding

> A := matrix(C_A) :

> `det(A, 3 <– K)` := copyinto(K, A, 1, 3) : Det(A, `3 <– K`) = Det(matrix(`det(A, 3 <– K)`)) ;

> `det(A, 3 <– K)` := det(`det(A, 3 <– K)`) : Det(A, `3 <– K`) = `det(A, 3 <– K)` ;

> x[3] := `det(A, 3 <– K)`/`det(A)` : 'x[3]' = x[3] ; x[3] := 'x[3]' :

(d) Finding

> A := matrix(C_A) :

> `det(A, 4 <– K)` := copyinto(K, A, 1, 4) : Det(A, `4 <– K`) = Det(matrix(`det(A, 4 <– K)`)) ;

> `det(A, 4 <– K)` := det(`det(A, 4 <– K)`) : Det(A, `4 <– K`) = `det(A, 4 <– K)` ;

> x[4] := `det(A, 4 <– K)`/`det(A)` : 'x[4]' = x[4] ; x[4] := 'x[4]' :

* * *

N.B. The solution vector [ X ] may be obtained in a simpler way using the linsolve function, viz.

> A := genmatrix(`seq(Eq[i])`, `seq(x[i])`) : X := linsolve(A, K) : X = matrix(X) ;

Extracting the numerical value of each unknown from the solution vector [ X ] gives

> for i to No_Eq do x[i] := X[i,1] : print(evaln(x[i]) = x[i]) : od :

* * *

Thus, the solved system of the equations may be written in the matrix form as follows:

> matrix(A) * matrix(X) = matrix(K) ;

Verification of the solution may be performed by evaluating the product of matrices [ A ] and [ X ], viz.

> `AX` := multiply(A, X) : A * X = matrix(`AX`) ;

which is equal to the inhomogeneous vector [ K ].

* * *

N.B. The simplest method of solving systems of such equations in Maple is by using the solve (or fsolve ) function. For instructive purposes, this is done hereunder for the above system of equations.

> for i to No_x do x[i] := evaln(x[i]) : od :

> solution := solve({Eq[1], Eq[2], Eq[3], Eq[4]}, {x[1], x[2], x[3], x[4]}) : solution ;

or

> i := 'i' : `seq(Eq[i])` := {seq(Eq[i], i=1..No_Eq)} : `seq(x[i])` := {seq(x[i], i=1..No_x)} :

> solution := solve(`seq(Eq[i])`, `seq(x[i])`) : solution ;

Extracting each solution from the unordered solution set gives

> assign(solution) : for i to No_x do x[i] := x[i] : print(evaln(x[i]) = x[i]) : od :

> for i to No_x do x[i] := evaln(x[i]) : od :

* * *

N.B. The solution of a system of simultaneous equations may, on occassion, become difficult due to round-off errors if a calculator is used for computations based on Cramer’ s rule. In general, this may take place in solving systems for which the characteristic determinant of the coefficient matrix is close to . Such systems are said to be ill-conditioned . A simple example is provided by the system

> Eq[1] := x[1] + x[2] = 4 :

> Eq[2] := x[1] + 1.0001*x[2] = 1 :

> Eq[1] ; Eq[2] ;

which, however, does not pose a problem when tackled using Maple and any of the methods presented earlier. This is demonstrated hereunder.

Method 1 . Using Cramer’ s rule:

> A := genmatrix([Eq[1], Eq[2]], [x[1], x[2]]) : A = matrix(A) ; C_A := copy(A) :

> X := matrix(2, 1, [x[1], x[2]]) : X = matrix(X) ;

> K := matrix(2, 1, [rhs(Eq[1]), rhs(Eq[2])]) : K = matrix(K) ;

> `det(A)` := det(A) : Det(A) = `det(A)` ;

(a) Finding

> `det(A, 1 <– K)` := copyinto(K, A, 1, 1) : Det(A, `1 <– K`) = Det(matrix(`det(A, 1 <– K)`)) ;

> `det(A, 1 <– K)` := det(`det(A, 1 <– K)`) : Det(A, `1 <– K`) = `det(A, 1 <– K)` ;

> x[1] := `det(A, 1 <– K)`/`det(A)` : 'x[1]' = x[1] ; x[1] := 'x[1]' :

(b) Finding

> A := matrix(C_A) :

> `det(A, 2 <– K)` := copyinto(K, A, 1, 2) : Det(A, `2 <– K`) = Det(matrix(`det(A, 2 <– K)`)) ;

> `det(A, 2 <– K)` := det(`det(A, 2 <– K)`) : Det(A, `2 <– K`) = `det(A, 2 <– K)` ;

> x[2] := `det(A, 2 <– K)`/`det(A)` : 'x[2]' = x[2] ; x[2] := 'x[2]' :

Method 2 . Using the linsolve function:

> A := genmatrix([Eq[1], Eq[2]], [x[1], x[2]]) : A = matrix(A) ;

> X := linsolve(A, K) : X = matrix(X) ;

Extracting the numerical value of either unknown from the solution vector [ X ] gives

> for i to 2 do x[i] := X[i,1] : print(evaln(x[i]) = x[i]) : od :

Method 3 . Using the solve function:

> for i to 2 do x[i] := evaln(x[i]) : od :

> solution := solve({Eq[1], Eq[2]}, {x[1], x[2]}) : solution ;

Extracting either solution from the unordered solution set gives

> assign(solution) : for i to 2 do x[i] := x[i] : print(evaln(x[i]) = x[i]) : od :

> for i to 2 do x[i] := evaln(x[i]) : od :

* * *

Yet another example of an ill-conditioned system are the simultaneous equations

> Eq[1] := 0.00001*x[1] + x[2] + 0.00001*x[3] = 0.00002 :

> Eq[2] := x[1] + 2*x[2] + x[3] = 1 :

> Eq[3] := 0.00001*x[1] + x[2] - 0.00001*x[3] = 0.00001 :

> Eq[1] ; Eq[2] ; Eq[3] ;

The coefficient matrix [ A ] is:

> A := genmatrix([Eq[1], Eq[2], Eq[3]], [x[1], x[2], x[3]]) : A = matrix(A) ;

The characteristic determinant of the coefficient matrix [ A ] has the value:

> `det(A)` := det(A) : Det(A) = `det(A)` ;

The inhomogeneous vector [ K ] is:

> K := matrix(3, 1, [rhs(Eq[1]), rhs(Eq[2]), rhs(Eq[3])]) : K = matrix(K) ;

With the default environment variable Digits ( i.e. ) , the solution vector obtained using the linsolve function is:

> X := linsolve(A, K) : X = matrix(X) ;

Extracting the numerical value of each unknown from the solution vector [ X ] gives

> for i to 3 do x[i] := X[i,1] : print(evaln(x[i]) = x[i]) : od :

Conclusion : The notion of ill-conditioned systems of linear inhomogeneous equations does not apply at all if Maple is used for solving such systems.

* * *

Proceed to Unit (20) for " Trace of a matrix ".

-------------------------------------------------------------------