Here is my contribution to this topic... I've decided to change your (Eigenvalue) sub as QR-decomposition is not good for me, and Jakobi and Power method also gave wrong results... I hope you don't mind I hijacked your topic, but you said that any suggestions are welcomed...
;;;-------------------------------------------------------------------------------;;;
;;; Get EigenValues and EigenVectors ;;;
;;;-------------------------------------------------------------------------------;;;
(defun [E][V] (mat / Eigenvalue EigenValue->EigenVectors eea:gauss detm eqmember [I] [¦Ë] [k*] [-] [trp] vslip vzero vrefzero ;; subs
el vl res ;; vars
)
;; get Eigenvalues and Eigenvectors
;;
;; First Use QR-decomposition to get all EigenValues , - bad - checked by M.R. - using old (Eigenvalue) sub instead...
;; Then Use Gaussian-elimination to get all EigenVectors .
;;
;; mat ---- a nonsingular matrix
;;
;; returns a list of all Eigenvalues and Eigenvectors , like '((Value1 Vector1) (Value2 Vector2) ... (Valuen Vectorn))
;;
;; e.g.
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.))) ;_([E][V] a)--> ((1.0 1.0 -1.0 1.0) (2.0 -1.0 0.0 1.0) (4.0 1.0 2.0 1.0))
;;
;; (setq a '((7. 3. -2.) (3. 4. 1.) (-2. -1. 3.))) ;_([E][V] a)--> ((1.72508 2.51661 -3.75831 1.0) (3.0 -1.0 2.0 1.0) (9.27492 -2.51661 -1.24169 1.0))
;;
;; (setq a '((4. 1. 0.) (1. 0. -1.) (1. 1. -4.))) ;_([E][V] a)--> ((-3.7601 -0.0354877 0.275388 1.0) (-0.442931 -1.03315 4.59022 1.0) (4.20303 6.81864 1.38439 1.0))
;;
;; (setq a '((1. 1. 0.5) (1. 1. 0.25) (0.5 0.25 2.))) ;_([E][V] a)--> ((-0.0166473 -7.69469 7.32278 1.0) (1.48012 -0.63687 -0.805775 1.0) (2.53653 0.748221 0.649661 1.0))
;;
;; (setq a '((4.0 2.0 6.0 8.0 2.0) (2.0 10.0 9.0 13.0 7.0) (6.0 9.0 14.0 20.0 12.0) (8.0 13.0 20.0 54.0 35.0) (2.0 7.0 12.0 35.0 43.0)))
;; ;_([E][V] a)--> ((0.0660769 6.41687 3.11375 -4.83309 -0.559052 1.0) (4.18334 4.56872 -4.5105 3.56226 -1.68936 1.0) (7.50996 -0.156679 1.91471 1.45052 -1.88531 1.0) (18.2577 -0.391138 -0.491509 -0.584425 -0.385896 1.0) (94.9829 0.171168 0.328923 0.502966 1.23721 1.0))
;;
;; (setq a '((2.0 3.0 8.0) (1.0 0.0 5.0) (-1.0 5.0 0.0))) ;_([E][V] a)--> ((-4.19258 -0.807418 -1.0 1.0) (1.19258 -6.19258 -1.0 1.0) (5.0 4.58333 1.91667 1.0))
;;
;; (setq a '((-3. 1. -1.) (-7. 5. -1.) (-6. 6. -2.))) ;_([E][V] a)--> ((-2.0 -9.93344e+07 -9.93344e+07 1.0) (-2.0 1.0 1.0 0.0) (4.0 0.0 1.0 1.0))
;;
;; (setq a '((1. 0. 0.) (0. 2. 0.) (0. 0. 3.))) ;_([E][V] a)--> ((1.0 1.0 0.0 0.0) (2.0 0.0 1.0 0.0) (3.0 0.0 0.0 1.0))
;;
;; (setq a '((1. 2.) (2. 1.))) ;_([E][V] a)--> ((-1.0 -1.0 1.0) (3.0 1.0 1.0))
;;
;; (setq a '((1 2) (2 1))) ;_([E][V] a)--> nil - all matrix elements must be reals
;;
;; get: (setq vl ([E][V] a) v (nth 1 vl))
;;
;; checkEg: (equal (detm ([-] ([¦Ë] (length (car a)) (car v)) a)) 0.0 1e-8)
;;
;; checkEv: (equal (mxv a (cdr v)) (mxv ([¦Ë] (length a) (car v)) (cdr v)) 1e-8)
;;
;; by GSLS(SS) 2013-1-2
;; mod by M.R. (Eigenvalue) sub.
;;
;; Any suggestions are welcome , and please post them to chlh_jd@126.com or 56454300@qq.com
;;
;; Thanks Robert for method suggestions .
;; Thanks HighFlyingBird's help .
;; Thanks ElpanovEvgeniy's Gaussian-elimination function , it's really a pretty codes .
;; SUB FUNCTIONS ;;
;; old sub - not QR-decomposition
(defun Eigenvalue (mat / f1 f2 eigr imat breakdown eps mi ma st si er res)
;; To get all Eigenvalues
(defun f1 (mat / i v a r rr)
;; Eigenvalue interval analysis
(setq i 0)
(repeat (length (car mat))
(setq v (nth i mat)
a (apply (function +) (mapcar (function abs) v))
r (cons (list (- (apply (function min) v) a) (+ (apply (function max) v) a)) r)
i (1+ i)
)
)
(setq r (apply (function append) r))
(setq rr (list (cons (apply (function min) r) (apply (function max) r))))
rr
)
;; Iteration
(defun f2 (mi si / chk ei ri mil r11 r21 r12 r22 rr)
(setq chk mi)
(setq ei (list mi)
mil mi
)
(repeat (fix (/ breakdown 2))
(setq mi (+ mi si)
ei (cons mi ei)
)
)
(setq ei (reverse ei))
(repeat (fix (/ breakdown 2))
(setq mil (- mil si)
ei (cons mil ei)
)
)
(setq ri
(mapcar
(function
(lambda (e)
(cons e (abs (detm ([-] ([k*] e imat) mat))))
)
) ei
)
)
(foreach x ri
(if (> (cdr (nth (1+ (vl-position x ri)) ri)) (cdr x))
(setq r11 (cons x r11))
)
)
(foreach x r11
(if (> (cdr (nth (1+ (vl-position x r11)) r11)) (cdr x))
(setq r21 (cons x r21))
)
)
(foreach x (setq ri (reverse ri))
(if (> (cdr (nth (1+ (vl-position x ri)) ri)) (cdr x))
(setq r12 (cons x r12))
)
)
(foreach x r12
(if (> (cdr (nth (1+ (vl-position x r12)) r12)) (cdr x))
(setq r22 (cons x r22))
)
)
(setq rr (append r21 r22))
(if (and (equal (car (last r11)) (car (cadr (reverse r11))) (* 1.1 si)) (< (cdr (last r11)) (cdr (cadr (reverse r11)))) (not (vl-position (last r11) rr)))
(setq rr (cons (last r11) rr))
)
(if (and (equal (car (last r12)) (car (cadr (reverse r12))) (* 1.1 si)) (< (cdr (last r12)) (cdr (cadr (reverse r12)))) (not (vl-position (last r12) rr)))
(setq rr (cons (last r12) rr))
)
(if (null rr)
(setq rr (cons (car (vl-sort ri (function (lambda (x y) (< (cdr x) (cdr y)))))) rr))
)
(if (equal chk (caar rr) 1e-5)
(vl-sort rr (function (lambda (x y) (< (cdr x) (cdr y)))))
rr
)
)
;; MAIN Eigenvalue ;;
(setq eigr (f1 mat)
imat ([I] (length (car mat)))
breakdown 1000 ;_here is a main arg .
eps 1e-8 ;_Tolerance value
)
(foreach a eigr
(setq mi (car a)
ma (cdr a)
st (- ma mi)
mi (/ (+ mi ma) 2.)
si (/ st breakdown)
er (f2 mi si)
)
(while (and (vl-some (function (lambda (x) (> (cdr x) eps))) er) (> si 1e-20))
(setq er
(mapcar
(function
(lambda (x)
(car (f2 (car x) (* 1.1 (/ si breakdown))))
)
) er
)
)
(setq si (/ si breakdown))
)
(foreach a (vl-sort er (function (lambda (x y) (< (cdr x) (cdr y)))))
(if (and (<= (cdr a) 1e-5) (not (eqmember (car a) res Eps)))
(setq res (cons (car a) res))
)
)
)
res
)
(defun EigenValue->EigenVectors (¦Ë mat / n imat matt imatt ires)
;_(setq mat a ¦Ë -2. ires nil)
;_(if (equal ¦Ë 4. 1e-4)(princ))
;_(princ (rtos ¦Ë 2 8))
(setq n (length mat)
imat ([I] n)
mat ([trp] ([-] ([¦Ë] n ¦Ë) mat))
mat (mapcar
(function
(lambda (x y)
(append x y)
)
) mat imat
)
mat (eea:gauss mat)
imat nil
)
(setq mat
(mapcar
(function
(lambda (x / y)
(setq x (reverse x))
(repeat n
(setq y (cons (car x) y)
x (cdr x)
)
)
(setq imat (cons (reverse y) imat))
x
)
) mat
)
)
(setq mat (reverse mat))
(setq matt mat)
(setq imatt imat)
(while (car mat)
(if (<= (apply (function +) (mapcar (function *) (car mat) (car mat))) 5e-7)
(setq ires (cons (car imat) ires))
)
(setq mat (cdr mat)
imat (cdr imat)
)
)
(if (null ires)
(setq mat matt imat imatt)
)
(while (and (cadr mat) (not (vslip (car mat) (cadr mat)))) ;_determine Linealy Independent
(setq
ires (cons (vrefzero (car imat) (cadr imat) (car mat) (cadr mat)) ires)
imat (cdr imat)
mat (cdr mat)
)
)
(mapcar (function reverse)
(mapcar
(function
(lambda (b)
(if (equal b 0.0 5e-7)
0.0
b
)
)
) ires
)
)
)
(defun eea:gauss (lst)
;; Gauss eliminationby
;; by ElpanovEvgeniy
;; http://www.theswamp.org/index.php?action=post;quote=163485;topic=13505.13505.msg163485#msg163485
;; Implementation Gaussian elimination
;; For text:
;; 1x+2y+3z=2
;; 10x+1y+8z=17
;; 7z+2y=5
;; (eea:gauss '((1.0 2.0 3.0) (10.0 1.0 8.0) (0.0 2.0 7.0)))
;; =>
;; ((1.0 2.0 3.0 2.0) (0.0 -19.0 -22.0 -3.0) (0.0 0.0 4.68421 4.68421))
;;(gauss lst)
(if (car lst)
(if ;|(zerop (caar lst))|; (equal (caar lst) 0.0 1e-8) ;_here changed for [E]
(if (vl-every ;|(function zerop)|;
(function (lambda (x) (equal x 0.0 1e-8))) ;_here changed for [E]
(mapcar (function car) lst)
)
(if (cdr lst)
(cons
(car lst)
(mapcar
(function (lambda (x) (cons 0. x)))
(eea:gauss (mapcar (function cdr) (cdr lst)))
) ;_ mapcar
) ;_ cons
lst
) ;_ if
(eea:gauss
(cons
(mapcar
(function +)
(car lst)
(car (vl-remove-if (function (lambda (x) (zerop (car x)))) (cdr lst)))
) ;_ mapcar
(cdr lst)
) ;_ cons
) ;_ gauss
) ;_ if
(cons
(car lst)
(mapcar
(function (lambda (x) (cons 0. x)))
(eea:gauss
(mapcar
(function
(lambda (x / i)
(if (and (equal (car x) 0.0 1e-6)
(equal (caar lst) 0.0 1e-6)
) ;_here add '0 / 0 = 1.'
(setq i 1.)
(setq i (/ (car x) (caar lst)))
)
;|(setq i (/ (car x) (caar lst)))|;
(mapcar
(function -)
(cdr x)
(mapcar (function (lambda (x1) (* x1 i))) (cdar lst))
);_ mapcar
);_ lambda
);_ function
(cdr lst)
);_ mapcar
);_ eea:gauss
);_ mapcar
);_ cons
);_ if
);_ if
);_ defun
(defun detm (m / d)
;; Matrix Determinant - ElpanovEvgeniy
;; Last edit 2013.11.13
;; Args: m - nxn matrix
;; (detm '((0 1) (1 0)))
;; (detm '((1.4 2.1 5.4 6.5) (4.1 9.3 4.5 8.5) (1.2 4.1 6.2 7.5) (4.7 8.5 9.3 0.1)))
(cond ((null m) 1)
((and (zerop (caar m))
(setq d (car (vl-member-if-not (function (lambda (a) (zerop (car a)))) (cdr m))))
)
(detm (cons (mapcar (function +) (car m) d) (cdr m)))
)
((zerop (caar m)) 0)
((* (caar m)
(detm (mapcar (function (lambda (a / d)
(setq d (/ (car a) (float (caar m))))
(mapcar (function (lambda (b c) (- b (* c d)))) (cdr a) (cdar m))
)
)
(cdr m)
)
)
)
)
)
)
(defun eqmember (p l eps)
(cond ((not l) nil)
((if (numberp eps)
(equal p (car l) eps)
(equal p (car l))
)
l
)
((eqmember p (cdr l) eps))
)
)
(defun [I] (d / i r n m)
(setq i d)
(while (<= 0 (setq i (1- i)))
(setq n d
r nil
)
(while (<= 0 (setq n (1- n)))
(setq r (cons (if (= i n) 1.0 0.0) r))
)
(setq m (cons r m))
)
)
(defun [¦Ë] (d ¦Ë / i r n m)
(setq i d)
(while (<= 0 (setq i (1- i)))
(setq n d
r nil
)
(while (<= 0 (setq n (1- n)))
(setq r (cons (if (= i n) ¦Ë 0.0) r))
)
(setq m (cons r m))
)
)
(defun [k*] (k a)
(mapcar
(function
(lambda (x)
(mapcar
(function
(lambda (y)
(* k y)
)
) x
)
)
) a
)
)
(defun [-] (a b)
(mapcar
(function
(lambda (x y)
(mapcar (function -) x y)
)
) a b
)
)
(defun [trp] (a)
(apply (function mapcar) (cons (function list) a))
)
(defun vslip (v1 v2 / z1 z2 e1 sc v)
;; Determine two vectors Linearly Independent .
;; (vslip '(0. 1. -1.) '(0. -1. 1.)) --> NIL
;; (vslip '(0. 1. 1.) '(0. -1. 1.))--> T
(setq z1 (vl-every (function (lambda (x) (equal x 0.0 1e-8))) v1))
(setq z2 (vl-every (function (lambda (x) (equal x 0.0 1e-8))) v2))
(cond
((and z1 z2) nil)
((or z1 z2))
(t
(setq v (vzero v2 v1))
(not (vl-every (function (lambda (x) (equal x 0.0 1e-8))) v))
)
)
)
(defun vzero (vec ref / e1 sc v)
;; Vector Zero
;; (vzero '(0. 1. -1.) '(0. -1. 1.)) --> '(0. 0. 0.)
;; (vzero '(0. 1. 1.) '(0. -1. 1.)) --> '(0. 0. 2.)
(setq e1 (car (vl-sort ref (function (lambda (a b) (> (abs a) (abs b)))))))
(if (not (zerop e1))
(progn
(setq sc (/ (nth (vl-position e1 ref) vec) e1))
(if (zerop sc)
vec
(mapcar
(function
(lambda (a b)
(- a (* b sc))
)
) vec ref
)
)
)
)
)
(defun vrefzero (vec v1 ref ref1 / e1 sc)
;; Vector Zero Refer another two Vectors zero .
(if (vl-every (function (lambda (x) (equal x 0.0 1e-8))) ref)
vec
(progn
(setq e1 (car (vl-sort ref1 (function (lambda (a b) (> (abs a) (abs b))))))
sc (/ (nth (vl-position e1 ref1) ref) e1)
)
(if (zerop sc)
vec
(mapcar
(function
(lambda (a b)
(- a (* b sc))
)
) vec v1
)
)
)
)
)
;; MAIN FUNCTION ;;
(setq el (vl-sort (Eigenvalue mat) (function >)))
(foreach a el
(setq vl (EigenValue->EigenVectors a mat))
(if vl
(foreach c vl
(setq res (cons (cons a c) res))
)
(setq res (cons (cons a nil) res))
)
)
(mapcar
(function
(lambda (x)
(mapcar
(function
(lambda (y)
(if (equal y 0.0 1e-9)
0.0
y
)
)
) x
)
)
)
res
)
)
Regards, M.R.
