For example, the augmented matrix in the circuit problem is

Each row of the augmented matrix represents an equation in the algebraic system.

The next step was to add or subtract multiples of rows to get all zeros in the lower triangular part of the matrix. There are three basic row operations:

(i). interchange the order of two rows or equations,

(ii). multiply a row or equation by a nonzero constant and

(iii). add or subtract rows or equations.

We used a combination of (ii) and (iii). Each row operation is equivalent to a multiplication by an elementary matrix.

Example. Consider the above problem.

First, subtract 1/2 of row 1 from row 3 to get a zero in the (3,1) position:

Second, add 2/3 of row 2 to row 3 to get a zero in the (3,2) position:

It is customary to write this as

Each elementary row operation can be reversed, and this has the form of a matrix inverse of each elementary matrix:

Note that A = LU where L = E₁^-1E₂^-1 because by repeated use of the associative property

(E1-1E2-1)(EA)	= (E1-1E2-1)((E2E1)A)
	= ((E1-1E2-1)(E2E1))A
	= (E1-1(E2-1(E2E1)))A
	= (E1-1((E2-1E2)E1))A
	= (E1-1E1)A
	= A.

The product L = E₁E₂ is a lower triangular matrix and A = LU is called an LU factorization of A.

Definition. An n by n matrix, A, has an inverse n by n matrix, A^-1 , if and only if

A^-1A = AA^-1 = I = the n by n identity matrix.

Basic Properties. Let A be an n by n matrix which has an inverse.

1. A^-1 is unique.

2. x = A^-1d is a solution to Ax = d.

3. (AB)^-1 = B^-1A^-1 provided B also has an inverse.

4. A^-1 = [c₁,...,c_n] has column vectors that are solutions of Ac_j = e_j where
   e_j are unit column vectors with all zero components except the j^th which is equal to one.

We will later discuss these properties in more detail. Note, given an inverse matrix one can solve the associated linear system. Conversely, if one can solve the linear problems in property 4 via Gaussian elimination, then one can find the inverse matrix. Also one can use the elementary matrices to find the LU factorizations and the inverses of L and U; now, apply property 3 to find the inverse of A.

1.6.2.   Applied Area.

Steady state solutions have models which are also derived from Fourier's heat law. The difference now is that the change, with respect to time, in the heat content is zero. Also the temperature is a function of just space so that u_i ~ u(ih).

change in heat content = 0 ~ (heat from the source)

+ (heat from diffusion from the left side)

+ (heat from diffusion from the right side).

Figure: Heat Diffusion in a Thin Wire

The heat diffusing in the right face (when (u_i+1 - u_i)/h > 0) is

A dt K (u_i+1 - u_i)/h.

The heat diffusing out the left face (when (u_i - u_i-1)/h > 0) is

A dt K (u_i - u_i-1)/h.

Therefore, the heat from diffusion is

A dt K (u_i+1 - u_i)/h - A dt K (u_i - u_i-1)/h.

The heat from the source is

A h dt f.

By combining these we have the following approximation of the change in the heat content for the small volume Ah:

0 = A h dt f + A dt K (u_i+1 - u_i)/h - A dt K (u_i - u_i-1)/h.

Now, divide by Ah dt , let ß = K/h², and we have the following n-1 unknowns and n-1 equations.

Finite Difference Equations for Heat Transfer.

0 = f + ß (ui+1 + ui-1) - 2ßui	(1)
where
i = 1,...,n-1 and ß = K/h2,
u0 = un = 0 .	(2)

Equation (2) is the temperature at the left and right ends set equal to zero.

1.6.3.   Model.

The finite difference model may be written in matrix form where the matrix is a tridiagonal matrix. For example, if n = 4, then we are dividing the wire into four equal parts and there will be 3 unknowns with the end temperatures set equal to zero.

Tridiagonal Algebraic System: n = 4.

Suppose the length of the wire is 1 so that h = 1/4, and the thermal conductivity is .001. Then ß = .016 and if f_i = 1, then upon dividing all rows by ß and using the augmented matrix notation we have

Forward Sweep (put into upper triangular form):

Add 1/2(row 1) to (row 2),

Add 2/3(row 2) to (row 3),

Backward Sweep (solve the triangular system):

u₃ = (2)62.5 (3/4) = 93.75,

u₂ = ((3/2)62.5 + 93.75) (2/3) = 125 and

u₁ = (62.5 + 125)/2 = 93.75.

Tridiagonal matrices often factor into lower and upper triangular matrices with two diagonals. Later we will learn how to take advantage of this nice structure.

1.6.4.   Method.

The general Gaussian elimination method requires the augmented matrix, forward sweep to convert the problem to upper triangular form, and the backward sweep to solve this upper triangular system. The row operations needed to form the upper triangular system must be done in a systematic way:

(i). Start with column 1 and row 1 of the augmented matrix. Use an appropriate multiple of row 1 and subtract it from row i to get a zero in the (i,1) position in column 1 with i > 1.

(ii). Move to column 2 and row 2 of the new version of the augmented matrix. In the same way use row operations to get zero in each (i,2) position of column 2 with i > 2.

(iii). Repeat this until all the components of altered augmented matrix are zero.

This is depicted in the following figure where the (i,j) component is about to be set to zero.

Figure: Gaussian Elimination

Gaussian Elimination Algorithm.

        define the augmented matrix [A d]
        for j = 1,n-1
            for i = j+1,n  	(forward sweep)
                add multiple of (row j) to (row i) to get a zero in the (i,j) position
            endloop
        endloop
	for i = n,1             (backward sweep)
            solve for xi using row i
	endloop.

The above description is not very complete. In the forward sweep more details and special considerations with regard to roundoff errors are essential. The backward sweep is just the upper triangular solve step, and two versions of this were studied in the previous section.

1.6.5.   Implementation.

Maple has a number of intrinsic procedures which are useful for illustration of Gaussian elimination. These include linsolve, augment, gausselim, backsub, pivot, addrow and others.

First, we use linsolve which is a combination of the forward and backward sweeps.

Maple Code and Intrinsic Procedures.

   Input Data:
        >n:=4.:
        > h:=1/(n+1):
        > K:=.001:
        > beta:=K/(h*h):
        > with(linalg):
        > d:=vector([1/beta,1/beta,1/beta,1/beta]);

        

        > A:=matrix(4,4,[2,-1,0,0, -1,2,-1,0, 0,-1,2,-1, 0,0,-1,2]);

        
  
   Call to linsolve:
        > temp:=linsolve(A,d);

   Output Data:
        

Second, we use a combination of intrinsic procedures. The procedure augment forms the augmented matrix. The procedure gausselim does the forward elimination sweep which is the row operations to get zeros in the lower triangular part of the augmented matrix. And, the procedure backsub does the backward sweep to solve the upper triangular problem.

  Define Augmented Matrix:
        > Aug:=augment(A,d);
      
        

  Forward Sweep:
        > UD:=gausselim(Aug);
 
       

  Backward Sweep:
        > temp:=backsub(UD);

        .

Matlab Code for Linear Solves (heatgelm.m) with Variable n.

        n = 4
        h = 1./(n+1)
        K = .001;
        beta = K/(h*h);
        A= zeros(n,n);
        for i=1:n
            d(i) = 1/beta;
            A(i,i) = 2;
            if i<n 
               A(i,i+1) = -1;
            end;
            if i>1 
               A(i,i-1) = -1;
            end;
        end
        d = d'
        A
        temp = A\d

   Output for n = 4:
        temp =
           80.0000
          120.0000
          120.0000
           80.0000

   Output for n = 8:
        temp =
           49.3827
           86.4198
          111.1111
          123.4568
          123.4568
          111.1111
           86.4198
           49.3827

   Output for n = 16:
        temp =
           27.6817
           51.9031
           72.6644
           89.9654
          103.8062
          114.1869
          121.1073
          124.5675
          124.5675
          121.1073
          114.1869
          103.8062
           89.9654
           72.6644
           51.9031
           27.6817

1.6.6.   Assessment.

The above model for heat conduction depends upon the mesh size, h, but as the mesh size h goes to zero there will be little difference in the computed solutions. For example in the Matlab output, the i^th component of temp is the approximate temperature at ih where h = 1/(n+1). The approximate temperatures at the center of the wire are 120.0000 for n = 4, 123.4568 for n = 8 and 124.5675 for n = 16. An analysis of this will be given in a later section.

The four basic properties of inverse matrices needs some justification.

Proof that inverse are unique:

Let B and C be inverses of A so that AB = BA = I and AC = CA = I.
Subtract these matrix equations and use the distributive property

AB - AC = I - I = 0

A(B - C) = 0.

Use B is an inverse of A and use the associative property

B(A(B - C)) = B0 = 0

(BA)(B - C) = 0

B - C = 0.

Proof that A^-1d is a solution of Ax = d:

Let x = A^-1d and again use the associative property

A(A^-1d) = (AA^-1)d = I d = d.

Proofs of properties 3 and 4 are also a consequence of the associative property.

1.6.7.   Homework.

1. Find the LU factors for the matrix in the circuit problem.

2. Experiment with the mesh size in the heat conduction problem and verify that the computed solution converges as the mesh goes to zero.

3. Prove property 3.

4. Prove property 4.

5. Prove that the solution of Ax = d is unique if A^-1 exists.