(defun [E] (mat / v0 f0 vk f1)
;;
;; mat : the given matrix , EPS : 1e-8
;; 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))
;; ([E] a)
(setq v0 (mapcar (function (lambda (x) 1.0)) (car mat))
f0 1.0
vk (mxv mat v0)
f1 (apply (function max) (mapcar (function abs) vk))
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 (apply (function max) (mapcar (function abs) vk))
vk (pt* vk (1/ f1))))
(list f1 vk))
;;----------------
(defun pt* (v c)
(mapcar (function (lambda (x)
(* x c)
)
)
v
)
)
;;----------------
(defun mxv(m v)(mapcar(function(lambda (r)(vxv r v)))m))
(defun vxv (v1 v2)
(apply
(function +)
(mapcar
(function *)
v1
v2
)
)
)
;;
(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
)
;;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))
If has got all Eigenvalues , how get their Eigenvectors ?
If just consider 3 dimension ,then it will get a better way:Thanks highflyingbird ! :-)Code - Auto/Visual Lisp: [Select]Other code see the attachment.
;;;----------------------------------------------------; ;;;EigenValues in n^3 ; ;;;Input: mat -- a 3x3 matrix ; ;;;Output: the EigenValues of this matrix ; ;;;----------------------------------------------------; (Math:Cubic_Equation 1 m n (- z)) )
I have use Matlab to solve SVD because my lisp skills are to low...
http://www.theswamp.org/index.php?topic=40495.msg457915#msg457915
if that will help you here some more information about calculating eigenvalues and eigenvectors of a non-symmetric matrix
http://www.alglib.net/eigen/nonsymmetric/nonsymmetricevd.php
"Solving a non-symmetric problem of finding eigenvalues is performed in some steps. In the first step, the matrix is reduced to upper Hessenberg form by using an orthogonal transformation. In the second step, which takes the most amount of time, the matrix is reduced to upper Schur form by using an orthogonal transformation. If we only have to find the eigenvalues, this step is the last because the matrix eigenvalues are located in the diagonal blocks of a quasi-triangular matrix from the canonical Schur form. If we have to find the eigenvectors as well, it is necessary to perform a backward substitution with Schur vectors and quasi-triangular vectors (in fact - solving a system of linear equations; the process of backward substitution itself takes a small amount of time, but the necessity to save all the transformations makes the algorithm twice as slow)."
Best Regards
Robert
I think the next step is to get the Eigenvectors.Thanks , I 'll search and study 'use Gaussian elimination ' . :-)
Maybe use Gaussian elimination。
Thanks , I 'll search and study 'use Gaussian elimination ' . :-)
Thanks highflyingbird ! :-)Thanks , I 'll search and study 'use Gaussian elimination ' . :-)
See here (http://en.wikipedia.org/wiki/Gaussian_elimination).
Or
http://www.theswamp.org/index.php?topic=13505.msg163485#msg163485
http://www.theswamp.org/index.php?topic=13505.msg163485#msg163485
(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))
)
)
(setq mat
'(
(2.0 3.0 8.0)
(1.0 0.0 5.0)
(-1.0 5.0 0.0)
)
)
(setq
a (car (nth 0 mat)) ;00
b (cadr (nth 0 mat)) ;01
c (caddr (nth 0 mat)) ;02
d (car (nth 1 mat)) ;10
e (cadr (nth 1 mat)) ;11
f (caddr (nth 1 mat)) ;12
g (car (nth 2 mat)) ;20
h (cadr (nth 2 mat)) ;21
i (caddr (nth 2 mat)) ;22
)
(setq
am -1
bm (+ a e i)
cm (- (+ (* d b) (* g c) (* f h)) (* a e) (* a i) (* e i))
dm (- (+ (- (* a e i) (* a f h) (* d b i)) (* d c h) (* g b f)) (* g c e))
)
(progn
(setq x (/ (- (/ (* 3.0 cm) am) (/ (expt bm 2) (expt am 2))) 3.0))
(setq y (/ (+ (- (/ (* (expt bm 3) 2.0) (expt am 3))
(/ (* 9.0 bm cm) (expt am 2))
)
(/ (* 27.0 dm) am)
)
27.0
)
)
(setq z (+ (/ (expt y 2) 4.0) (/ (expt x 3) 27.0)))
(setq i2 (sqrt (- (/ (expt y 2) 4.0) z)))
;; (setq i (expt (- (/ (expt y 2) 4.0) z) 0.5))
(setq j (expt i2 (/ 1.0 3.0)))
;; ArcCosine - Lee Mac
;; Args: -1 <= x <= 1
;; Return: Degree in radians
(defun acos (x)
(if (<= -1.0 x 1.0)
(atan (sqrt (- 1.0 (expt x 2))) x)
)
)
(setq k (acos (* (/ y (* 2.0 i2)) -1.0)))
(setq m (cos (/ k 3.0)))
(setq n (* (sqrt 3.0) (sin (/ k 3.0))))
(setq p (* (/ bm (* 3.0 am)) -1.0)) ;(setq p (/ (* b -1.0) (* 3.0 a)))
;;; (setq l (* j -1.0))
;;; (setq eig1 (- (* (* 2.0 j) (cos (/ k2 3.0))) (/ b (* 3 a))))
;;; (setq eig2 (+ (* l (+ m n)) p))
;;; (setq eig3 (+ (* l (- m n)) p))
(setq eig1 (+ (* (* 2.0 j) m) p))
(setq eig2 (+ (* (* j -1.0) (+ m n)) p))
(setq eig3 (+ (* (* j -1.0) (- m n)) p))
)
(defun chp (eig /)
(setq mat
(list
(list (- a eig) b c)
(list d (- e eig) f)
(list g h (- i eig))
)
)
)
(setq result 0.0)
(defun m_elim (A / T1 T2 X U)
(cond
((= (length (car A)) 2) A)
(t
(setq T1 (apply 'max (mapcar '(lambda (X) (abs (car X))) A))
T1 (if (assoc T1 A) (assoc T1 A) (assoc (* -1 T1) A))
T2 (car A)
)
(if (not (equal T1 T2))
(setq A (cdr A)
A (subst T2 T1 A)
A (cons T1 A)
)
)
(foreach U (cdr A)
(setq T2 (/ (car U) (car T1))
T2 (mapcar '(lambda (T3 T4) (- T4 (* T3 T2))) T1 U)
A (subst T2 U A)
)
)
(cons T1
(mapcar
'(lambda (T2)
(cons 0.0 T2))
(m_elim (mapcar 'cdr (cdr A)))))
)
)
)
(defun m_backsub (A / U V)
(cond
((= (length (car A)) 2)
(list (/ (cadr (car A)) (car (car A))))
)
(t
(setq U (car A)
V (m_backsub (mapcar 'cdr (cdr A)))
W (/ (- (last U) (apply '+ (mapcar '* V (cdr U)))) (car U))
)
(cons W V)
)
)
)
(defun gauss (A / B C)
(foreach tmp A
(prompt (strcat "\n"
(rtos (car tmp))
))
(foreach tmp2 (cdr tmp)
(prompt (strcat "\t" (rtos tmp2))))
)
(setq B (m_elim A))
(prompt "\n\nElimination Matrix result")
(foreach tmp B
(prompt (strcat "\n"
(rtos (car tmp))
))
(foreach tmp2 (cdr tmp)
(prompt (strcat "\t" (rtos tmp2))))
)
(setq C (m_backsub B))
(prompt "\nBack substitution result")
(print C)
(setq result C)
(princ)
)
(setq eigm
(list
(+ eig1 0)
(+ eig2 0)
(+ eig3 0)
)
)
(setq
eigvec1 nil
eigvec2 nil
eigvec3 nil
)
(defun vl-sort-all (lst func)
(mapcar '(lambda (x) (nth x lst)) (vl-sort-i lst func))
)
(setq seigm (vl-sort-all eigm '>))
;; Matrix x Scalar - Lee Mac
;; Args: m - nxn matrix, n - real scalar
(defun mxs ( m s )
(mapcar '(lambda ( r ) (mapcar '(lambda ( n ) (* n s)) r)) m)
)
(chp (nth 0 seigm))
(gauss mat)
(setq eigvec1 (append (list (append result (list 1.0)))))
(chp (nth 1 seigm))
(gauss mat)
(setq eigvec2 (append (list (append result (list 1.0)))))
(chp (nth 2 seigm))
(gauss mat)
(setq eigvec3 (append (list (append result (list 1.0)))))
;norm1
(setq temp1 (/ 1 (sqrt (+ (expt (car (nth 0 eigvec1)) 2) (expt (cadr (nth 0 eigvec1)) 2) (expt (caddr (nth 0 eigvec1)) 2)))))
(setq v1 (mxs (list (nth 0 eigvec1)) temp1))
;end norm1
;norm2
(setq temp2 (/ 1 (sqrt (+ (expt (car (nth 0 eigvec2)) 2) (expt (cadr (nth 0 eigvec2)) 2) (expt (caddr (nth 0 eigvec2)) 2)))))
(setq v2 (mxs (list (nth 0 eigvec2)) temp2))
;end norm2
;norm3
(setq temp3 (/ 1 (sqrt (+ (expt (car (nth 0 eigvec3)) 2) (expt (cadr (nth 0 eigvec3)) 2) (expt (caddr (nth 0 eigvec3)) 2)))))
(setq v3 (mxs (list (nth 0 eigvec3)) temp3))
;end norm3
(print v1)
(print v2)
(print v3)
Hi, here some code i have made last year, its compute eigenvector of 3x3 matrixCode: [Select](setq mat
Result:
'(
(2.0 3.0 8.0)
(1.0 0.0 5.0)
(-1.0 5.0 0.0)
)
)
((-0.904439 -0.37822 0.197332))
((0.974901 0.15743 0.15743))
((0.495813 0.614072 0.614072))
Matlab Result:
>> [V,D] = eig(A)
V =
0.9749 -0.4958 0.9044
0.1574 -0.6141 0.3782
-0.1574 0.6141 0.1973
D =
1.1926 0 0
0 -4.1926 0
0 0 5.0000
(Eigenvalue mat)-->(5.0 1.19258)
([E] mat) --> ((5.0 4.58333 1.91667 1.0) (1.19258 -6.19258 -1.0 1.0))
Ei 5.0 --> EV (4.58333 1.91667 1.0) --> (0.9044 0.3782 0.19732)(setq mat
'(
(2.0 3.0 8.0)
(1.0 0.0 5.0)
(-1.0 5.0 0.0)
)
)
([E] mat) -->((5.0 0.904439 0.37822 0.197332) (-4.19258 -0.495813 -0.614072 0.614072) (1.19258 -0.974901 -0.15743 0.15743))
[/ code]
(defun [E] (mat / el b mb vl res)
;; 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 ?
;; (setq a '((2.0 3.0 8.0) (1.0 0.0 5.0) (-1.0 5.0 0.0)))
;;
;; (setq a '((-3. 1. -1.) (-7. 5. -1.) (-6. 6. -2.)))
;;
;; get: (setq vl ([E] a)) (setq v (nth 3 vl))
;; (mapcar (function (lambda (x) (rtos x 2 10) ))(Eigenvalue a 25 1e-7)) --> ("18.2577195811" "7.5099592408" "4.1833434321")
;; checkEg: (equal (gc:determ ([-] ([λ] (length (car a)) (car v)) a)) 0.0 1e-6)
;;
;; checkEv (equal (mxv a (cdr v)) (mxv ([λ] (length a) (car v)) (cdr v)) 1e-6)
;;
;; ([E] a)
(setq el (Eigenvalue mat 20 1e-7))
(foreach a el
(setq vl (TUM->BS a mat))
(foreach c vl
(setq res (cons (cons a c) res)))
)
(reverse res)
)
;;----
(defun UTM->BS (λ mat / n imat)
(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)
imat (reverse imat))
(mapcar
(function reverse)
;|
(mapcar
(function v2u)|;
(vl-remove
nil
(mapcar
(function
(lambda (x y)
(if (vl-every (function (lambda (a)
(equal a 0.0 1e-7)))
x)
y)))
mat
imat))) ;|)|;
)
;;;-------------------------------------------------------------------------------;;;
;;; Get EigenValues and EigenVectors ;;;
;;;-------------------------------------------------------------------------------;;;
(defun [E] (mat / em i ai aii el vl res)
;; get Eigenvalues and Eigenvectors
;;
;; First Use QR-decomposition to get all EigenValues ,
;; 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] a)--> ((4.0 1.0 2.0 1.0) (2.0 -1.0 0.0 1.0) (1.0 1.0 -1.0 1.0))
;;
;; (setq a '((7. 3. -2.) (3. 4. 1.) (-2. -1. 3.))) ;_([E] a)--> ((9.27492 -2.51661 -1.24169 1.0) (3.0 -1.0 2.0 1.0) (1.72508 2.51661 -3.75831 1.0))
;;
;; (setq a '((4. 1. 0.) (1. 0. -1.) (1. 1. -4.))) ;_([E] a)--> ((4.20303 6.81864 1.38439 1.0) (-3.7601 -0.0354877 0.275388 1.0) (-0.442931 -1.03315 4.59022 1.0))
;;
;; (setq a '((1. 1. 0.5) (1. 1. 0.25) (0.5 0.25 2.))) ;_([E] a)--> ((2.53653 0.748221 0.649661 1.0) (1.48012 -0.63687 -0.805775 1.0) (-0.0166473 -7.69468 7.32278 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] a)--> ((94.9829 0.171168 0.328923 0.502966 1.23721 1.0) (18.2577 -0.391138 -0.491509 -0.584425 -0.385896 1.0) (7.50996 -0.156679 1.91471 1.45052 -1.88531 1.0)
;; (4.18334 4.56872 -4.5105 3.56226 -1.68936 1.0) (0.0660769 6.41687 3.11375 -4.83309 -0.559052 1.0))
;;
;; (setq a '((2.0 3.0 8.0) (1.0 0.0 5.0) (-1.0 5.0 0.0))) ;_([E] a)--> ((5.0 4.58333 1.91667 1.0) (-4.19258 -0.807418 -1.0 1.0) (1.19258 -6.19258 -1.0 1.0))
;;
;; (setq a '((-3. 1. -1.) (-7. 5. -1.) (-6. 6. -2.))) ;_([E] a)--> ((4.0 0.0 1.0 1.0) (-2.0 -1.0 -1.0 0.0) (-2.0 -1.0 -1.0 0.0))
;;
;; get: (setq vl ([E] a) v (nth 1 vl))
;;
;; checkEg: (equal (gc:determ ([-] ([λ] (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
;;
;; 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 .
(setq Em (QR->EV mat 1e-14));_1e--14 The tolerance between the matrix modules with the product of it's Trace ,
;_You can set a larger value like 1e-9 to reduce computations .
(setq i 0)
(while (and (setq ai (nth i Em)) (setq aii (nth i ai)))
(setq el (cons aii el)
i (1+ i)))
(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)
)
;;
(defun EigenValue->EigenVectors (λ mat / n imat 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))
(while (<= (apply (function +)
(mapcar (function *) (car mat) (car mat)))
5e-7)
(setq mat (cdr mat)
ires (cons (car imat) ires)
imat (cdr imat)))
(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 householder (mat / imat i m n d a b u hk q)
;; Use Householder transformation to Get QR matrix for QR-decomposition
;; (householder mat)
;; (setq mat '((2. 3.) (1. 2.) (5. 7.)))
;; (setq mat '((1. 2. -3.) (1. 3. 10.) (1. -1. 6.)))
(setq i 0
n (length (car mat))
m (length mat)
imat ([i] m)
q imat)
(while (< i n)
(setq d (mapcar (function (lambda (x) (nth i x))) mat))
(repeat i
(setq d (cdr d)))
(setq a (* (sign (car d))
(sqrt (apply (function +)
(mapcar (function (lambda (x)
(* x x)))
d))))
b (* a (+ a (car d)))
d (cons (+ (car d) a)
(cdr d)))
(repeat i
(setq d (cons 0. d)))
(setq u (vtxv d d))
(if (zerop b)
nil
(setq u ([k*] (1/ b) u)))
(setq hk ([-] imat u)
q ([*] q hk)
mat ([*] hk mat)
i (1+ i)))
(list q mat))
;;
(defun QR->EV (mat eps / d i b t1)
;; QR-decomposition
;; by GSLS(SS) 2013-1-2
(setq t1 (getvar "millisecs")
d (gc:determ mat)
i 0)
(while (and (not (equal d ([tr] mat) eps)) (< i 1000))
(setq b (householder mat)
mat ([*] (cadr b) (car b))
i (1+ i)))
(princ (strcat "\nQR-decomposition of iterations is : "
(rtos i 2)
" times ."))
(princ (strcat "\nTotel time of QR-decomposition : "
(rtos (- (getvar "millisecs") t1) 2 0)
" m.s."))
mat)
;;
(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
);_ test
);_ mapcar
);_ cons
);_ if
);_ if
);_ defun
;;;
Nice to see that you made it and thanks for share :)Thanks Robert ! You are always welcome. :-)
;;;---------------------------------------幂法-----------------------------------;;;
(defun PowerEigen (mat eps / v0 f0 vk f1)
;; 求矩阵第一特征值与特征向量求解
;; Get First EigenValue and it's EigenVector
;; mat 为方阵 ; eps 为精度、该精度不宜太小否收敛很慢, 宜为1e-6~1e-8
;; by GSLS(SS) 2012-12-24
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.))) -->(4. (0.5 1.0 0.5)) 即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))
;; (PowerEigen a 1e-9)
(setq v0 (mapcar (function (lambda (x) 1.0)) (car mat))
f0 1.0
vk (mxv mat v0)
f1 (apply (function max) (mapcar (function abs) vk))
vk (pt* vk (1/ f1)))
(while
(not
(and
(<= (abs (- f1 f0)) eps)
(<=
((lambda (x)
(sqrt (apply (function +) (mapcar (function *) x x))))
(pt- vk v0))
eps)))
(setq v0 vk
f0 f1
vk (mxv mat v0)
f1 (apply (function max) (mapcar (function abs) vk))
vk (pt* vk (1/ f1))))
(cons f1 vk))
;;
(defun Power->EV (mat eps / m v v1 eig res)
;; Get all EigenValues
;; by GSLS(SS) 2013-1-3
;; (PowerEigenValues a)
(while (cadr mat)
(setq eig (powereigen mat eps)
res (cons eig res))
(if res
(setq v (car mat)
v1 (cdar res)
v1 (pt* v1 (1/ (car v1)))
m (vtxv v1 v)
mat ([-] mat m)
mat (mapcar (function cdr) (cdr mat))
))
)
(reverse (mapcar (function car) (cons (car mat) res)))
)
;;--------------------------------Givens变换法-----------------------------------;;;
(defun Givens (mat Q / i j aij aii ajj aji ang a b m n0 n r Q1)
;;(Givens a)
;; Givens 变换求上三角阵
(setq i (length mat)
m0 i
n0 (length (car mat))
j -1)
(while (< j (1- n0))
(setq i m0
j (1+ j))
(while (> (setq i (1- i)) j)
(setq aij (nth j (nth i mat)))
(if (not (equal aij 0.0 1e-14))
(progn
(setq aii (nth i (nth i mat))
ajj (nth j (nth j mat))
aji (nth i (nth j mat)))
(if (= aii ajj)
(if (> aij 0)
(setq ang _pi4)
(if (< aij 0)
(setq ang (- _pi4))))
(setq ang (/ (atan (/ aij (- ajj aii) 0.5)) 2.))
)
(setq ang (max (min ang _pi4) (- _pi4)))
(setq m m0
Q1 nil)
(while (<= 0 (setq m (1- m)))
(setq n n0
r nil)
(while (<= 0 (setq n (1- n)))
(setq r (cons (cond ((and (= m i) (= n j))
(- (sin ang)))
((and (= m j) (= n i))
(sin ang))
((and (= m i) (= n i))
(cos ang))
((and (= m j) (= n j))
(cos ang))
((= m n)
1.)
(0.0))
r)))
(setq Q1 (cons r Q1))
)
(setq Q ([*] Q1 Q))
(setq mat ([*] ([*] Q1 mat) ([trp] Q1)))
)
)
)
)
(list mat Q))
;;-------------------
(defun Givens->EV (mat eps / d i b t1 q )
;; Givens method main function
;; by GSLS(SS) 2013-1-4
;; (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)))
;; (Givens->EV a 1e-12) --> (0.0660769 4.18334 18.2577 94.9829 7.50996)
(setq t1 (getvar "millisecs")
d (gc:determ mat)
q ([i] (length (car mat)))
i 0 ;_Set the maximum number of iterations to prevent stack overflow
)
(while (and (not (equal d ([tr] mat) eps)) (< i 1000))
(setq b (Givens mat Q)
mat (car b)
q (cadr b))
(setq i (1+ i)))
(print mat)
(princ (strcat "\nJacobi iterations is : "
(rtos i 2)
" times ."))
(princ (strcat "\nTotel time of Jacobi solving : "
(rtos (- (getvar "millisecs") t1) 2 0)
" m.s."))
(setq i -1)
(mapcar (function (lambda (r1 )
(nth (setq i (1+ i)) r1)))
mat)
)
(defun Givens (mat Q / i j aij aii ajj aji ang a b m n0 n r Q1)
;;(Givens a)
;; Givens 变换求上三角阵
(setq i (length mat)
m0 i
n0 (length (car mat))
j -1)
(while (< j (1- n0))
(setq i m0
j (1+ j))
(while (> (setq i (1- i)) j)
(setq aij (nth j (nth i mat)))
(if (not (equal aij 0.0 1e-14))
(progn
(setq aii (nth i (nth i mat))
ajj (nth j (nth j mat))
aji (nth i (nth j mat)))
(if (= aii ajj)
(if (> aij 0)
(setq ang _pi4)
(if (< aij 0)
(setq ang (- _pi4))))
(setq ang (/ (atan (/ aij (- ajj aii) 0.5)) 2.))
)
(setq ang (max (min ang _pi4) (- _pi4)))
(setq m m0
Q1 nil)
(while (<= 0 (setq m (1- m)))
(setq n n0
r nil)
(while (<= 0 (setq n (1- n)))
(setq r (cons (cond ((and (= m i) (= n j))
(- (sin ang)))
((and (= m j) (= n i))
(sin ang))
((and (= m i) (= n i))
(cos ang))
((and (= m j) (= n j))
(cos ang))
((= m n)
1.)
(0.0))
r)))
(setq Q1 (cons r Q1))
)
(setq Q ([*] Q1 Q))
(setq mat ([*] ([*] Q1 mat) ([trp] Q1)))
)
)
)
)
(list mat Q))
;;-------------------------------------------------------------------------------;;;
(defun Jacobi (mat eps / d i b t1 q)
;; Jacobi method
;; by GSLS(SS) 2013-1-4
;;
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.)))--> ((1.0 0.57735 -0.57735 0.57735) (4.0 0.408248 0.816497 0.408248) (2.0 -0.707107 -2.59319e-010 0.707107))
;;
;; (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))) -->
;; -->((0.0660769 0.738289 0.358251 -0.556068 -0.0643214 0.115054) (4.18334 -0.60114 0.59348 -0.468712 0.222282 -0.131577) (18.2577 0.284885 0.35799 0.425666 0.281067 -0.72835)
;; (94.9829 0.100149 0.192449 0.29428 0.723881 0.585089) (7.50996 -0.0487036 0.595184 0.450894 -0.586047 0.310849))
;;
;; (Jacobi a 1e-12)
(setq t1 (getvar "millisecs")
d (gc:determ mat)
q ([i] (length (car mat)))
i 0 ;_Set the maximum number of iterations to prevent stack overflow
)
(while (and (not (equal d ([tr] mat) eps)) (< i 1000))
(setq b (Givens mat Q)
mat (car b)
q (cadr b))
(setq i (1+ i)))
(princ (strcat "\nJacobi iterations is : "
(rtos i 2)
" times ."))
(princ (strcat "\nTotel time of Jacobi solving : "
(rtos (- (getvar "millisecs") t1) 2 0)
" m.s."))
(setq i -1)
(mapcar (function (lambda (a b)
(cons a b)))
(mapcar (function (lambda (r1)
(nth (setq i (1+ i)) r1)))
mat)
q)
)
Great job!thank you,chlh_jd!You are always welcome , and thank you for help ! :-)
;;;----------------------------------Householder变换------------------------------;;;
(defun householder (q mat / i m n d a b u hk)
;; Use Householder transformation to Get QR matrix for QR-decomposition
;; (householder mat)
;; (setq mat '((2. 3.) (1. 2.) (5. 7.)))
;; (setq mat '((1. 2. -3.) (1. 3. 10.) (1. -1. 6.)))
(setq i 0
n (length (car mat))
)
(while (< i n)
(setq d (mapcar (function (lambda (x) (nth i x))) mat))
(repeat i
(setq d (cdr d)))
(setq a (* (sign (car d))
(sqrt (apply (function +)
(mapcar (function (lambda (x)
(* x x)))
d))))
b (* a (+ a (car d)))
d (cons (+ (car d) a)
(cdr d)))
(repeat i
(setq d (cons 0. d)))
(setq u (vtxv d d))
(if (not (zerop b))
(setq u ([k*] (1/ b) u)))
(setq hk ([-] imat u)
q ([*] q hk)
mat ([*] ([*] hk mat) ([trp] hk))
i (1+ i)))
(list q mat))
;;;------------------------------------QR分解迭代---------------------------------;;;
;;
(defun QR->EV (mat eps / imat q d i b t1)
;; QR-decomposition
;; by GSLS(SS) 2013-1-2
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.)))
;; (QR->EV a 5e-11)
(setq t1 (getvar "millisecs")
d (gc:determ mat)
imat ([i] (length (car mat)))
q imat
i 0) ;_Set the maximum number of iterations to prevent stack overflow
(while (and (not (equal d ([tr] mat) eps)) (< i 1000))
(setq b (householder q mat)
mat (cadr b)
q (car b)
i (1+ i)))
(princ (strcat "\nQR-decomposition of iterations is : "
(rtos i 2)
" times ."))
(princ (strcat "\nTotel time of QR-decomposition : "
(rtos (- (getvar "millisecs") t1) 2 0)
" m.s."))
(setq i -1)
(mapcar (function (lambda (a b)
(cons a b)))
(mapcar (function (lambda (r1)
(nth (setq i (1+ i)) r1)))
mat)
([trp] q))
)
...[*]...
example
.....
(setq hk ([-] imat u)
q ([*] q hk)
mat ([*] ([*] hk mat)
......
;;;-------------------------------------------------------------------------------;;;
;;; 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
)
)
Maybe I have to firstly normalize matrix, like this user wanted to explain :
https://mymathforum.com/threads/how-do-i-normalize-a-matrix.18218/#post-159640
And then multiply matrix ~ 40 times to get all rows the same (row=eigenvector)... Of course eigenvalue should have it's purpose, so I suppose that matrix for multiplication should be : mat-mat[E] - just guessing
Something like this was shown at the bottom of page I already mentioned :
https://www.intmath.com/matrices-determinants/8-applications-eigenvalues-eigenvectors.php
Now the question : How to normalize matrix with elements - complex numbers?