Newversion , can get more Eigenvalues and it's Eigenvectors , but has bug yet .
(defun [E] (mat / el b mb vl res d)
;; get Eigenvalues and Eigenvectors
;; by GSLS(SS) 2012-12-28
;;
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.))) -->(4. (0.5 1.0 0.5)) --> ((4.0 (1.0 2.0 1.0)) (1.0 (1.0 -1.0 1.0)) (2.0 (-1.0 0.0 1.0)))
;;
;; (setq a '((7. 3. -2.) (3. 4. 1.) (-2. -1. 3.))) ---> ((9.27492 -2.51661 -1.24169 1.0) (3.0 -1.0 2.0 1.0))
;;
;; (setq a '((4. 1. 0.) (1. 0. -1.) (1. 1. -4.))) --> ((-0.442931 -1.03315 4.59022 1.0) (-3.7601 -0.0354877 0.275388 1.0)) ;_here lost one ?
;; (setq v ([PE] a)) -->(4.20303 (1.0 0.20303 0.146657))
;; ;_ if set EPS = 1e-4 , Only get 1 ? if set EPS = 1e-7 ,Only get 2 ?
;;
;; (setq a '((1. 1. 0.5) (1. 1. 0.25) (0.5 0.25 2.))) --> ((2.53653 0.748221 0.649661 1.0) (1.48012 -0.63687 -0.805775 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)))
;; --> ((94.9829 0.171168 0.328923 0.502966 1.23721 1.0) (18.2577 -0.391138 -0.491509 -0.584425 -0.385896 1.0) (0.0660769 6.41687 3.11375 -4.83309 -0.559052 1.0))
;; ;_ if set EPS = 1e-4 , Only get 3 ? if set EPS = 1e-7 ,Only get 2 ?
;;
;; get: (setq vl ([E] a) v (nth 1 vl))
;; checkEg: (equal (gc:determ ([-] ([k*] (car v) ([I] (length (car a)))) a)) 0.0 1e-7)
;; checkEv (equal (mxv a (cdr v)) (mxv ([λ] (length a) (car v)) (cdr v)) 1e-7)
;;
;; ([E] a)
(setq el (Eigenvalue mat))
(foreach a el
(setq b ([λ] (length (car mat)) a)
mb ([-] mat b)
mb (eea:gauss mb)
vl (UTM->BS mb))
(foreach c vl
(if (equal (mxv mat c) (mxv b c) 1e-7)
(setq res (cons (cons a c) res))
(if (equal (mxv mat (setq d (reverse c))) (mxv b d) 1e-7)
(setq res (cons (cons a d) res)))))
)
(reverse res)
)
;;
(defun UTM->BS (mat / row n i a m)
(setq mat (reverse mat)
row (car mat)
n (length row)
mat (mapcar (function (lambda (r)
(mapcar (function (lambda (a b)
(- a b)))
r
row)))
(cdr mat))
;|
row (car mat)
mat (mapcar (function (lambda (r)
(mapcar (function (lambda (a b)
(- a b)))
r
row)))
(cdr mat))|;
)
(defun f1 (r i n mat / row a)
(foreach row mat
(repeat (- n 1 i)
(setq row (cdr row)))
(setq a (car row)
row (cdr row))
(if (equal a 0.0 1e-13)
(setq r (cons 0.0 r))
(setq row (mapcar (function (lambda (x y)
(* x y)))
r
row)
r (cons (/ (apply (function +) row) (- a)) r))
)
(setq i (1+ i))
)
r)
(vl-remove-if
(function
(lambda (x)
(vl-every (function (lambda (y) (equal y 0.0 1e-9)))
x)))
(list (f1 (list 0.0) 1 n mat)
(f1 (list 1.0) 1 n mat)))
)
(defun [PE] (mat / v0 f0 vk f1)
;; get main Eigenvalue and it's Eigenvector for Real_positive matrix
;; by GSLS(SS) 2012-12-24
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.))) -->(4. (0.5 1.0 0.5))
;; (setq a '((7. 3. -2.) (3. 4. 1.) (-2. -1. 3.))) --> (9.27492 (1.0 0.493399 -0.39736))
;; (setq a '((4. 1. 0.) (1. 0. -1.) (1. 1. -4.))) --> (4.20303 (1.0 0.20303 0.146657))
;; (setq a '((1. 1. 0.5) (1. 1. 0.25) (0.5 0.25 2.))) --> (2.53652587 (0.74822116 0.64966115 1.0))
;; (setq b ([E] a) b (list (rtos (car b) 2 8) (mapcar (function (lambda (x) (rtos x 2 8))) (cadr b))))
;; ([PE] a)
(setq v0 (mapcar (function (lambda (x) 1.0)) (car mat))
f0 1.0
vk (mxv mat v0)
f1 (car (vl-sort vk
(function (lambda (e1 e2)
(> (abs e1) (abs e2))))))
vk (pt* vk (1/ f1)))
(while
(not
(and
(<= (abs (- f1 f0)) 1e-8)
(<=
((lambda (x)
(sqrt (apply (function +) (mapcar (function *) x x))))
(pt- vk v0))
1e-8)))
(setq v0 vk
f0 f1
vk (mxv mat v0)
f1 (car (vl-sort vk
(function (lambda (e1 e2)
(> (abs e1) (abs e2))))))
vk (pt* vk (1/ f1))))
(cons f1 vk))
;;;-----------------
(defun Eigenvalue (mat / f1 f3 eigr imat mi ma st si er res i t1)
;; 矩阵特征值求解
;; To get all Eigenvalues
;; (Eigenvalue a)
(defun f1 (mat / i v a r b c r1)
;; 矩阵特征值区间分析
;; Eigenvalue interval analysis
;; (f1 a)
(setq i 0)
(repeat (length (car mat))
(setq v (nth i mat)
a (nth i v)
v (- (apply (function +) (mapcar (function abs) v))
(abs a))
r (cons (cons (- a v) (+ a v)) r)
i (1+ i)))
(setq r (vl-sort r
(function (lambda (e1 e2)
(if (= (car e1) (car e2))
(< (cdr e1) (cdr e2))
(< (car e1) (car e2))
)))))
(setq a (car r)
r (cdr r))
(while r
(setq b (car r)
r (cdr r))
(if (> (cdr a) (car b))
(progn (if (and r1 (< (car b) (cdar r1)))
(setq c (/ (+ (cdr a) (cdar r1)) 2.))
(setq c (/ (+ (cdr a) (car b)) 2.)))
(setq a (cons (car a) c)
b (cons c (cdr b))
)))
(setq r1 (cons a r1)
a b))
(setq r1 (cons a r1))
(reverse r1))
;;-----------------Iteration
(defun f3 (mi si / ei ri mil)
(setq ei (list mi)
mil mi)
(repeat 10
(setq mi (+ mi si)
ei (cons mi ei)))
(setq ei (reverse ei))
(repeat 10
(setq mil (- mil si)
ei (cons mil ei)))
(setq ri
(mapcar
(function
(lambda (e)
(cons e (abs (gc:determ ([-] ([k*] e imat) mat))))))
ei))
(setq ri (vl-sort ri
(function (lambda (x y)
(< (cdr x) (cdr y))))))
(setq ri (list (car ri)
(cadr ri)
(caddr ri)
(cadddr ri)
(car (cddddr ri))))
(setq si (/ si 10.))
(setq ri
(mapcar
(function
(lambda (x / mi ei)
(setq mi (car x)
mi (- mi (* 5. si))
ei (list mi))
(repeat 10
(setq mi (+ mi si)
ei (cons mi ei)))
(car
(vl-sort
(mapcar
(function
(lambda (e)
(cons
e
(abs
(gc:determ ([-] ([k*] e imat) mat))))))
ei)
(function (lambda (x y)
(< (cdr x) (cdr y))))))))
ri))
(setq ri (vl-sort ri
(function (lambda (x y)
(< (cdr x) (cdr y))))))
(list (car ri) (cadr ri)))
(setq t1 (getvar "millisecs"))
(setq eigr (f1 mat)
imat ([I] (length (car mat))))
(foreach a eigr
(setq mi (car a)
ma (cdr a)
st (- ma mi)
mi (/ (+ mi ma) 2.)
si (/ st 20.))
(setq er (f3 mi si))
(setq i 0)
(while (and (not (< (abs (cdar er)) 1e-7)) (< i 100))
(setq st (* 2.0 si)
mi (caar er)
si (/ st 20.))
(setq er (f3 mi si))
(setq i (1+ i)))
(if (and (<= (cdar er) 1e-7)
(not (eqmember (caar er) res 1e-7)))
(setq res (cons (caar er) res)))
(if (and (<= (cdadr er) 1e-7)
(not (eqmember (caadr er) res 1e-7)))
(setq res (cons (caadr er) res)))
)
(princ
(strcat (rtos (- (getvar "millisecs") t1) 2 0) "msecs"))
res
)
;;
(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))
(if (vl-every (function zerop) (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)
(setq i (/ (car x) (caar lst)))
(mapcar
(function -)
(cdr x)
(mapcar (function (lambda (x1) (* x1 i)))
(cdar lst))
) ;_ mapcar
) ;_ lambda
) ;_ function
(cdr lst)
) ;_ mapcar
) ;_ test
) ;_ mapcar
) ;_ cons
) ;_ if
) ;_ if
) ;_ defun
;;;
;;Use function--------------------------------------------------------------------;;;
;;------------------------------------------Determinant
(defun gc:determ (mat / col piv row res d)
;; get the determinant of a matrix
;; Use Gauss method
;; by Gile
;;
(setq res 1.)
(while (< 2 (length mat))
(setq col (mapcar (function (lambda (x) (abs (car x)))) mat))
(repeat (vl-position (apply (function max) col) col)
(setq mat (append (cdr mat) (list (car mat))))
)
(if (equal (setq piv (float (caar mat))) 0.0 1e-13)
(setq mat nil
res 0.0)
(setq row (mapcar (function (lambda (x) (/ x piv))) (car mat))
mat (mapcar
(function
(lambda (r / e)
(setq e (car r))
(cdr (mapcar
(function (lambda (x n) (- x (* n e))))
r
row))
))
(cdr mat))
res (* piv res)
)))
(setq d (- (* (caar mat) (cadadr mat)) (* (caadr mat) (cadar mat))))
(* res d))
;;;
(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))
))
;;; Unit matrix
(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 [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 [λ] (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 [0 (n / r)
(repeat n
(setq r (cons 0.0 r))
)
)