UCES
Course: Methods and Analysis in UCES
Applied Area: inventory management, electrical circuits
Model: triangular algebraic system
Method: ij and ji upper triangular solve
Assessment: nonsingular triangular matrices
Prerequisites:
Objectives:
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 7-28-94
Copyright ©1994, R. E. White. All rights reserved.
The next four lectures will be concerned with solving systems of algebraic equations that have the form
where A is a given matrix, d is a given column vector and x is a column vector to be found. In this lecture we will focus on the special case where A is square and triangular. We will consider applications to inventory management, electrical circuits and later to the steady state heat conduction in a wire.
There several classes of problems which can be easily illustrated by the following examples.
Example 1. The algebraic system may not have a solution. Consider
If d = [1,2]T, then there are an infinite number of solutions
given by points on the line l1 in the figure below. If
then there will be exactly one solution given by the intersection of lines l1 and l3.
 |
| Figure: Solutions, infinite or none or one |
Example 2. This example illustrates a system with three equations with either no solutions or a set of solutions that is a straight line in 3D space.
If d2 
3, then the
second row or equation implies
3 1.
The vector form of the solution is
This is a straight line in 3D space containing the point [0,0,1]T and going in the direction [-1,1,0]T.
The easiest algebraic systems to solve are either ones where the matrix is a diagonal or a triangular matrix.
Example 3. Consider the case where A is a diagonal matrix.
Example 4. Consider the case where A is a lower triangular matrix.
The first row or equation gives x1 = 1. Use this in the second row or equation to get
Put these two into the third row or equation to get
This is known as a forward sweep. If the matrix were an upper triangular matrix, then one would have to start with the bottom row or equation, solve for the last component and substitute into the next to the last row or equation, and so forth. This is called a backward sweep.
In this section we present two applications. The first is to management of inventories that are required to produce a number of products. Here we consider a small number of products, but the size of such problems can be very large. The second is an application is the electrical circuits where we also restrict the size of the problem.
Application to Inventories. A small business produces three items:
If they have orders for x1 = 10 wood boxes, x2 = 4 work benches and x3 = 2 metal cabinets, then the
can be computed as follows:
Or, in matrix form
The above was resolved by a simple matrix multiplication. The related
inverse problem requires a solution of an upper triangular algebraic system.
Suppose the same business has a current inventory of
This is equivalent to the last rows or bottom equations listed first
Thus, with the present inventory, 16 wood boxes, 5 work benches and 5 metal cabinets can be built. This is an example of a backward sweep which is used to compute the solution when the matrix is upper triangular.
Application to Circuits. Consider a circuit with only two loops, 12 volt power source and three devices that have resistances of 2 ohms, 6 ohms and 1 ohm, see the figure below. The problem is to determine if the current in each device will meet the required amounts to operate the devices.
 |
