Author Topic: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?  (Read 9574 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]

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

  • Water Moccasin
  • Posts: 2226
  • 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

  • Water Moccasin
  • Posts: 2226
  • 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

  • Water Moccasin
  • Posts: 2226
  • 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] (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: December 23, 2019, 10:29:10 AM by ribarm »
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Water Moccasin
  • Posts: 2226
  • 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

ribarm

  • Water Moccasin
  • Posts: 2226
  • Marko Ribar, architect
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #30 on: December 26, 2019, 05:25:12 AM »
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?
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Water Moccasin
  • Posts: 2226
  • Marko Ribar, architect
Re: How to Get Eigenvalues ​​and Eigenvectors of the matrix ?
« Reply #31 on: December 26, 2019, 05:53:01 AM »

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?


I just realized that this path is wrong, but I searched for answers, so nothing wrong...
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube