練習 Gauss-Jordanの掃き出し法
練習としてLispで連立1次方程式を解いてみた.
Gauss-Jordanの掃き出し法では,三つ行基本変形を繰り返してA=bをE=cに変形して連立1次方程式を解く.ここでEは単位行列で,A=bとE=cは等価な連立一次方程式でる.行基本変形は以下の三つである.
- 1つの行をk倍する(k!=0)
- ある行に別の行のk倍を加える(k!=0)
- 2つの行を入れ換える
1次元行列を第一引数にとり,全要素を第二引数倍する関数arr*kを作成する.
(defun arr*k (arr k)
(let* ((arr-row-num (array-dimension arr 0))
(result-arr (make-array arr-row-num)))
(dotimes (i arr-row-num result-arr)
(setf (aref result-arr i)
(* (aref arr i) k)))))
実行してみる.
> (arr*k #1a(1 2 3) 2) #(2 4 6)arr*kではlet*を使用している.letでは,局所変数を設定する部分で,同じ部分の他の局所変数を参照できない.これはletを入れ子にすることで解決する方法と,let*を使用する方法がある.
> (let ((a 1) (b a)) (format t "a[~A]:b[~A]~%" a b)) *** - EVAL: variable A has no value > (let ((a 1)) (let ((b a)) (format t "a[~A]:b[~A]~%" a b))) a[1]:b[1] NIL > (let* ((a 1) (b a)) (format t "a[~A]:b[~A]~%" a b)) a[1]:b[1] NIL
前節で行をk倍する関数を作成したので,ここではある行に別の行を加える関数add-arrayを作成する.この関数は二つの1次元配列を引数にとり,その和となる1次元行列を返す.
(defun add-array (arr1 arr2)
(let ((arr1-row-num (array-dimension arr1 0))
(arr2-row-num (array-dimension arr2 0)))
(if (/= arr1-row-num arr2-row-num)
(progn (format t "arr1-row-num[~A]:arr2-row-num[~A]~%"
arr1-row-num arr2-row-num)
(return-from add-array)))
(let ((result-arr (make-array arr1-row-num)))
(dotimes (i arr1-row-num result-arr)
(setf (aref result-arr i)
(+ (aref arr1 i) (aref arr2 i)))))))
第一引数と第二引数の要素数の入力チェックを入れてみた.次元の誤りはエラーになるので検出できる.実行してみる.
> (add-array #1a(1 2 3) #1a(4 5 6)) #(5 7 9) > (add-array #1a(1 2 3 4) #1a(4 5 6)) arr1-row-num[4]:arr2-row-num[3] NIL > (add-array #1a(1 2 3) #1a(4 5 6 7)) arr1-row-num[3]:arr2-row-num[4] NIL > (add-array #2a((1 2 3)) #1a(4 5 6)) arr1-row-num[1]:arr2-row-num[3] NIL > (add-array #2a((1) (2) (3))) #1a(4 5 6)) *** - EVAL/APPLY: too few arguments given to ADD-ARRAYこれである行に別の行のk倍を加えるには,以下のようにできる.
> (add-array #1a(1 2 3) (arr*k #1a(4 5 6) 2)) #(9 12 15)
まず2次元配列の指定した行に指定した1次元配列を代入する関数set-rowを作成する.この関数は行をk倍した結果を代入する際など,二つの行の入れ換え以外でも使用する.set-rowは第一引数に代入先の2次元行列,第二引数に代入元となる1次元行列,第三引数に代入先の行を指定する整数を引数にとる.この関数は破壊的関数とし,代入が成功した場合の返り値は(破壊後の)第一引数とする.
(defun set-row (arr1 arr2 i)
(let ((arr1-col-num (array-dimension arr1 1))
(arr2-col-num (array-dimension arr2 0)))
(if (<= arr1-col-num i)
(progn (format t "arr1-col-num[~A]:i[~A]~%" arr1-col-num i)
(return-from set-row)))
(if (/= arr1-col-num arr2-col-num)
(progn (format t "arr1-col-num[~A]:arr2-col-num[~A]~%" arr1-col-num arr2-col-num)
(return-from set-row)))
(dotimes (j arr1-col-num arr1)
(setf (aref arr1 i j) (aref arr2 j)))))
実行してみる.
> (setf a1 (make-array '(3 3))) #2A((NIL NIL NIL) (NIL NIL NIL) (NIL NIL NIL)) > (set-row a1 #1a(1 2 3) 0) #2A((1 2 3) (NIL NIL NIL) (NIL NIL NIL)) > (set-row a1 #1a(4 5 6) 1) #2A((1 2 3) (4 5 6) (NIL NIL NIL)) > (set-row a1 #1a(7 8 9) 2) #2A((1 2 3) (4 5 6) (7 8 9)) > (set-row a1 #1a(10 11) 0) arr1-col-num[3]:arr2-col-num[2] NIL > (set-row a1 #1a(10 11 12) 3) arr1-col-num[3]:i[3] NILこれを使用して,二つの行を入れ換える関数swap-rowを作成した.
(defun swap-row (arr row1 row2)
(let ((temp-row (aref-row arr row1)))
(set-row arr (aref-row arr row2) row1)
(set-row arr temp-row row2)))
この関数は,第一引数のarrのrow1行目とrow2行目を入れ換える破壊的関数である.関数aref-rowは行列計算の頁で作成した関数で,第一引数の二次元行列から第二引数で指定した行を取り出す.実行してみる.
> (setf arr1 #2a((1 0) (0 1))) #2A((1 0) (0 1)) > (swap-row arr1 0 1) #2A((0 1) (1 0)) > arr1 #2A((0 1) (1 0))
いよいよ掃き出し法による連立1次方程式を解く関数gauss-jordanを作成する.
(defun gauss-jordan (arr)
(let ((arr-row-num (array-dimension arr 0))
(arr-col-num (array-dimension arr 1)))
;; 入力チェック
(if (/= (+ arr-row-num 1) arr-col-num)
(progn (format t "arr-row-num[~A]:arr-col-num[~A]~%"
arr-row-num arr-col-num)
(return-from gauss-jordan)))
(let ((result-arr (make-array (list arr-row-num arr-col-num))))
;; 返り値行列領域初期化
(dotimes (i arr-row-num)
(set-row result-arr (aref-row arr i) i))
;; 返り値行列の対角成分が0でないように行を入れ換える
(dotimes (i arr-row-num)
(if (= 0 (aref result-arr i i))
(do* ((j (+ i 1) (+ j 1)) (swap-target (mod j arr-row-num) (mod j arr-row-num )))
((/= (aref result-arr swap-target i) 0)
(swap-row result-arr i swap-target))
(if (= (mod j arr-row-num) i)
(return-from gauss-jordan)))
(set-row result-arr (aref-row result-arr i) i)))
;; 対角成分を除く(i,j)成分を0にする
(dotimes (i arr-row-num)
(dotimes (j arr-row-num)
(if (/= i j)
(let* ((flag (if (> (* (aref result-arr j i)
(aref result-arr i i)) 0)
-1 1)))
(if (and (/= (aref result-arr i i) 0)
(/= (aref result-arr j i) 0))
(set-row result-arr
(add-array (arr*k (aref-row result-arr j) (abs (aref result-arr i i)))
(arr*k (aref-row result-arr i) (* (abs (aref result-arr j i)) flag)))
j))))))
;; 対角成分を1にする
(dotimes (i arr-row-num result-arr)
(if (/= (aref result-arr i i) 0)
(set-row result-arr
(arr*k (aref-row result-arr i) (/ 1 (aref result-arr i i)))
i)))
;; 単位行列か否かをチェックして,解が出たかどうかをチェック
(dotimes (i arr-row-num)
(dotimes (j arr-row-num)
(if (not (or (and (= i j)
(= (aref result-arr i j) 1))
(and (/= i j)
(= (aref result-arr i j) 0))))
(return-from gauss-jordan))))
(return-from gauss-jordan result-arr))))
若干長くなった.この関数はn行(n+1)列の2次元配列を引数にとり,行基本変形を行なって0〜n列目の部分を単位行列とする関数である.破壊的関数ではない.また,この関数では,行列計算の頁で作成した関数aref-rowも使用している.result-arrを初期化後,その対角成分が0にならないようにswap-rowを呼び出している部分は,あまり自身がないが,いくつかのテストをしたところ,問題ないようだ.実行結果のかわりにここまでの関数といくつかのテスト関数を含むファイルを置いておく.