| Figure: Circuit with Two Loops |
In each loop the voltage drops must sum to zero. Also Ohm's law requires voltage to be equal to current times resistance.
| Loop 1: i12 + i26 = 12, | (1) |
| Loop 2: +i26 + i31 = 0. | (2) |
| Node 1: i1 - i2 + i3 = 0, | (3) |
| Node 2: -i1 + i2 - i3 = 0. | (4) |
Equation (4) is equivalent to equation (3), and so, we delete it. The first three equations can be written in matrix form as
Unfortunately, the matrix is not in either lower or upper triangular form. However, by manipulating the equations (1)-(3) we can reformulate the problem so that the new matrix will have upper triangular form.
First, multiply equation (1) by 1/2 and subtract it form equation (3) to get
| 0i1 -4i2 + 1i3 = -6. | (3') |
| 0i1 + 0i2 + (5/3)i3 = -6. | (3'') |
Now, a backward sweep easily generates the desired currents
The systems of equations (1), (2),(3) and (1),(2),(3'') have the same solution because the equation manipulation to obtain (3'') is reversible. The process of putting the problem in upper triangular for is called Gaussian elimination, and we shall study this in more detail in the next two lectures.
The general model will be an algebraic system of n equations and n unknowns. We will assume the matrix has upper triangular form
i,j n.
The row numbers of the matrix are associated with i, and the column numbers are given by j. The component form of Ax = d when A is upper triangular is
aijxj = di.
The last equation is annxn = dn, and hence,
The (n-1)th equation is
and hence, we can solve for
This can be repeated, provided each aii is nonzero, until all xi have been computed. We have just proved the first part of the following theorem.
Upper Triangular Existence Theorem. Consider Ax = d where A is upper triangular (aij = 0 for i > j) and an n by n matrix. If each aii is not zero, then Ax = d has a solution. Moreover, this solution is unique.
In order to execute this on a computer, there must be two loops: one for the equation (i loop) and one for the summation (j loop). As in the case of matrix-vector products there are two versions: the ij version with the i loop on the outside, and the ji version with the j loop on the outside.
The ij version is a reflection of the backward sweep as previously presented. Note the inner loop retrieves data from the array by jumping from one column to the next. In Fortran this is in stride n and is not good for effective use of vector pipelines.
Upper Triangular Solve Algorithm: ij Version.
x(n) = d(n)/a(n,n)
for i = n-1,1
sum = d(i)
for j = i+1,n
sum = sum - a(i,j)*x(j)
endloop
x(i) = sum/a(i,i)
endloop.
In order to understand the ji version, consider the algebraic system from the inventory problem.
Recall that matrix-vector products can also be viewed as linear combinations of the columns of the matrix. For this example, this is
First, solve for x3 = 20/4 = 5.
Second, subtract the last column times x3 from both sides to reduce the dimension of the problem
Third, solve for x2 = 5/1.
Fourth, subtract the second column times x2 from both sides
Fifth, solve for x1 = 65/4.
By subtracting multiples of the columns of A, one changes the order of the loops. And, one retrieves components of A moving down the columns of A.
Upper Triangular Solve: ji version.
x(n) = d(n)/a(n,n)
for j = n,2
for i = 1,j-1
d(i) = d(i) - a(i,j)*x(j)
endloop
x(j-1) = d(j-1)/a(j-1,j-1)
endloop.
We illustrate some of the above computations in Maple. First, we have made use of Maple procedure call linsolve which solves algebraic systems. Next, we implement the ij version of the upper triangular solve algorithm.
Input Data:
> with(linalg):
> A:=matrix(3,3,[4,6,1, 0,1,1, 0,0,4]);
> d:=vector([100,10,20]);
Call to linsolve:
> x:=linsolve(A,d);
Output Data:
Maple Code for ij Version of Upper Triangular Solve.
Input Data:
> n:=3:
> x:=vector(n,0):
Algorithm is Defined:
> x[n] := d[n]/A[n,n]:
> for i from n-1 to 1 do
> sum:=d[i]:
> for j from i+1 to n by -1 do
sum:= sum - A[i,j]*x[j]:
> od;
> x[i]:=sum/A[i,i]:
> od;
Output Data:
> eval(x);
The Matlab version of the linear solve problems is a little less complicated. Also, we give the ji version of the upper triangular solve algorithm.
Matlab Code for Linear Solve.
>>A = [4 6 1;0 1 1;0 0 4]
A =
4 6 1
0 1 1
0 0 4
>>d = [100 10 20]'
d =
100
10
20
>>x = inv(A)*d
x =
16.2500
5.0000
5.0000
>>x = A\d
x =
16.2500
5.0000
5.0000
The last method of solving for x is the most efficient because it does not require the computation of the inverse matrix of A (more on this later). Another method is to use the ji upper triangular solve algorithm, jiutrisol.m.
Matlab Code for ji Upper Triangular Solve.
A = [4 6 1;0 1 1;0 0 4]
d = [100 10 20]'
n = 3
x(n) = d(n)/A(n,n)
for j = n:-1:2
for i = 1:j-1
d(i) = d(i) - A(i,j)*x(j);
end
x(j-1) = d(j-1)/A(j-1,j-1);
end
Output:
x =
16.2500
5.0000
5.0000
Both applications to inventory and circuits are very small in the number of unknowns. Moreover, it is likely that these applications will have multiple solutions.
One problem with the upper triangular solve algorithm may occur if aiiare very small. In this case the floating point approximation may induce significant error. An other instance is two equations which are nearly the same. For example, in 2D, suppose the lines associated with the two equations are almost parallel. Then small changes in the slopes, given by either floating point or empirical data approximations, will induce big changes in the location of the intersection, that is, the solution.
The existence theorem had a second result that any solution is unique. In order to prove this, let x and y be two solutions
Subtract these two and use the basic properties of matrix products.
Now apply the upper triangular solve algorithm with d replaced by 0 and x
replaced by
This implies x = y which is what we wanted to show.
The results of this lecture also hold for lower triangular matrices.
and additional construction material
The new matrix and d vector are
.
What are the meanings of the new entries? Solve for the number of items that can be produced with the present inventory of supplies d.
