自作数値計算プログラム ほげ
私はmatlab環境がないため、schemeで数値計算環境を作り始めた。
一方、世間にはmatlab以外にもscilabやpythonのnumpy,scipy等
すぐれた数値計算ライブラリがある。それらは、すぐれている。
それらを使わずに、数値計算環境を自作するのは、とてつもなく非生産的だ
という自覚はある。
初歩的で恥ずかしいのだが、とりあえず、できたところまで
schemeのコードをブログに書くことにした。
今日のブログも、主な目的は私の備忘録である。
;;行列を表示。デバッグ用に使用。 (define (disp-mat a) (begin (newline) (for-each (lambda (x) (display x) (newline)) a) (newline))) (define (my-map-item2 op a b) (define (my-map-item2-p op a b tmp) (if (and (null? a) (null? b)) tmp (my-map-item2-p op (cdr a) (cdr b) (append tmp (list (op (car a) (car b))))))) (my-map-item2-p op a b '() )) (define (l-accum2 op1 op2 ini a b) (define (l-accum2-p op1 op2 ini a b) (if (and (null? a) (null? b)) ini (l-accum2 op1 op2 (op1 ini (op2 (car a) (car b))) (cdr a) (cdr b)))) (l-accum2-p op1 op2 ini a b)) (define (do-n-times op ini n) (define (do-n-times-p op ini n) (if (= n 0) ini (do-n-times-p op (op ini) (- n 1)))) (do-n-times-p op ini n)) (define (find-1stitem-from-list filter? list-a) (define (find-1stitem-from-list-p filter? counter list-a) (cond ((null? list-a) #f) ;; 見つからなかったときはfalseを出力 ((filter? (car list-a)) (cons counter (car list-a))) ;;(位置. 非ゼロの数)を出力 (else (find-1stitem-from-list-p filter? (+ counter 1) (cdr list-a))))) (find-1stitem-from-list-p filter? 1 list-a)) ;;先頭の位置を1としている (define (zero-line n) (do-n-times (lambda (x) (append x '(0))) '() n)) ;;mXnの零行列 (define (zero-mat m n) (do-n-times (lambda (x) (append x (list (zero-line n)))) '() m)) (define (do-n-times-with-filter-counter op ini n filter? op-c counter) (define (do-n-times-with-filter-counter-p op ini n filter? op-c counter) (if (= n 0) ini (do-n-times-with-filter-counter-p op (op ini filter? counter) (- n 1) filter? op-c (op-c counter)))) (do-n-times-with-filter-counter-p op ini n filter? op-c counter)) (define (eye-line n loc1) (do-n-times-with-filter-counter (lambda (ini filter? counter) (if (filter? counter) (append ini '(1)) (append ini '(0)))) '() n (lambda (counter) (= counter 1)) (lambda (counter) (- counter 1)) loc1)) ;;単位行列 (define (eye-mat n) (do-n-times-with-filter-counter (lambda (ini filter counter) (append ini (list (filter counter)))) '() n (lambda (counter) (eye-line n counter)) (lambda (counter) (+ counter 1)) 1)) ;;和 (define (add-mat a1 a2) (my-map-item2 (lambda (s t) (my-map-item2 + s t)) a1 a2)) ;;スカラー倍 (define (scalar-mat a AA) (define (scalar-each a a1) (map (lambda (x) (* a x)) a1)) (map (lambda (x) (scalar-each a x)) AA)) ;;引き算 (define (sub-mat A1 A2) (add-mat A1 (scalar-mat -1 A2))) (define (sum-of-mul-each a1 a2) (l-accum2 + * 0 a1 a2)) ;;転置行列 (define (transposed-mat a) (apply map list a)) ;;行列の積 (define (mul-mat A1 A2) (define (mul-oneline-mat oneline A2) (map (lambda (x) (sum-of-mul-each oneline x)) (transposed-mat A2))) (define (mul-mat-pre A1 A2 tmp) (if (null? A1) tmp (mul-mat-pre (cdr A1) A2 (append tmp (list (mul-oneline-mat (car A1) A2)))))) (mul-mat-pre A1 A2 '() )) ;;行列のサイズ ;;(行数. 列数) (define (size-mat A) (define (N-mat A j) (if (null? (car A)) j (N-mat (list (cdr (car A))) (+ j 1)))) (define (M-mat A i) (if (null? A) i (M-mat (cdr A) (+ i 1)))) (cons (M-mat A 0) (N-mat A 0))) ;;行列の左右連結 (define (connect-mat-side A B) (define (connect-mat-side-p A B tmp) (if (null? A) tmp (connect-mat-side-p (cdr A) (cdr B) (append tmp (list (append (car A) (car B))))))) (if (= (car (size-mat A)) (car (size-mat B))) (connect-mat-side-p A B '()) (error "CONNECT-MAT-SIDE: (connect-mat-side A B) WRONG SIZE MATRIXES !"))) ;;行列から行列を抽出 (define (extract-mat A i j k l) (define (extract-from-list x1 x2 list-a) (define (extract-from-list-p x1 x2 list-a tmp count) (cond ((> count x2) tmp) ((and (>= count x1) (<= count x2)) (extract-from-list-p x1 x2 (cdr list-a) (append tmp (list (car list-a))) (+ count 1))) ((< count x1) (extract-from-list-p x1 x2 (cdr list-a) tmp (+ count 1))))) (if (and (> x1 0) (>= x2 x1) (>= (length list-a) x2)) (extract-from-list-p x1 x2 list-a '() 1) (error "EXTRACT-FROM-LIST: (extract-from-list x1 x2 list-a) WRONG X1 X2!"))) (define (extract-mat-p D k l tmp) (if (null? D) tmp (extract-mat-p (cdr D) k l (append tmp (list (extract-from-list k l (car D))))))) (extract-mat-p (extract-from-list i j A) k l '())) ;;リストのx番目の要素 (define (x-th-item x list-a) (if (null? list-a) (error "x-th-item") (if (= x 1) (car list-a) (x-th-item (- x 1) (cdr list-a))))) (define (x-th-scaled-list-to-each x list-a A) (define (x-th-scaled-list-to-each-p x list-a A tmp) (if (null? A) tmp (x-th-scaled-list-to-each-p x list-a (cdr A) (append tmp (scalar-mat (/ (x-th-item x (car A)) (x-th-item x list-a)) (list list-a)))))) (x-th-scaled-list-to-each-p x list-a A '())) ;;cdrのdをn個繰り返す (define (cdr*n x n) (define (cdr*n-p x n) (if (= n 0) x (cdr*n-p (cdr x) (- n 1)))) (cdr*n-p x n)) ;;nondiag-biggest-abs ;;入力:対称行列A ;;非対角成分の絶対値が最大となる成分をa_{ij}としたとき ;;出力:(a_{ij} i . j) (define (nondiag-biggest-abs A) ;;内部定義 (define (lbai l r) (define (lbai-p l tmp r c) (if (null? l) tmp (lbai-p (cdr l) (if (> (abs (car l)) (abs (car tmp))) (cons (car l) (cons r c)) tmp) r (+ c 1)))) (lbai-p l (cons 0 0) r (+ r 1))) (define (nondiag-biggest-abs-p A tmp row) (if (= (car (size-mat A)) 1) tmp (nondiag-biggest-abs-p (cdr A) (if (> (abs (car (lbai (cdr*n (car A) row) row))) (abs (car tmp))) (lbai (cdr*n (car A) row) row) tmp) (+ row 1)))) (nondiag-biggest-abs-p A (cons 0 "all-nondiag") 1)) ;;ヤコビ法 ;;入力:実対称行列A, 正の実数tolerance ;;出力:(対角化された行列 . 固有ベクトルを並べた行列) (define (jacobi-method A tolerance) ;;内部定義 (define (norm-plane a b) (sqrt (+ (* a a) (* b b)))) ;;ギブンスの回転行列(準備) ;;入力:行列A ;;出力: '(k m cos(psi) sin(psi)) ;;ただしk,m は、いずれも回転する行と列 (define (Pdat-from-mat A) (let ((nba (nondiag-biggest-abs A))) (let ((k (cadr nba)) (m (cddr nba)) (akm (car nba))) (let ((akk (caar (extract-mat A k k k k))) (amm (caar (extract-mat A m m m m)))) (let ((s (/ (- akk amm) 2.0))) (if (>= s 0) (let ((c-psi (sqrt (/ (+ 1 (/ s (norm-plane s akm))) 2)))) (let ((s-psi (/ akm (* 2 (norm-plane s akm) c-psi)))) (list k m c-psi s-psi))) (let ((c-psi (sqrt (/ (- 1 (/ s (norm-plane s akm))) 2)))) (let ((s-psi (/ (* -1.0 akm) (* 2 (norm-plane s akm) c-psi)))) (list k m c-psi s-psi))))))))) ;;入力: '(k m cos(psi) sin(psi)), 行列A ;;出力: ギブンスの回転行列 (define (make-P-from-Pdat-A Pdat A) (let ((k (car Pdat)) (m (car (cdr*n Pdat 1))) (c-psi (car (cdr*n Pdat 2))) (s-psi (car (cdr*n Pdat 3))) (n (car (size-mat A)))) (let ((m1 (- k 1)) (m2 (- m k 1)) (m3 (- n m))) (let ((els (lambda (lines counter-ini) (do-n-times-with-filter-counter (lambda (ini filter counter) (append ini (list (filter counter)))) '() lines (lambda (counter) (eye-line n counter)) (lambda (counter) (+ counter 1)) counter-ini)))) (let ((em1 (els m1 1)) (em2 (els m2 (+ k 1))) (em3 (els m3 (+ m 1)))) (let ((linek (map + (map (lambda (x) (* c-psi x)) (eye-line n k)) (map (lambda (x) (* -1 s-psi x)) (eye-line n m)))) (linem (map + (map (lambda (x) (* s-psi x)) (eye-line n k)) (map (lambda (x) (* c-psi x)) (eye-line n m))))) (append em1 (list linek) em2 (list linem) em3))))))) (define (jacobi-method-p A tolerance tmpA tmpP) (if (< (abs (car (nondiag-biggest-abs tmpA))) tolerance) (cons tmpA tmpP) ; (disp-mat tmpA) ;;デバッグ用 (jacobi-method-p A tolerance (mul-mat (mul-mat (transposed-mat (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA)) tmpA) (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA)) (mul-mat tmpP (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA)))));) ;;displayデバック時 ;;ここから本体 ヤコビ法 (if (equal? (transposed-mat A) A) (jacobi-method-p A tolerance A (eye-mat (car (size-mat A)))) (error "jacobi-method: (jacobi-method A tolerance) A should be a symetric matrix!!"))) ;べき乗法 (define (power-iter A ini-v counter-limit) (define (power-iter-p A v1 v2 counter counter-limit) (if (>= counter counter-limit) (let ((l1 (car (transposed-mat v1))) (l2 (car (transposed-mat v2)))) (let ((max-eiglist (map (lambda (x1 x2) (/ x2 x1)) l1 l2))) (begin (display "\n絶対値最大の固有値の近似値\n") (display max-eiglist) (display "\n\nそれに対応する固有ベクトル\n") (disp-mat (transposed-mat (list (map (lambda (x y) (/ x (expt y counter))) l2 max-eiglist))))))) (power-iter-p A v2 (mul-mat A v2) (+ counter 1) counter-limit))) (power-iter-p A '() ini-v 0 counter-limit)) (define (first-nonzero list-a) (let ((u (find-1stitem-from-list (lambda (x) ;(not (= x 0))) (> (abs x) 0.00000000000001)) list-a))) (if (not (equal? #f u)) u (cons (+ (length list-a) 1) 0)))) ;;非ゼロが見つからなかったとき、位置をlength+1とする ;;ガウス消去法における前進消去 ;n=列数 (define (gauss-elim-forward tmp-mat-up tmp-mat-low b-line zeros n option) (if (null? tmp-mat-low) (begin (disp-mat (append tmp-mat-up zeros)) (cond ;;出力 (rank . 行操作された行列) ((equal? option "solv") (cons (length tmp-mat-up) (append tmp-mat-up))) ((equal? option "rank") (cons (length tmp-mat-up) (append tmp-mat-up zeros))) (else (error "unknown option")))) ;ここからがelse (if (not (equal? (cdr (first-nonzero (car tmp-mat-low))) 0)) (begin (set! tmp-mat-low (sort tmp-mat-low (lambda (x y) (< (car (first-nonzero x)) (car (first-nonzero y)))))) (set! b-line (map (lambda (a) (/ a (cdr (first-nonzero (car tmp-mat-low))))) (car tmp-mat-low))) ;; normalize (cond ((equal? option "solv") (if (= (car (first-nonzero b-line)) (length b-line)) (error "GAUSS-ELIM-FORWARD: No Solution!"))) ((equal? option "rank") (newline)) (else (error "unknown option"))) (disp-mat (append tmp-mat-up tmp-mat-low zeros)) ;デバッグ用 (set! tmp-mat-up (append tmp-mat-up (list b-line))) (if (not (null? (cdr tmp-mat-low))) (set! tmp-mat-low (sub-mat (cdr tmp-mat-low) (x-th-scaled-list-to-each (car (first-nonzero (car tmp-mat-low))) b-line (cdr tmp-mat-low)))) (set! tmp-mat-low (cdr tmp-mat-low))) (gauss-elim-forward tmp-mat-up tmp-mat-low b-line zeros n option)) (gauss-elim-forward tmp-mat-up (cdr tmp-mat-low) b-line (append zeros (zero-mat 1 n)) n option)))) ;;行列のランク (define (rank A) (car (gauss-elim-forward '() A '() '() (cdr (size-mat A)) "rank"))) ;;連立1次方程式Ax=bをガウスの消去法で解く (define (div-gauss-mat A b) (define (gauss-elim-back tmp-mat-up tmp-mat-low) (if (null? tmp-mat-low) (reverse tmp-mat-up) (begin (set! tmp-mat-up (append tmp-mat-up (list (car tmp-mat-low)))) (if (not (null? (cdr tmp-mat-low))) (set! tmp-mat-low (sub-mat (cdr tmp-mat-low) (x-th-scaled-list-to-each (car (first-nonzero (car tmp-mat-low))) (car tmp-mat-low) (cdr tmp-mat-low)))) (set! tmp-mat-low (cdr tmp-mat-low))) (gauss-elim-back tmp-mat-up tmp-mat-low)))) (let ((Ab (connect-mat-side A b))) (gauss-elim-back '() (reverse (cdr (gauss-elim-forward '() Ab '() '() (cdr (size-mat Ab)) "solv")))))) ;;逆行列 ;;ガウスの消去法で求める ;;2018/01/14 inv 前進消去部 訂正/ 行交換していなかったので行交換するようにした。gauss-elim-forwardを利用 (define (inv A) (define (inv-hakidashi AE) (define (hakidashi-back i tmp-mat-up tmp-mat-low size-mat-AE) (if (> i size-mat-AE) (reverse tmp-mat-up) (begin (set! tmp-mat-up (append tmp-mat-up (list (car tmp-mat-low)))) (if (not (null? (cdr tmp-mat-low))) ;;or (< i size-mat-AE) (set! tmp-mat-low (sub-mat (cdr tmp-mat-low) (x-th-scaled-list-to-each (- (+ size-mat-AE 1) i) (car tmp-mat-low) (cdr tmp-mat-low))))) (hakidashi-back (+ i 1) tmp-mat-up tmp-mat-low size-mat-AE)))) (let ((N (car (size-mat AE)))) (hakidashi-back 1 '() (reverse (cdr (gauss-elim-forward '() AE '() '() (cdr (size-mat AE)) "solv"))) N))) (let ((N (car (size-mat A)))) (extract-mat (inv-hakidashi (connect-mat-side A (eye-mat N))) 1 N (+ N 1) (* 2 N))))
和add-mat , スカラー倍scalar-mat , 引き算sub-mat , 積mul-mat , 転置行列transposed-mat
単位行列eye-mat
この部分の使い方に関しては、ほとんど変わっていない。これが一番最初の記事。作成日記 Schemeで数値計算(1) - nibosiiwasi’s blog
ただ、転置行列に関しては下記の変更を行った。
転置行列のコード - nibosiiwasi’s blog
その後、下記の数値計算法のコードを追加した。
ヤコビ法jacobi-method（実対称行列の固有値計算）
ガウスの消去法
（rank, 連立1次方程式の求解div-gauss-mat , 逆行列 inv）
※ガウスの消去法における、前進消去の際、
first-nonzeroでxがゼロでないことを
(> (abs x) 0.00000000000001)
とした。
浮動小数点の演算の際、厳密にxが0にならない場合があるためである。
逆行列のコードについては、以前のブログ
http://nibosiiwasi.hatenablog.com/entry/2016/10/02/230615
と基本的に同じである。
2018/01/14追記ここから
逆行列の計算に、誤りがありましたので、コードを訂正しました。
申し訳ございません。
前進消去部 訂正/ 行交換していなかったので行交換するようにしました。gauss-elim-forwardを利用しています。
正確数と非正確数の利用状況に応じて
first-nonzeroの
(> (abs x) 0.00000000000001)
を変えた方がいいかもしれません。
2018/01/14追記ここまで
今回、行列Aの逆行列は
(inv A)
と書くことにした。
以下、rankと連立1次方程式の求解の例を挙げる。
行列のランク
例
(define A '((1 2 3 4) (2 1 3 0) (-2 3 1 1))) (rank A)
実行結果
前進消去の計算過程も表示。
(1 2 3 4) (2 1 3 0) (-2 3 1 1) (1 2 3 4) (0 -3 -3 -8) (0 7 7 9) (1 2 3 4) (0 1 1 8/3) (0 0 0 -29/3) (1 2 3 4) (0 1 1 8/3) (0 0 0 1) 3
連立1次方程式の求解
線形方程式Ax=bを解く。ガウスの消去法を使い、行列の形で答えとした。
例
(define A '((1 -1 0) (3 -1 1) (2 -1 2) (0 1 -1))) (define b (transposed-mat '((-2 -2 -1 1)))) (disp-mat (div-gauss-mat A b))
実行結果
前進消去の計算過程も表示。（後退代入は表示させていない。）
(1 -1 0 -2) (3 -1 1 -2) (2 -1 2 -1) (0 1 -1 1) (1 -1 0 -2) (0 2 1 4) (0 1 2 3) (0 1 -1 1) (1 -1 0 -2) (0 1 1/2 2) (0 0 3/2 1) (0 0 -3/2 -1) (1 -1 0 -2) (0 1 1/2 2) (0 0 1 2/3) (0 0 0 0) (1 0 0 -1/3) (0 1 0 5/3) (0 0 1 2/3) #<undef>
本ページのコードのアルゴリズムは全て、下記の本を参考にしました。
・洲之内著,石渡改訂:数値計算,サイエンス社
コードはご自由にお使い下さい。
コードの使用により受けた損害について、著者は一切の責任を負いません。
