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

0 Members and 1 Guest are viewing this topic.

chlh_jd

  • Guest
How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« on: December 26, 2012, 04:16:03 AM »
Hi all , as the title .
How to get Eigenvalues and Eigenvectors of the matrix by Vlisp ?

Here's Power method to get one Eigenvalue and it's Eigenvector , but how to get all ?
Code: [Select]
(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
    )
  )
)
« Last Edit: December 26, 2012, 04:22:38 AM by chlh_jd »

Roben25

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #1 on: December 26, 2012, 09:24:29 PM »
Hi, power method gets only the dominant eigenpair, you have to use other method like SVD, QR or Jacobi. First of all do you have symmetric or non symmetric matrix?

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #2 on: December 27, 2012, 07:51:44 AM »
Thanks Roben25 !
As I knew :
      Power method just can get the Real_Matrix's One main Eigenvalue and it's Eigenvector .
      QR method just can get all Eigenvalues .
      Svd method , I just know it can use Matlab solove the problem ?
Can you  explain more about this ?

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #3 on: December 27, 2012, 08:10:43 AM »
If has got all Eigenvalues , how get their Eigenvectors ?

Use Interation get all Eigenvalues .
Code: [Select]
;;
(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))

Roben25

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #4 on: December 27, 2012, 09:29:44 AM »
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

highflyingbird

  • Bull Frog
  • Posts: 415
  • Later equals never.
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #5 on: December 27, 2012, 10:10:36 AM »
If has got all Eigenvalues , how get their Eigenvectors ?

If just consider 3 dimension ,then it will get a better way:
Code - Auto/Visual Lisp: [Select]
  1. ;;;----------------------------------------------------;
  2. ;;;EigenValues in n^3                                  ;
  3. ;;;Input: mat -- a 3x3 matrix                          ;
  4. ;;;Output: the EigenValues of this matrix              ;
  5. ;;;----------------------------------------------------;
  6. (defun MAT:EigenValues (mat / a b c d e f g h i m n k z)
  7.   (setq s (apply 'append mat))
  8.   (mapcar 'set '(a b c d e f g h i) s)
  9.   (setq m (- (+ a e i)))
  10.   (setq n (- (+ (* a e) (* e i) (* i a)) (* b d) (* c g) (* f h)))
  11.   (setq z (apply 'MAT:Det3 s))
  12.   (Math:Cubic_Equation 1 m n (- z))
  13. )
  14.  

Other code see the attachment.
I am a bilingualist,Chinese and Chinglish.

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #6 on: December 27, 2012, 10:19:44 AM »
If just consider 3 dimension ,then it will get a better way:
Code - Auto/Visual Lisp: [Select]
  1. ;;;----------------------------------------------------;
  2. ;;;EigenValues in n^3                                  ;
  3. ;;;Input: mat -- a 3x3 matrix                          ;
  4. ;;;Output: the EigenValues of this matrix              ;
  5. ;;;----------------------------------------------------;
  6. (defun MAT:EigenValues (mat / a b c d e f g h i m n k z)
  7.   (setq s (apply 'append mat))
  8.   (mapcar 'set '(a b c d e f g h i) s)
  9.   (setq m (- (+ a e i)))
  10.   (setq n (- (+ (* a e) (* e i) (* i a)) (* b d) (* c g) (* f h)))
  11.   (setq z (apply 'MAT:Det3 s))
  12.   (Math:Cubic_Equation 1 m n (- z))
  13. )
  14.  
Other code see the attachment.
Thanks highflyingbird !  :-)

highflyingbird

  • Bull Frog
  • Posts: 415
  • Later equals never.
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #7 on: December 27, 2012, 10:23:45 AM »
I think the next step is to get the Eigenvectors.
Maybe use Gaussian elimination。
I am a bilingualist,Chinese and Chinglish.

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #8 on: December 27, 2012, 10:58:13 AM »
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

Thanks Robert , thank you very much !  :-)

I think the next step is to get the Eigenvectors.
Maybe use Gaussian elimination。
Thanks , I 'll search and study 'use Gaussian elimination ' . :-)

highflyingbird

  • Bull Frog
  • Posts: 415
  • Later equals never.
« Last Edit: December 27, 2012, 11:25:27 AM by highflyingbird »
I am a bilingualist,Chinese and Chinglish.


chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #11 on: December 27, 2012, 02:46:05 PM »
Newversion , can get more Eigenvalues and it's Eigenvectors , but has bug yet .
Code: [Select]
(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))
  )
)

Roben25

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #12 on: December 27, 2012, 03:33:43 PM »
Hi, here some code i have made last year, its compute eigenvector of 3x3 matrix

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

Result:
((-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

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #13 on: December 28, 2012, 08:00:45 AM »
Hi, here some code i have made last year, its compute eigenvector of 3x3 matrix

Code: [Select]
(setq mat
   '(
        (2.0 3.0 8.0)
        (1.0 0.0 5.0)
        (-1.0 5.0 0.0)
    )
)
Result:
((-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


Thanks Robort for share it ! :-)

MATLAB is not installed on my computer now , so I can't check my result .

The matrix your provided , test by my routine take the following result :
Code: [Select]
(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)
Ei 1.9258 -->EV (-6.19258 -1.0 1.0) --> (0.9749  0.15743 -0.15743)
« Last Edit: December 28, 2012, 09:09:29 AM by chlh_jd »

chlh_jd

  • Guest
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #14 on: December 28, 2012, 08:55:20 AM »
Fix some bug in Eigenvalues function , see the upload Lsp file .
Code: [Select]
(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]