Gauss was a project that I settled upon to test my programming
skills in a language to which I had just recently been introduced: LISP.
I also find it useful in linear algebra. Gauss cannot help with
everything that Gaussian elimination is good for (for example LU
decompositions), but it is general enough to have merit in many
applications.
Top-Level
* (load "gauss.lsp")
; Reading file C:\GCLISP\WORK\GAUSS.LSP
#.(PATHNAME "C:\\GCLISP\\WORK\\GAUSS.LSP")
The main body of the code supports a function called gauss which
uses Gaussian elimination to reduce a matrix to row echelon form. For
reasons I'll mention later, it does so in a manner that does not change
the determinant of the matrix. It expects one argument, the matrix. In
LISP, a matrix can be implemented as a list of rows, where each row is a
list of elements. In this version of gauss, elements must be numbers.
Perhaps in the future I'll update it to allow variables as well, so that
eigenvalues can be computed.
* (gauss '((1 3 5) (2 4 6) (7 8 9)))
((1 3 5) (0 -2.0 -4.0) (0 0 0.0))
The flat notation for a matrix is a little difficult to read, so I
whipped up display to print a matrix row by row. Notice that display
returns the matrix it was given as an argument.
* (display (gauss '((1 3 5) (2 4 6) (7 8 9))))
(1 3 5)
(0 -2.0 -4.0)
(0 0 0.0)
((1 3 5) (0 -2.0 -4.0) (0 0 0.0))
That a little better. Notice that the bottom row is contains all
zeros. As you know, this shows that the determinant is zero. Let's
check.
* (det '((1 3 5) (2 4 6) (7 8 9)))
0.0
Det can be used to find the determinant of any size matrix. All it
does is upper triangulate the matrix with gauss and multiply the
diagonal elements. Since gauss does not change the determinant, this
answer will be correct.
* (det '((2 4 6 1 3)
(3 5 2 3 9)
(4 5 2 1 1)
(3 0 0 5 7)
(1 2 3 4 5)))
1540.0
Also included is a function called gauss_jordan, which completes the
solution finding process by converted the matrix to the identity. It
does so in a very nieve fashon: it flips the matrix over, runs it
through gauss again, then turns it rightside-up. After that dizzying
process, all rows are divided by their first non-zero number to finnish
the conversion. The flipping function only operates on as many columns
as you have rows, so augmented matricies are treated correctly, in most
cases. Let's use this to find the inverse of the monster matrix shown
above.
* (display (gauss_jordan '((2 4 6 1 3 1 0 0 0 0)
(3 5 2 3 9 0 1 0 0 0)
(4 5 2 1 1 0 0 1 0 0)
(3 0 0 5 7 0 0 0 1 0)
(1 2 3 4 5 0 0 0 0 1))))
(1.0 0.0 0.0 0.0 0.0 1.85714F-01 -0.1 0.1 0.3 -3.71428F-01)
(0.0 1.0 0.0 0.0 0.0 -0.21948 1.18182F-01 1.54545F-01 -2.63636F-01
2.57143F-01)
(0.0 0.0 1.0 0.0 0.0 2.48701F-01 -9.54545F-02 -8.63636F-02 5.90909F-02
-4.28571F-02)
(0.0 0.0 0.0 1.0 0.0 -2.21429F-01 -0.15 0.15 -0.05 4.42857F-01)
(0.0 0.0 0.0 0.0 1.0 7.85714F-02 0.15 -0.15 0.05 -1.57143F-01)
((1.0 0.0 0.0 0.0 0.0 1.85714F-01 -0.1 0.1 0.3 -3.71428F-01) (0.0 1.0
0.0 0.0 0.0 -0.21948 1.18182F-01 1.54545F-01 -2.63636F-01 2.57143F-01)
(0.0 0.0 1.0 0.0 0.0 2.48701F-01 -9.54545F-02 -8.63636F-02 5.90909F-02
-4.28571F-02) (0.0 0.0 0.0 1.0 0.0 -2.21429F-01 -0.15 0.15 -0.05
4.42857F-01) (0.0 0.0 0.0 0.0 1.0 7.85714F-02 0.15 -0.15 0.05
-1.57143F-01))
Yeeeach!! Maybe that wasn't such a good idea. If you care to
verify, I'm confident that this answer is correct. Notice that the left
side has been converted to the identity, while the right just tagged
allong step by step, hence changing into the inverse. Let's try with a
nicer matrix, i.e. one I took from the book.
* (display (gauss_jordan '((1 2 3 1 1 0 0 0)
(1 3 3 2 0 1 0 0)
(2 4 3 3 0 0 1 0)
(1 1 1 1 0 0 0 1))))
(1.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0)
(0.0 1.0 0.0 0.0 1.0 -2.0 2.0 -3.0)
(0.0 0.0 1.0 0.0 3.97364F-08 1.0 -1.0 1.0)
(0.0 0.0 0.0 1.0 -2.0 3.0 -2.0 3.0)
((1.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0) (0.0 1.0 0.0 0.0 1.0 -2.0 2.0 -3.0)
(0.0 0.0 1.0 0.0 3.97364F-08 1.0 -1.0 1.0) (0.0 0.0 0.0 1.0 -2.0 3.0
-2.0 3.0))
That's better. To verify, the determinant of the matrix should be
the recipricol of the determinant of the original. Let's assume that
the incredibly small number should actually be a zero.
* (det '((1 2 3 1)
(1 3 3 2)
(2 4 3 3)
(1 1 1 1)))
-1.0
* (det '(( 1 -2 1 0)
( 1 -2 2 -3)
( 0 1 -1 1)
(-2 3 -2 3)))
-1.0
Indeed this is true, so I assume nothing has gone wrong. Now I
shall demonstrate how to use gauss_jordan to solve a linear system.
Give it the augmented matrix for the system:
2x + 4y - 4z = 12
x - 4y + 3z =-21
-6x - 9y +10z =-24.
* (display (gauss_jordan '(( 2 4 -4 12)
( 1 -4 3 -21)
(-6 -9 10 -24))))
(1.0 0.0 0.0 -4.0)
(0.0 1.0 0.0 2.0)
(0.0 0.0 1.0 -3.0)
((1.0 0.0 0.0 -4.0) (0.0 1.0 0.0 2.0) (0.0 0.0 1.0 -3.0))
And behold! The solution x = -4, y = 2, z = -3 is correct. After
all, that's what's in the book. Let's try another inverse.
* (display (gauss_jordan '(( 2 1 0 1 0 0)
(-4 -1 -3 0 1 0)
( 6 2 3 0 0 1))))
(0.0 0.0 0.0 1.0 -1.0 -1.0)
(1.0 0.5 0.0 0.5 0.0 0.0)
(0.0 1.0 -3.0 2.0 1.0 0.0)
((0.0 0.0 0.0 1.0 -1.0 -1.0) (1.0 0.5 0.0 0.5 0.0 0.0) (0.0 1.0 -3.0
2.0 1.0 0.0))
What happened? I told you that my implementation of gauss_jordan
was nieve. When the matrix is singular, this sort of thing will happen.
The top row of the right-hand square matrix, which was the bottom row
for a little while, was zeroed by the second gauss. Hence, we can't use
my gauss_jordan for systems with free variables. The alternatives are
to just gauss the augmented matrix and use manual back-substitution, or
write a better gauss_jordan. Let's verify that the matrix is not
invertable.
* (det '(( 2 1 0)
(-4 -1 -3)
( 6 2 3)))
0.0
Sho'nuf. Thus endeth the demonstration. Please send any comments,
updates, tesimonials on how this code changed you're life, or other neat
stuff to ·······@TAMVENUS.
Mikey Perry.
* (dribble)
___________________________________Cut here____________________________________
(defun insert_column (matrix)
(cond ((null matrix) nil)
(t (cons (cons 0 (car matrix))
(insert_column (cdr matrix))))))
(defun remove_column (matrix)
(cond ((null matrix) nil)
(t (cons (cdar matrix)
(remove_column (cdr matrix))))))
(defun pivot (matrix)
(cond ((or (not (zerop (caar matrix)))
(null (cdr matrix))) matrix)
(t (let ((p (pivot (cdr matrix))))
(cons (negate_row (car p))
(cons (car matrix)
(cdr p)))))))
(defun negate_row (row)
(cond ((null row) nil)
(t (cons (- (car row))
(negate_row (cdr row))))))
(defun add_multiple (m row target)
(cond ((null target) nil)
(t (cons (+ (* m (car row)) (car target))
(add_multiple m (cdr row) (cdr target))))))
(defun reduce_all (matrix row)
(cond ((null matrix) nil)
(t (cons (add_multiple (- (/ (caar matrix) (car row)))
row (car matrix))
(reduce_all (cdr matrix) row)))))
(defun eliminate (matrix)
(cond ((zerop (caar matrix)) (insert_column (slide_right
(remove_column matrix))))
(t (cons (car matrix)
(insert_column (gauss (remove_column (reduce_all
(cdr matrix) (car matrix)))))))))
(defun slide_right (matrix)
(cond ((null (car matrix)) matrix)
(t (eliminate (pivot matrix)))))
(defun gauss (matrix)
(cond ((< (length matrix) 2) matrix)
(t (slide_right matrix))))
(defun flip_square (matrix)
(reverse (flip_n matrix (length matrix))))
(defun flip_n (matrix n)
(cond ((null matrix) nil)
(t (cons (reverse_n n (car matrix))
(flip_n (cdr matrix) n)))))
(defun first_n (n row)
(cond ((or (null row) (< n 1)) nil)
(t (cons (car row) (first_n (- n 1) (cdr row))))))
(defun reverse_n (n row)
(append (reverse (first_n n row)) (nthcdr n row)))
(defun gauss_jordan (matrix)
(unit_diag (flip_square (gauss (flip_square (gauss matrix))))))
(defun det (matrix)
(prod_diag (gauss matrix)))
(defun prod_diag (matrix)
(cond ((null matrix) 1)
(t (* (caar matrix)
(prod_diag (remove_column (cdr matrix)))))))
(defun first_non (row)
(cond ((null row) 0)
((zerop (car row)) (first_non (cdr row)))
(t (car row))))
(defun mul_row (row m)
(cond ((null row) nil)
(t (cons (* (car row) m)
(mul_row (cdr row) m)))))
(defun unit_diag (matrix)
(cond ((null matrix) nil)
(t (cons (mul_row (car matrix) (/ (first_non (car matrix))))
(unit_diag (cdr matrix))))))
(defun display (matrix)
(when (consp matrix)
(format t "~&~S" (car matrix))
(display (cdr matrix))
(format t "~&")
matrix))