% FILE anton1.tex

% Note re I for Identity matrix versus sqrt(-1) is at end

% Chapter 1. 
\section{Systems of Linear Equations and Matrices}

\subsection{Introduction to Linear Systems and Matrices}
The central problem of linear algebra is the solution of linear equations.
You should already have some experience of these, and in particular you should
know 
(see [AR] \S1.1)
that any system of linear equations can have 0, 1 or an infinite number of
solutions.
A common illustration of this involves the graphical presentation
of two linear equations in the two unknowns $x$ and $y$.
Each equation can be thought of as defining a line in the
$(x,y)$-plane.
When the 2 lines are\\
 distinct and parallel, the equations have 0 solutions;\\
 non-parallel, the equations have 1 solution;\\
 parallel and overlapping, the equations have infinitely many
solutions.

Maple has a command {\tt linsolve} 
\index{linsolve}
for solving linear systems.
Realistic engineering applications occur, for example,
in connection with engineering statics ([AR] \S11.3),
and in connection with electrical networks ([AR] \S11.19).
Maple Worksheets are available %on the disk supplement 
for both of these.

\smallskip
We begin by illustrating {\tt linsolve}
for the system given in [AR] \S1.1 Example 3:

% does one need with(linalg,linsolve):     ??????
\begin{maple}
# AR6 and AR7  Section 1.1 Example 3 p5-
# 
# Consider the linear system of equations
#   x +   y + 2*z = 9 ,
# 2*x + 4*y - 3*z = 1 ,
# 3*x + 6*y - 5*z = 0 .
#
# If the preceding system is written A x = b, the matrix A
# can be defined to Maple
with(linalg,linsolve):
A  :=array([[1,1, 2],
            [2,4,-3],
            [3,6,-5]]);
b  :=array([9,1,0]);
# The solution is obtained
linsolve(A,b);
                                  [ 1, 2, 3 ]
\end{maple}

The theory for solving linear systems, i.e. the remainder
of Chapter 1 is worth knowing.
However, the Maple of the rest of this chapter is
presented merely to indicate
the algorithms used in solving linear systems,
and can be skipped while you are content to use
{\tt linsolve} as a black box.
The basic idea is that of back-substitution from a system
in row-echelon form. In this approach, the key Maple commands are
{\tt augment}\index{augment}
for forming the augmented matrix of a system of equations, and 
{\tt backsub}\index{backsub}
for solving a system in row-echelon form:

\begin{maple}
# AR6 and AR7 Section 1.1 Example 3, leads to the following 
# system in row-echelon form.
> with(linalg,matrix,vector,augment,backsub,col):
# The coefficient and augmented matrices 
# of the following row-echelon system
#          x + y +  2  z =     9  ,
#              y - 7/2 z =  -17/2 ,
#                - 1/2 z =   -3/2 ,
#                    5 z = 15
# are given by
Coeff:= matrix([ [1, 1, 2],
                 [0, 1, -7/2],
                 [0, 0, -1/2] ]):
\end{maple}
\goodbreak
\begin{maple}
linSystem:= augment( Coeff, vector([9,-17/2,-3/2]) );
                                     [ 1  1   2    9   ]
                                     [                 ]
                        linSystem := [ 0  1 -7/2 -17/2 ]
                                     [                 ]
                                     [ 0  0 -1/2  -3/2 ]

# and the answer is given by
> backsub(");
                                  [ 1, 2, 3 ]
\end{maple}

We next move on to how
row operations enable us to transform the original system
$Ax=b$ into the easy-to-solve row-echelon form.
Maple has higher-level functions to do this, but we
begin by treating the underlying lower-level row operations.

\subsection{Gaussian Elimination ([J] \S1.2)}

In fact, for floating point numeric computation on ordinary serial
computers, Gaussian elimination is the best method for large systems.
Gaussian elimination involves the systematic use of
row-operations.
Maple has the row-operations defined in the functions
{\tt addrow}\index{addrow}, 
{\tt mulrow}\index{mulrow} and 
{\tt swaprow}\index{swaprow}.
We can use these to treat the same Example as in \cite{AR}
as follows.

\begin{maple}
# This concerns the "Gaussian elimination" method of solving linear systems. 
> with(linalg,matrix,vector,augment,addrow):
# Consider the linear system of equations as above
# This system can be written using an augmented matrix A_b = [A | b] where
> A:=matrix([[1,1,2],[2,4,-3],[3,6,-5]]):
> b:= vector( [9,1,0] ):                                 A_b:= augment(A,b);
                                  [ 1  1   2  9 ]
                                  [             ]
                           A_b := [ 2  4  -3  1 ]
                                  [             ]
                                  [ 3  6  -5  0 ]

# The example begins by reducing A_b to row-echelon form with row operations.
# It should be emphasised that these operations are likely to be more useful
# in learning some matrix manipulative processes than they are in actually
# solving equations, as the latter is now built into Maple's  linsolve.
# 
# The command to add -A[2,1]/A[1,1] times the first row to the second row is
> A1:=addrow(A_b,1,2,-A[2,1]/A[1,1]);
                                 [ 1  1   2   9  ]
                                 [               ]
                           A1 := [ 0  2  -7  -17 ]
                                 [               ]
                                 [ 3  6  -5   0  ]

# and similarly to add  -A[3,1]/A[1,1] times the first row to the third row:
\end{maple}
\goodbreak
\begin{maple}
> A2:=addrow(A1,1,3,-A[3,1]/A[1,1]);
                                [ 1  1   2    9  ]
                                [                ]
                          A2 := [ 0  2   -7  -17 ]
                                [                ]
                                [ 0  3  -11  -27 ]

> Ech:=addrow(A2,2,3,-A2[3,2]/A2[2,2]);
                                [ 1  1    2     9  ]
                                [                  ]
                         Ech := [ 0  2   -7    -17 ]
                                [                  ]
                                [ 0  0  -1/2  -3/2 ]

# If you have had some experience with Maple procedures,
# you could put the above row operations into a procedure.
# (If you are just beginning your Maple you should ignore this next part.)
#
> noswaps:=proc(A)
   local A1,A2,Ech;
   A1:=addrow(  A,1,2,-A[2,1]/A[1,1] ):
   A2:=addrow( A1,1,3,-A1[3,1]/A[1,1]):
   Ech:= addrow( A2, 2,3,-A2[3,2]/A2[2,2]); 
  end:
#
> with(linalg,equal,backsub,pivot,col);
> equal( noswaps(A_b),Ech );
                                      true

# showing that the result of  noswaps  is  A3.  The answer to the system is
> backsub(Ech);
                                  [ 1, 2, 3 ]
\end{maple}

{\tt gausselim(Aug)}\index{gausselim}
gives the same numerical result as {\tt Ech} above.
\smallskip
\noindent A case where row swaps are required is given in 
the following example.

\begin{maple}
# AR6 and AR7 Section 1.2 Example 3, p13
> with(linalg):
# Consider the linear system of equations
#   x1 + 3*x2 - 2*x3        + 2*x5        = 0 ,
# 2*x1 + 6*x2 - 5*x3 - 2*x4 + 4*x5 - 3*x6 = -1 ,
#               5*x3 +10*x4        +15*x6 = 5 ,
# 2*x1 + 6*x2        + 8*x4 + 4*x5 +18*x6 = 6 .
# The associated augmented matrix is
> A_b:=matrix([[1,3,-2, 0,2, 0, 0],
               [2,6,-5,-2,4,-3,-1],
               [0,0, 5,10,0,15, 5],
               [2,6, 0, 8,4,18, 6]]):

> addrow(A_b,1,2,-A_b[2,1]/A_b[1,1]):
\end{maple}
\goodbreak
\begin{maple}
> A2:=addrow(",1,4,-A_b[4,1]/A_b[1,1]);
                            [ 1  3  -2   0  2   0   0 ]
                            [                         ]
                            [ 0  0  -1  -2  0  -3  -1 ]
                      A2 := [                         ]
                            [ 0  0   5  10  0  15   5 ]
                            [                         ]
                            [ 0  0   4   8  0  18   6 ]

> addrow(A2,2,3,-A2[3,3]/A2[2,3]):
A3:=addrow(",2,4,-A2[4,3]/A2[2,3]);
                            [ 1  3  -2   0  2   0   0 ]
                            [                         ]
                            [ 0  0  -1  -2  0  -3  -1 ]
                      A3 := [                         ]
                            [ 0  0   0   0  0   0   0 ]
                            [                         ]
                            [ 0  0   0   0  0   6   2 ]

# In this case, the pivot is in the fourth row, and so proceed by
# interchanging the third and fourth rows:
# (All the zero rows of an Echelon form must be at the bottom.)
> with(linalg,swaprow,gausselim):
> Ech:=swaprow(A3,3,4);
                            [ 1  3  -2   0  2   0   0 ]
                            [                         ]
                            [ 0  0  -1  -2  0  -3  -1 ]
                     Ech := [                         ]
                            [ 0  0   0   0  0   6   2 ]
                            [                         ]
                            [ 0  0   0   0  0   0   0 ]


# The  gausselim  function in the linalg package also gives this matrix
> equal( gausselim(A_b),Ech );
                                      true

# There are numerous items of this general kind in AR.
# A treatment of issues involved in floating point approximate
# computations in this is given in AR Chapter 9.
\end{maple}

(Another item in Maple V which can be mentioned here is
{\tt pivot}\index{pivot}.)

A treatment of LU factorization is given in 
\cite{AR} Chapter \ChNM.
% 1987 edition A&R S8.2.



\par\noindent{\it An optional remark.} There are some theoretical
(but not practical numerical) uses in continuing the elimination process
until the matrix is in a still simpler {\it row-reduced echelon form}
(defined in \cite{AR} at beginning of \S1.2).
The difference is that all pivots are normalized to $+1$, by dividing each row
by a constant, and zeros are produced not only below but also above each pivot.
The Maple command for this is\quad 
{\tt rref}\index{rref}
or 
{\tt gaussjord}\index{gaussjord}.
\bigskip

\subsection*{\ARsixlplll\ Homogeneous Systems of Linear Equations}
\ifodd\ARed\setcounter{subsection}{2}\else\setcounter{subsection}{3}\fi

\setcounter{theorem}{0}
\begin{theorem} A homogeneous system of linear equations with more
unknowns than equations has infinitely many solutions.
\end{theorem}

The Maple functions {\tt kernel} or {\tt nullspace} are
useful for solving homogeneous systems.
\bigskip

\subsection{Matrices and Matrix Operations ([J] \S1.3,1.4)}
These have been considered in the opening section on the  {\tt linalg} package.

In this subsection we mostly present some exercises.
%\medskip
\goodbreak

\headq{EXERCISES.}
\begin{enumerate}
\item Let

  \verb$
  R:= matrix(2,2);                        S:= matrix(2,2);$
  \begin{enumerate}
    \item Use Maple to calculate the products $RS$ and $SR$. \\
      Check that you could find them without Maple too. 
    \item Use the function {\tt equal} to discover if these products are the
      same.
  \end{enumerate}
  (The systematic study of matrices, their noncommutative multiplication, 
  etc., started to interest people relatively recently, only last century.)

\item
%EXERCISE. 
Let
\begin{verbatim}
> R:=array([[-1,0],
            [ 0,1]]); 
> S:=array([[ 0,1],
            [-1,0]]); 
\end{verbatim}
(Geometrically $R$ is associated with a reflection and $S$ with a rotation
through a right angle. Think about this.) Use Maple to calculate the
products $RS$ and $SR$. Check that you could find them without Maple too.
What is $(RS - SR)$? 
%(The systematic study of matrices, their
%noncommutative multiplication, etc., started to interest people
%relatively recently, only last century.)




\item Use the {\tt linalg} multiplication to CHECK the answers in 
[AR6]  \S1.4 Example 6 (p29),
[AR7]  \S1.3 Example 5 (p29).
\end{enumerate}

\noindent The following example illustrates the different ways of viewing matrix
multiplication:

\begin{maple}
# AR6 Section 1.4 p30; AR7 Section 1.3 p30 Example 6
> with(linalg,matrix,vector,col):
> A:= matrix(2,3, (i,j) -> a.i.j):
> X:= vector(3, j -> x.j): AX:= evalm( A&*X ):  B:= vector(2, j -> b.j):
> print( AX[1] = B[1] );                  print( AX[2] = B[2] );
                        a11 x1 + a12 x2 + a13 x3 = b1
                        a21 x1 + a22 x2 + a23 x3 = b2
# Thus the system may be represented by the matrix equation  A*X = B.
# Now make the substitution
#
# The relationship between matrix multiplication and the columns of a matrix
> A1:= col(A,1):              A2:= col(A,2):              A3:= col(A,3):
> evalm( x1*A1 + x2*A2 + x3*A3 - AX );
                                    [ 0, 0 ]

# This shows that  AX  is the linear combination  x1*A1 + x2*A2 + x3*A3.
\end{maple} 


\subsection{Rules of Matrix Arithmetic; Inverses}

\headq{EXERCISE.}\quad Use the {\tt linalg} {\tt inverse} command to
\begin{enumerate}
\item[]
  \begin{enumerate}
  \item check the matrix inverses given in 
[AR6] \S1.5 Example 8 (p. 41),
[AR7] \S1.4 Example 7 (p. 44);
  \item verify, for the matrices $R$ and $S$ of the previous exercises, that
    \quad $(RS)^{-1} = S^{-1}R^{-1}$\quad but\quad $(RS)^{-1} \ne R^{-1}S^{-1}$.
  \end{enumerate}
\end{enumerate}
\medskip

\goodbreak
\subsection{Elementary Matrices and a Method for Finding $A^{-1}$}

\headq{Definition.} An $n$ by $n$ matrix is called an
{\it elementary matrix} if it can be obtained from the
$n$ by $n$ identity matrix by performing a single elementary row operation.

Thus the elementary matrices can be obtained in Maple 
by using the row operations
{\tt addrow}\index{addrow}, 
{\tt mulrow}\index{mulrow}, 
{\tt swaprow}\index{swaprow}
on the identity matrix:
\smallskip

\begin{maple}
# 
> with(linalg):
> I3:= diag( 1,1,1):                                   M1:= mulrow(I3,1,c):
> P23:= swaprow(I3,2,3):                             E13:= addrow(I3,1,3,m):
> print( M1,P23,E13 );
                     [ c  0  0 ]  [ 1  0  0 ]  [ 1  0  0 ]
                     [         ]  [         ]  [         ]
                     [ 0  1  0 ], [ 0  0  1 ], [ 0  1  0 ]
                     [         ]  [         ]  [         ]
                     [ 0  0  1 ]  [ 0  1  0 ]  [ m  0  1 ]
\end{maple}
\smallskip

\headq{EXERCISE.} Perform the same row operations on the matrix
\begin{verbatim}
A:= matrix(3,4, (i,j) -> a.i.j);
\end{verbatim}
and show that the results are the same as $M1*A$, $P23*A$ and $E13*A$
respectively. 

\smallskip
\headq {Definition.} An {\it elementary permutation matrix} is a matrix
$P_{ij}$ corresponding to the Maple command {\tt swaprow}, and a
{\it permutation matrix} is a product of such matrices.
\smallskip

\headq{EXERCISE.} (a) Show that the product of two permutation matrices 
is a permutation matrix. \\
(b) Show that any permutation matrix is invertible
and its inverse is also a permutation matrix.
\medskip

Returning to the more general case, treating elementary matrices
without regard to their kind, we have the following.
\smallskip

\setcounter{theorem}{1}
\begin{theorem}
%\par\noindent{\bf Theorem (AR 1.6.2).} {\it
  Each elementary matrix is invertible, and its inverse is an elementary matrix
  of the same type.
\end{theorem}

\noindent{\bf EXERCISES.}
\begin{enumerate}
\item Verify the preceding theorem for the matrices $M1$, $P23$ and $E13$.
\item Verify, for the matrix $P23$ above, that\quad $P23^2 = I_3$\quad and so
  \quad $P23^{-1} = P23$. 
\end{enumerate}

\goodbreak

%\setcounter{theorem}{2}
\begin{theorem}
%\par\noindent{\bf Theorem (AR 1.6.3).} {\it
  If $A$ is an $n\times n$ matrix, then the following 
  statements are equivalent:
  \begin{enumerate}
  \item[]
    \begin{enumerate}
    \item $A$ is invertible;
%    \item $AX=B$ has a unique solution;
    \item $AX=\b 0$ has the unique solution $X=\b 0$;
    \item $A$ is row-equivalent to the identity matrix $I_n$;
    \item $A$ is a product of elementary matrices.
    \end{enumerate}
  \end{enumerate}
\end{theorem}
\bigskip

The fact that any matrix can be reduced to Echelon form by pre-multiplying it
by a sequence of (invertible) elementary matrices suggests a method for
inverting a matrix by hand illustrated by 
(although Maple has a built-in function to do this).
\smallskip
\goodbreak

\begin{maple}
# Example 4  - AR6 Section 1.6, p51; AR7 Section 1.5, p55
> with(linalg,matrix,diag,augment,gausselim):
A:= matrix([ [ 1, 2, 3],
               [ 2, 5, 3],
               [ 1, 0, 8] ]):
# Use the function  augment  to construct the matrix  [A : I]
#
> Id:= diag( 1,1,1):                       A_I:= augment(A,Id):
#
# Reduce this matrix to Echelon form
AI1:= gausselim(A_I);

                               [ 1  2   3   1  0  0 ]
                               [                    ]
                        AI1 := [ 0  1  -3  -2  1  0 ]
                               [                    ]
                               [ 0  0  -1  -5  2  1 ]

# Do a formal back-substitution on this using Maple's row operations 
> with(linalg,addrow,mulrow,submatrix,equal):
addrow(AI1,3,2,-3): AI2:=addrow(",3,1,3);
                             [ 1  2   0  -14   6   3 ]
                             [                       ]
                      AI2 := [ 0  1   0   13  -5  -3 ]
                             [                       ]
                             [ 0  0  -1   -5   2   1 ]

> AI3:=addrow(AI2,2,1,-2);
                             [ 1  0   0  -40  16   9 ]
                             [                       ]
                      AI3 := [ 0  1   0   13  -5  -3 ]
                             [                       ]
                             [ 0  0  -1   -5   2   1 ]

# Now divide each row by the diagonal element.
# The inverse of the matrix A is now the right hand half of this matrix:
> AI4:=mulrow(AI3,3,-1):  
  AInverse:=submatrix(AI4,1..3,4..6);
                                     [ -40  16   9 ]
                                     [             ]
                         AInverse := [  13  -5  -3 ]
                                     [             ]
                                     [  5   -2  -1 ]
> equal( evalm(AInverse&*A), Id );
                                      true
\end{maple}
\medskip

\headq{EXERCISE.} At each stage of Example 4 above, let\par
  \smallskip
  \verb$ A.i:= submatrix(AI.i,1..3,1..3);$\par
  \verb$ E.i:= submatrix(AI.i,1..3,4..6);$

  \smallskip
Verify that\verb$ evalm( E.i &* A ) = A.i$.

\vfill
\eject
%\goodbreak
\subsection{Further Results on Linear Systems}

[AR6] p56, [AR7] p60
treat the problem of solving multiple
systems $Ax=b_j$ with the same $A$ but many different
right hand sides.
By augmenting $A$ with the many right-hand sides and
then using Gaussion elimination the calculation is
arranged so that one does better than solving each
individual system independently.

This is related to the following fact.
In float numerical treatment of linear systems,
e.g. as in \cite{AR} Chapter \ChNM, it is useful to view
gaussian elimination as a matrix factorization,
a topic we mention in the next subsection.

\smallskip
 

\subsection*{Some special sorts of matrices}

\headq{Definition.} A square matrix in which all the entries off the
main diagonal are zero is called a {\it diagonal matrix}.
\par\noindent
A Maple function which can be used to construct diagonal matrices is
{\tt diag}.
\medskip

\headq{Definition.} A square matrix is said to be 
{\it lower triangular} if all
entries above and right of the main diagonal are zero. 
(See [AR6] \S2.2, [AR7] \S1.7.)
\par\noindent These matrices can be formed using Maple like:
\begin{maple}
with(linalg,matrix):  l:='l':
lower:= proc(i,j) if (i<j) then 0 else l.i.j fi end:
matrix(3,3,lower);
\end{maple}

%\par\noindent{\bf THEOREM.}
\theoremNN{\it
 The product of two $n$ by $n$ lower triangular matrices is 
lower triangular.
}
\theoremNNEnd
\headq{EXERCISE.}\quad 
Verify the preceding theorem for general 3 by 3 lower
triangular matrices.
\smallskip

\headq{Definition.} A lower triangular matrix is said to be {\it unit lower
triangular} if its diagonal entries are all one.

%\par\noindent{\bf THEOREM.}
\theoremNN{\it
 The product of two $n$ by $n$ unit lower triangular matrices is 
unit lower triangular.
\par\noindent
The inverse of a unit lower triangular matrix is unit lower
triangular.
}
\theoremNNEnd
\headq{EXERCISE.}\quad 
Verify the preceding theorem for general 3 by 3 unit lower
triangular matrices.

\smallskip
There are corresponding definitions and results for {\it upper
triangular} matrices. See [AR7] Theorem 1.7.1.
Later we will see uses of triangular matrices in connection
with the $LU$-factorisation method for numerical solution of
linear systems. Here is a typical theorem, very closely
related to the row-operations methods discussed in Anton's
treatment of the preceding subsection.

%\setcounter{theorem}{62}
%\begin{theorem} 
%\par\noindent{\bf THEOREM.} 
\theoremNN{\it
Let $A$ be an invertible matrix.
  Then there is a permutation matrix $P$ such that
    $$ P A = L U ,$$
  where $L$ is unit lower triangular, and $U$ is upper triangular
  (both depending on $P$).
}
\theoremNNEnd
%\end{theorem}

\medskip

\headq{Definition.} A square matrix is said to be 
{\it symmetric} if $A=A^T$.
\par\noindent
The importance of symmetric matrices will be seen 
in later chapters.
For the present, it is worth noting that the symmetry of a matrix
may be used when storing it, since only the values on and below the diagonal
need be entered.  Maple takes advantage of this by allowing you to type,
for example\index{transpose}
\begin{maple}
> A:= array( symmetric,1..3,1..3 ):
> with(linalg,transpose,equal):
> equal( A,transpose(A) );
                                      true
\end{maple}
\medskip

\headq{Definition.} A square matrix is said to be 
{\it skew-symmetric} if $A=-A^T$.
\par\noindent
See [AR6] p46, [AR7] p74.
Skew-symmetric matrices, although interesting,
are less important in practice than are symmetric ones.
\medskip

A generalisation of `diagonal' matrices is the following.
\headq{Definition.} A  $n$ by $n$ matrix is said to be 
{\it banded} if $A=(a_{ij})$ satisfies
$$ \vert i-j\vert > k \qquad\Longrightarrow a_{ij}=0 ,$$
and then its {\it bandwidth} is defined to be $2k+1$.
\par\noindent
Thus a diagonal matrix has bandwidth 1.
If the bandwidth is  3 the matrix is said to be {\it tridiagonal}.
To avoid conflict with Maple's {\tt band} function, we use the word
{\em banded} where other linear algebra texts use {\em band}.  
Later we will see that
symmetric banded matrices can have a nice factorization.
We have a Worksheet using this to solve for
the stresses in a truss (bridge) with given loads.


\medskip
Your texts and lectures will stress that, in solving large sets of
linear equations in floats, ONE DOES NOT USE {\tt det} OR {\tt inverse},
but uses matrix factorization techniques.
See the appropriate parts of \cite{AR} Chapter \ChNM.

This is especially appropriate in the case of banded matrices such as above.  The
number of operations for factoring such an $n\times n$ banded matrix is $O(n)$,
and for solving the pair of triangular equations 
is also $O(n)$.
The inverse matrix is generally not a banded matrix. Forming $A^{-1}\b b$ is
at least $O(n^2)$ -- i.e. many, many more operations.

\medskip
Interestingly, for exact solution in rationals the situation is different
in that rationals themselves have different complexities, and this means
that multiplications can have different costs.  At this stage,
let's agree to leave it to the package designers and other specialists.
It is rare in practice to want to solve very large systems exactly. 

\vfill
\eject
\endinput
\vfill \rightline{AR1.tex}
\end{document}
\end

\medskip
% \subsection{Numerical and Programming Considerations}
Maple is an example of a Computer Algebra System, which gives exact solutions
to problems.
However, most problems cannot be solved exactly, and even if they could,
in practice one only wants approximations to answers or graphics.
For this reason, the advances in computing have led to
`numerical analysis' taking an increasingly important role in science. 
In practice numerical libraries are very important.


VARIOUS APPLICATIONS NOT IN CHAPTER 1

We now give two examples which are related to using Maple for real problems, 
rather than illustrating mathematical algorithms.

\bigskip
\goodbreak
\subsection*{An application of systems of linear equations: 
DC electrical networks}


\begin{verbatim}
###########################################################################
#  A&R S10.20 Electrical networks                                         #
###########################################################################
#           R1             R2
#     .---/\/\/\----A----/\/\/\---.
#     !             !             !
#     !             \             !
#     !             /             !
#     !             \             !
#     !             / R3          !
#     !             \             !
#     !             /             !
#     !      !      !      !      !
#     !-----!!------B-----!! -----!
#            !             !
#           V1            V2
#
# Let I.j be the current flowing through resistor R.j
Z0:= I1-I2-I3: # Kirchoff's current law with Z0 standing for zero ..
Z1:= R1*I1+R3*I3 -V1: # Kirchoff's voltage law, for loop 1
Z2:= R2*I2-R3*I3 -V2: # Kirchoff's voltage law, for loop 2
solve({Z0=0,Z1=0,Z2=0},{I1,I2,I3});
#            - R2 V1 + R1 V2            R3 V1 + V2 R3 + R2 V1   
#{I3 = - -----------------------, I1 = -----------------------, 
#         R2 R1 + R2 R3 + R1 R3         R2 R1 + R2 R3 + R1 R3   
#
#         R3 V1 + R1 V2 + V2 R3  
#   I2 = -----------------------}
#         R2 R1 + R2 R3 + R1 R3  
#

\end{verbatim}

with(linalg,linsolve):
 A:= array([[1,-1,-1],[R1,0,R3],[0,R2,-R3]]):
  b:= array([0,V1,V2]):
  currvec:= linsolve(A,b);
#currvec := array ( 1 .. 3, [
#
# R3 V1 + V2 R3 + R2 V1    R3 V1 + R1 V2 + V2 R3         - R2 V1 + R1 V2     
#----------------------, -----------------------, - ----------------------
#R2 R1 + R2 R3 + R1 R3    R2 R1 + R2 R3 + R1 R3      R2 R1 + R2 R3 + R1 R3  
#])

 for j from 1 to 3 do
    subs(R1=7,R2=11,R3=3,V1=30,V2=50,currvec[j]);
    od;
#               570               590               20 
#              -----             -----          - -----
#               131               131              131
# These are the numerical values given in A&R S10.20 Example 1.
#
###########################################################################
#
# A&R S10.20 Q5 Wheatstone bridge
#
#                   1
#                .  .  .
#             R1    .    R3
#            .      .      .
#    I0->  0        R5        3 ->I0
#            .      .      .
#             R2    .    R4
#                .  .  .
#                   2
#
#
# Current I0 flows round an outer loop in which is located a battery
# voltage V. The currents I.j, j=0..5 are unknown. The I.j
# for j=1..5 are the currents through the resistors R.j shown in Fig.
# They are in dirn from node i to node k with i<k.
Aw:= array([[1,-1,-1,0,0,0],
            [0,1,0,-1,0,-1],
            [0,0,1,0,-1,1],
            [0,R1,-R2,0,0,R5],
            [0,0,0,R3,-R4,-R5],
            [0,R1,0,R3,0,0]]):

#
bw:= array([0,0,0,0,0,V]):
currvecw:= linsolve(Aw,bw):
solve(currvecw[6]=0,R4);
# When R4 is
#                                     R2 R3 
#                                    -------
#                                       R1  
# the current I5 across the central arc is zero.
###########################################################################
# Corresponding code for the Wein bridge, used for AC circuits, is given in
# D. Harper, C. Wooff and Hodgkinson "A guide to CA systems" (Wiley 1991).
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Mention truss example, refering to disk supplement for
% details. This is an example of banded matrices occuring
% in a practical problem.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

