Author Topic: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?  (Read 18808 times)

0 Members and 1 Guest are viewing this topic.

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #15 on: December 28, 2012, 09:00:15 AM »
I think Master ElpanovEvgeniy and Gile have all solved this problem .  :-)
I sincerely look forward to more guidance .  :-)

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #16 on: December 31, 2012, 03:08:43 AM »
New main function
Code: [Select]
(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))) ;|)|;
  )

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #17 on: January 02, 2013, 03:00:34 AM »
Last Version
First Use QR-decomposition to get all EigenValues , Then Use Gaussian-elimination to get all EigenVectors .

Thanks Robert and ElpanovEvgeniy for give me so much help .  :-)
Thanks ElpanovEvgeniy's Gaussian-elimination function , it's really a pretty codes .

Any suggestions are welcome , and please post them to chlh_jd@126.com or 56454300@qq.com

main functions see following codes , to see all functions please download the add file .

Code: [Select]
;;;-------------------------------------------------------------------------------;;;
;;;                         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
;;;


Roben25

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #18 on: January 02, 2013, 11:49:37 AM »
Nice to see that you made it and thanks for share :)

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #19 on: January 03, 2013, 12:19:06 PM »
Nice to see that you made it and thanks for share :)
Thanks Robert ! You are always welcome.  :-)

I add another method following .

Power & Contraction method  to solve EigenValue
Code: [Select]
;;;---------------------------------------幂法-----------------------------------;;;
(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 transformation method to solve EigenValues (only  Eigenvalues . In fact, if combined  with Jacobi method , it can solve Eigenvalues and Eigenvectors at the same time .) Givens method is higher than QR method in the same torrance .
Code: [Select]
;;--------------------------------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)
  )

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #20 on: January 04, 2013, 06:43:35 AM »
Use Givens transformation and Jacobi method .
Code: [Select]
(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)
  )

highflyingbird

  • Bull Frog
  • Posts: 415
  • Later equals never.
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #21 on: January 04, 2013, 07:43:04 AM »
Great job!thank you,chlh_jd!
I am a bilingualist,Chinese and Chinglish.

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #22 on: January 05, 2013, 04:42:55 AM »
Great job!thank you,chlh_jd!
You are always welcome , and thank you for help ! :-)

Change some codes of QR-decomposition , now it can solve the Eigenvalues and Eigenvectors at the same time .
However , QR method is lower efficient than Jacobi method ; To solving the stable Eigenvectors , QR method must Improve double accuracy than Jacobi .
Code: [Select]
;;;----------------------------------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))
  )

kosak

  • Mosquito
  • Posts: 8
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #23 on: March 27, 2013, 11:59:05 AM »
Hi chlh_jd,

i am trying to run the program  QR-Gauss-Evs.LSP but the function
Code: [Select]
...[*]...
example
.....
(setq hk  ([-] imat u)
  q   ([*] q hk)
  mat ([*] ([*] hk mat)
......

 is not included in the program code ?? 
Could you add/attach this function and check that the code functions properly ?

    Regards
    Konstantin
     

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #24 on: March 29, 2013, 09:17:46 AM »
@ kosak
  Please download Used functions from Reply 14#

kosak

  • Mosquito
  • Posts: 8
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #25 on: April 03, 2013, 07:46:58 AM »
Ok chlh_jd,

many thanks for your replay !

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #26 on: December 09, 2019, 03:03:18 AM »
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...

Code: [Select]
;;;-------------------------------------------------------------------------------;;;
;;;                         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.
« Last Edit: December 22, 2019, 11:08:24 AM by ribarm »
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #27 on: December 18, 2019, 06:33:55 PM »
Maybe now with my latest addition of matrix determinant and for complex numbers, someone will tackle this problem more deeply and find both real and complex eigenvalues and their eigenvectors...

Here is the link :
http://www.theswamp.org/index.php?topic=45666.msg597622#msg597622
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #28 on: December 22, 2019, 11:14:19 AM »
I did find eigenvalues - it's so slooow, but should work and for complex eigenvalue numbers... Now I don't have so much math knowledge, but I suppose that someone will be able to derive eigenvectors according to eigenvalues... So here is it (don't laugh), I believe this could be more elegant and faster - www solvers find this in a second, but my CAD need minute or more...

Code - Auto/Visual Lisp: [Select]
  1. (defun [E]-LST ( mat / detm detmc eqmember [I] [k*] [-] f1 f2 f2c process eigr imat matc imatc breakdown eps mi ma st si er res resc )
  2.  
  3. ;;; Command: ([e] '((1.0 2.0 0.5 3.0) (4.0 1.0 0.6 -2.7) (-1.2 2.3 4.1 -3.2) (-0.2 -0.6 1.4 -0.9)))
  4. ;;; ((-3.10709 0.0) (1.81132 -1.98745) (1.81132 1.98745) (4.68446 0.0))
  5. ;;; Command: ([e] '((0.5 -1.2 -3.4 0.7) (0.24 3.1 -2.6 -8.3) (-4.2 5.3 6.2 1.8) (-3.1 -7.2 -6.3 -2.4)))
  6. ;;; ((-4.01463 2.2902) (-4.01463 -2.2902) (5.95272 0.0) (9.47654 0.0))
  7.  
  8.   (defun detm ( m / d )
  9.     ;; Matrix Determinant  -  ElpanovEvgeniy
  10.     ;; Last edit 2013.11.13
  11.     ;; Args: m - nxn matrix
  12.     ;; (detm '((0 1) (1 0)))
  13.     ;; (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)))
  14.     (cond ((null m) 1)
  15.           ((and (zerop (caar m))
  16.                 (setq d (car (vl-member-if-not (function (lambda ( a ) (zerop (car a)))) (cdr m))))
  17.            )
  18.            (detm (cons (mapcar (function +) (car m) d) (cdr m)))
  19.           )
  20.           ((zerop (caar m)) 0)
  21.           ((* (caar m)
  22.               (detm (mapcar (function (lambda ( a / d )
  23.                                         (setq d (/ (car a) (float (caar m))))
  24.                                         (mapcar (function (lambda (b c) (- b (* c d)))) (cdr a) (cdar m))
  25.                                       )
  26.                             )
  27.                             (cdr m)
  28.                     )
  29.               )
  30.            )
  31.           )
  32.     )
  33.   )
  34.  
  35.   ;;;***********************************************************************************;;;
  36.   ;;; (detmc m) function calculates determinant of square martix                        ;;;
  37.   ;;; this version includes complex numbers                                             ;;;
  38.   ;;; Marko Ribar, d.i.a.                                                               ;;;
  39.   ;;; Args: m - nxn matrix of complex numbers                                           ;;;
  40.   ;;; (detmc '(((0 0) (1 0)) ((1 0) (0 0))))                                            ;;;
  41.   ;;; (detmc '(((1.4 0) (2.1 0) (5.4 0) (6.5 0)) ((4.1 0) (9.3 0) (4.5 0) (8.5 0)) ((1.2 0) (4.1 0) (6.2 0) (7.5 0)) ((4.7 0) (8.5 0) (9.3 0) (0.1 0))))                                                                                 ;;;
  42.   ;;;***********************************************************************************;;;
  43.  
  44.   ;;; Command: (detmc '(((1.2 1.0) (2.4 1.3) (0.5 0.6)) ((1.3 1.25) (1.5 2.0) (0.8 -1.4)) ((-2.3 -1.4) (0.6 -0.6) (-0.3 -0.3))))
  45.   ;;; (-15.3785 7.4325)
  46.  
  47.   ;;; det = (1.2 + 1.0i)(1.5 + 2.0i)(-0.3 - 0.3i)+(2.4 + 1.3i)(0.8 - 1.4i)(-2.3 - 1.4i)+(0.5 + 0.6i)(1.3 + 1.25i)(0.6 - 0.6i)-(-2.3 - 1.4i)(1.5 + 2.0i)(0.5 + 0.6i)-(0.6 - 0.6i)(0.8 - 1.4i)(1.2 + 1.0i)-(-0.3 - 0.3i)(1.3 + 1.25i)(2.4 + 1.3i)
  48.   ;;; det = (1.8 + 2.4i + 1.5i - 2)(-0.3 - 0.3i)+(1.92 - 3.36i + 1.04i + 1.82)(-2.3 - 1.4i)+(0.65 + 0.625i + 0.78i - 0.75)(0.6 - 0.6i)-(-3.45 - 4.6i - 2.1i + 2.8)(0.5 + 0.6i)-(0.48 - 0.84i - 0.48i - 0.84)(1.2 + 1.0i)-(-0.39 - 0.375i - 0.39i + 0.375)(2.4 + 1.3i)
  49.   ;;; det = (-0.2 + 3.9i)(-0.3 - 0.3i)+(3.74 - 2.32i)(-2.3 - 1.4i)+(-0.1 + 1.405i)(0.6 - 0.6i)-(-0.65 - 6.7i)(0.5 + 0.6i)-(-0.36 - 1.32i)(1.2 + 1.0i)-(-0.015 - 0.765i)(2.4 + 1.3i)
  50.   ;;; det = 0.06 + 0.06i - 1.17i + 1.17 - 8.602 - 5.236i + 5.336i - 3.248 - 0.06 + 0.06i + 0.843i + 0.843 + 0.325 + 0.39i + 3.35i - 4.02 + 0.432 + 0.36i + 1.584i - 1.32 + 0.036 + 0.0195i + 1.836i - 0.9945
  51.   ;;; det = -15.3785 + 7.4325i
  52.  
  53.   (defun detmc ( m / d cxc det )
  54.  
  55.     (defun d (k n / z)
  56.       (setq k (cdr k))
  57.       (setq k (apply (function mapcar) (cons (function list) k)))
  58.       (setq z -1)
  59.       (while (<= (setq z (1+ z)) (length k))
  60.         (if (= z n)
  61.           (setq k (cdr k))
  62.           (setq k (reverse (cons (car k) (reverse (cdr k)))))
  63.         )
  64.       )
  65.       (setq k (apply (function mapcar) (cons (function list) k)))
  66.       (if (= (length k) 1) (caar k) k)
  67.     )
  68.  
  69.     (defun cxc ( c1 c2 / a1 a2 b1 b2 r i )
  70.       (setq a1 (car c1) b1 (cadr c1))
  71.       (setq a2 (car c2) b2 (cadr c2))
  72.       (setq r (- (* a1 a2) (* b1 b2)))
  73.       (setq i (+ (* a1 b2) (* a2 b1)))
  74.       (list r i)
  75.     )
  76.  
  77.     (defun det ( q / i j r )
  78.       (if (not (vl-every (function numberp) q))
  79.         (progn
  80.           (setq i -1)
  81.           (setq j -1)
  82.           (setq r (list 0 0))
  83.           (foreach e (car q)
  84.             (setq i (1+ i))
  85.             (setq j (* j (- 1)))
  86.             (setq r (mapcar (function +) r (cxc (mapcar (function *) (list j j) e) (det (d q i)))))
  87.           )
  88.           r
  89.         )
  90.         q
  91.       )
  92.     )
  93.  
  94.     (det m)
  95.   )
  96.  
  97.   (defun eqmember ( p l eps )
  98.     (cond ((not l) nil)
  99.           ((if (numberp eps)
  100.              (equal p (car l) eps)
  101.              (equal p (car l))
  102.            )
  103.            l
  104.           )
  105.           ((eqmember p (cdr l) eps))
  106.     )
  107.   )
  108.  
  109.   (defun [I] ( d / i r n m )
  110.     (setq i d)
  111.     (while (<= 0 (setq i (1- i)))
  112.       (setq n d
  113.             r nil
  114.       )
  115.       (while (<= 0 (setq n (1- n)))
  116.         (setq r (cons (if (= i n) 1.0 0.0) r))
  117.       )
  118.       (setq m (cons r m))
  119.     )
  120.   )
  121.  
  122.   (defun [k*] ( k a )
  123.     (mapcar
  124.       (function
  125.         (lambda ( x )
  126.           (mapcar
  127.             (function
  128.               (lambda ( y )
  129.                 (cond
  130.                   ( (listp y)
  131.                     (mapcar (function *) k y)
  132.                   )
  133.                   ( (numberp y)
  134.                     (* k y)
  135.                   )
  136.                 )
  137.               )
  138.             ) x
  139.           )
  140.         )
  141.       ) a
  142.     )
  143.   )
  144.  
  145.   (defun [-] ( a b )
  146.     (mapcar
  147.       (function
  148.         (lambda ( x y )
  149.           (cond
  150.             ( (and (listp (car x)) (listp (car y)))
  151.               (mapcar (function (lambda (xx yy) (mapcar (function -) xx yy))) x y)
  152.             )
  153.             ( (and (numberp (car x)) (numberp (car y)))
  154.               (mapcar (function -) x y)
  155.             )
  156.           )
  157.         )
  158.       ) a b
  159.     )
  160.   )
  161.  
  162.   ;; To get all Eigenvalues
  163.   (defun f1 ( mat / i v a r rr )
  164.     ;; Eigenvalue interval analysis
  165.     (setq i 0)
  166.     (repeat (length (car mat))
  167.       (setq v (nth i mat)
  168.             a (apply (function +) (mapcar (function abs) v))
  169.             r (cons (list (- (apply (function min) v) a) (+ (apply (function max) v) a)) r)
  170.             i (1+ i)
  171.       )
  172.     )
  173.     (setq r (apply (function append) r))
  174.     (setq rr (list (cons (apply (function min) r) (apply (function max) r))))
  175.     rr
  176.   )
  177.   ;; Iteration
  178.   (defun f2 ( mi si / chk ei ri mil r11 r21 r12 r22 rr )
  179.     (setq chk mi)
  180.     (setq ei (list mi)
  181.           mil mi
  182.     )
  183.     (repeat (fix (/ breakdown 2))
  184.       (setq mi (+ mi si)
  185.             ei (cons mi ei)
  186.       )
  187.     )
  188.     (setq ei (reverse ei))
  189.     (repeat (fix (/ breakdown 2))
  190.       (setq mil (- mil si)
  191.             ei  (cons mil ei)
  192.       )
  193.     )
  194.     (setq ri
  195.       (mapcar
  196.         (function
  197.           (lambda (e)
  198.             (cons e (abs (detm ([-] ([k*] e imat) mat))))
  199.           )
  200.         ) ei
  201.       )
  202.     )
  203.     (foreach x ri
  204.       (if (> (cdr (nth (1+ (vl-position x ri)) ri)) (cdr x))
  205.         (setq r11 (cons x r11))
  206.       )
  207.     )
  208.     (foreach x r11
  209.       (if (> (cdr (nth (1+ (vl-position x r11)) r11)) (cdr x))
  210.         (setq r21 (cons x r21))
  211.       )
  212.     )
  213.     (foreach x (setq ri (reverse ri))
  214.       (if (> (cdr (nth (1+ (vl-position x ri)) ri)) (cdr x))
  215.         (setq r12 (cons x r12))
  216.       )
  217.     )
  218.     (foreach x r12
  219.       (if (> (cdr (nth (1+ (vl-position x r12)) r12)) (cdr x))
  220.         (setq r22 (cons x r22))
  221.       )
  222.     )
  223.     (setq rr (append r21 r22))
  224.     (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)))
  225.       (setq rr (cons (last r11) rr))
  226.     )
  227.     (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)))
  228.       (setq rr (cons (last r12) rr))
  229.     )
  230.     (if (null rr)
  231.       (setq rr (cons (car (vl-sort ri (function (lambda (x y) (< (cdr x) (cdr y)))))) rr))
  232.     )
  233.     (if (equal chk (caar rr) 1e-5)
  234.       (vl-sort rr (function (lambda (x y) (< (cdr x) (cdr y)))))
  235.       rr
  236.     )
  237.   )
  238.   ;; Iteration
  239.   (defun f2c ( mi1 mi2 si / ei1 ei2 eii ri mil sort l )
  240.     (setq ei1 (list mi1)
  241.           mil mi1
  242.     )
  243.     (repeat (fix (/ breakdown 2))
  244.       (setq mi1 (+ mi1 si)
  245.             ei1 (cons mi1 ei1)
  246.       )
  247.     )
  248.     (setq ei1 (reverse ei1))
  249.     (repeat (fix (/ breakdown 2))
  250.       (setq mil (- mil si)
  251.             ei1  (cons mil ei1)
  252.       )
  253.     )
  254.     (setq ei2 (list mi2)
  255.           mil mi2
  256.     )
  257.     (repeat (fix (/ breakdown 2))
  258.       (setq mi2 (+ mi2 si)
  259.             ei2 (cons mi2 ei2)
  260.       )
  261.     )
  262.     (setq ei2 (reverse ei2))
  263.     (repeat (fix (/ breakdown 2))
  264.       (setq mil (- mil si)
  265.             ei2  (cons mil ei2)
  266.       )
  267.     )
  268.     (foreach x ei1
  269.       (foreach y ei2
  270.         (setq eii (cons (list x y) eii))
  271.       )
  272.     )
  273.     (setq ri
  274.       (mapcar
  275.         (function
  276.           (lambda ( e )
  277.             (list e (apply (function +) (mapcar (function abs) (detmc ([-] ([k*] e imatc) matc)))))
  278.           )
  279.         ) eii
  280.       )
  281.     )
  282.     (setq l
  283.       (vl-remove nil
  284.         (append
  285.           (if (vl-remove-if-not (function (lambda ( x ) (> (cadar x) eps))) ri)
  286.             (vl-remove-if (function (lambda ( x ) (vl-member-if (function (lambda ( y ) (equal (caar x) y 0.5))) res))) (vl-remove-if (function (lambda ( x ) (> (vl-position x sort) (1- (- (length (car mat)) (length res)))))) (setq sort (vl-sort (vl-remove-if-not (function (lambda ( x ) (> (cadar x) eps))) ri) (function (lambda ( a b ) (< (cadr a) (cadr b))))))))
  287.           )
  288.           (if (vl-remove-if-not (function (lambda ( x ) (< (cadar x) (- eps)))) ri)
  289.             (vl-remove-if (function (lambda ( x ) (vl-member-if (function (lambda ( y ) (equal (caar x) y 0.5))) res))) (vl-remove-if (function (lambda ( x ) (> (vl-position x sort) (1- (- (length (car mat)) (length res)))))) (setq sort (vl-sort (vl-remove-if-not (function (lambda ( x ) (< (cadar x) (- eps)))) ri) (function (lambda ( a b ) (< (cadr a) (cadr b))))))))
  290.           )
  291.         )
  292.       )
  293.     )
  294.     (process l)
  295.   )
  296.  
  297.   (defun process ( l / unique )
  298.  
  299.     (defun unique ( l )
  300.       (if l (cons (car l) (unique (vl-remove-if (function (lambda ( x ) (equal (car x) (caar l) 0.5))) l))))
  301.     )
  302.  
  303.     (setq l (vl-remove-if-not (function (lambda ( x ) (vl-some (function (lambda ( y ) (and (equal (caar x) (caar y) 0.5) (equal (abs (cadar x)) (abs (cadar y)) 0.5)))) (vl-remove x l)))) l))
  304.     (setq l (vl-sort l (function (lambda ( a b ) (< (cadr a) (cadr b))))))
  305.     (unique l)
  306.   )
  307.  
  308.   ;; MAIN Eigenvalue ;;
  309.  
  310.   (setq eigr (f1 mat)
  311.         imat ([I] (length (car mat)))
  312.         breakdown 300 ;_here is a main arg .
  313.         eps 5e-3 ;_Tolerance value
  314.   )
  315.   (foreach a eigr
  316.     (setq mi (car a)
  317.           ma (cdr a)
  318.           st (- ma mi)
  319.           mi (/ (+ mi ma) 2.)
  320.           si (/ st breakdown)
  321.           er (f2 mi si)
  322.     )
  323.     (while (and (vl-some (function (lambda ( x ) (> (cdr x) eps))) er) (> si 1e-8))
  324.       (setq si (* 1.25 (/ si breakdown)))
  325.       (setq er
  326.         (mapcar
  327.           (function
  328.             (lambda ( x )
  329.               (car (f2 (car x) si))
  330.             )
  331.           ) er
  332.         )
  333.       )
  334.     )
  335.     (foreach a (vl-sort er (function (lambda ( x y ) (< (cdr x) (cdr y)))))
  336.       (if (and (<= (cdr a) eps) (not (eqmember (car a) res eps)))
  337.         (setq res (cons (car a) res))
  338.       )
  339.     )
  340.   )
  341.   (if (< (length res) (length (car mat)))
  342.     (progn
  343.       (setq imatc (mapcar (function (lambda ( x ) (mapcar (function (lambda ( y ) (list y y))) x))) imat)
  344.             matc (mapcar (function (lambda ( x ) (mapcar (function (lambda ( y ) (list y 0.0))) x))) mat)
  345.             breakdown 50 ;_here is a main arg .
  346.       )
  347.       (foreach a eigr
  348.         (setq mi (car a)
  349.               ma (cdr a)
  350.               st (- ma mi)
  351.               mi (/ (+ mi ma) 2.)
  352.               si (/ st breakdown)
  353.               er (f2c mi mi si)
  354.         )
  355.         (while (and (vl-some (function (lambda ( x ) (> (cadr x) eps))) er) (> si 1e-8))
  356.           (setq si (* 1.25 (/ si breakdown)))
  357.           (setq er
  358.             (apply
  359.               (function append)
  360.               (mapcar
  361.                 (function
  362.                   (lambda ( x )
  363.                     (f2c (caar x) (cadar x) si)
  364.                   )
  365.                 ) er
  366.               )
  367.             )
  368.           )
  369.           (setq er (process er))
  370.         )
  371.         (foreach a er
  372.           (if (and (<= (cadr a) eps) (not (eqmember (car a) resc eps)))
  373.             (setq resc (cons (car a) resc))
  374.           )
  375.         )
  376.       )
  377.     )
  378.   )
  379.   (vl-sort (append (mapcar (function (lambda ( x ) (list x 0.0))) res) resc) (function (lambda ( a b ) (< (car a) (car b)))))
  380. )
  381.  

Regards, M.R.
« Last Edit: August 17, 2020, 08:31:46 PM by ribarm »
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #29 on: December 24, 2019, 07:22:18 AM »
I tried something, but don't really know how to continue...

Code - Auto/Visual Lisp: [Select]
  1. (defun [V] ([E] mat / trp mxvc mxmc cxc c-c c+cm [E]mat v0 mat[E] mat-mat[E])
  2.  
  3.   (defun trp ( m )
  4.   )
  5.  
  6.   (defun mxvc ( m v )
  7.     (mapcar (function (lambda ( r ) (c+cm (mapcar (function cxc) r v)))) m)
  8.   )
  9.  
  10.   (defun mxmc ( m n )
  11.     ((lambda ( a ) (mapcar (function (lambda ( r ) (mxvc a r))) m)) (trp n))
  12.   )
  13.  
  14.   (defun cxc ( c1 c2 / a1 a2 b1 b2 r i )
  15.     (setq a1 (car c1) b1 (cadr c1))
  16.     (setq a2 (car c2) b2 (cadr c2))
  17.     (setq r (- (* a1 a2) (* b1 b2)))
  18.     (setq i (+ (* a1 b2) (* a2 b1)))
  19.     (list r i)
  20.   )
  21.  
  22.   (defun c-c ( c1 c2 )
  23.     (mapcar (function -) c1 c2)
  24.   )
  25.  
  26.   (defun c+cm ( clst )
  27.     (list (apply (function +) (mapcar (function car) clst)) (apply (function +) (mapcar (function cadr) clst)))
  28.   )
  29.  
  30.   (defun [E]mat ( [E] mat / k rown n matn )
  31.     (setq k -1)
  32.     (foreach row mat
  33.       (setq k (1+ k))
  34.       (setq rown (mapcar (function (lambda ( x ) (if (null n) (setq n 0) (setq n (1+ n))) (if (= n k) [E] (list 0.0 0.0)))) row))
  35.       (setq n nil)
  36.       (setq matn (cons rown matn))
  37.     )
  38.     (reverse matn)
  39.   )
  40.  
  41.   (defun v0 ( mat )
  42.     (mapcar (function (lambda ( x ) (list 0.0 0.0))) (car mat))
  43.   )
  44.  
  45.   (setq mat (mapcar (function (lambda ( x ) (mapcar (function (lambda ( y ) (list y 0.0))) x))) mat))
  46.   (setq mat[E] ([E]mat [E] mat))
  47.   (setq mat-mat[E] (mapcar (function (lambda ( a b ) (mapcar (function (lambda ( c d ) (c-c c d))) a b))) mat mat[E]))
  48.   ;;; (equal (mapcar (function c-c) (mxvc mat v) (mxvc mat[E] v)) (v0 mat) 1e-8) ;;; T ;;; v=? ;;; v /= ((0.0 0.0) (0.0 0.0) ... (0.0 0.0))
  49.   ;;; (equal (mxvc mat-mat[E] v) (v0 mat) 1e-8) ;;; T ;;; v=? ;;; v /= ((0.0 0.0) (0.0 0.0) ... (0.0 0.0))
  50.   ;;; https://www.intmath.com/matrices-determinants/8-applications-eigenvalues-eigenvectors.php
  51.   ;;; There is some example of using multiplication of matrices, but from this point I don't know how to continue to be able to get eigenvector...
  52. )
  53.  
« Last Edit: December 24, 2019, 03:35:52 PM by ribarm »
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube