明经CAD社区

 找回密码
 注册

QQ登录

只需一步,快速开始

搜索
查看: 11801|回复: 16

实对称矩阵的特征值与特征向量求解

    [复制链接]
发表于 2013-1-4 00:56:09 | 显示全部楼层 |阅读模式
本帖最后由 chlh_jd 于 2013-1-4 01:12 编辑

求非奇异实数矩阵A特征值与特征向量常用的求解方法有: 幂法,反幂法,QR迭代法, Givens变换收缩法,Jocobi法,同步求解法。

幂法用来求第一个特征值和特征向量,后面将给出第一次计算范例(详细步骤可在百度文库找到)。幂法收敛速度通常比较慢,如果提高计算速度就要损失精度,当前后两次迭代得到的特征值较为接近时,收敛很慢,因此还衍生出了加入瑞利(Rayleigh)商加速法、原位平移法等。幂法+收缩法就可以求全部特征值和特征向量,主要步骤为:
1.对矩阵A第一次求出特征值λ1和第1个特征向量x1;
2. 然后用0代替λ1所在行,将矩阵降阶再次使用幂法,就可以得到第2个特征值λ2与特征向量x2 ;
3.依次类推,最终得到所有的特征值λ
4. 对并列矩阵[ (λiI-A)^t | I ] 采用高斯消元法转为上三角矩阵,左阵0向量或相邻向量线性相关的向量对应的右阵向量即为特征向量,依次代入所有特征值求出所有特征向量。

  1. ;;;---------------------------------------幂法-----------------------------------;;;
  2. (defun PowerEigen  (mat eps / v0 f0 vk f1)
  3.   ;; 求矩阵第一特征值与特征向量求解
  4.   ;; Get First EigenValue and it's EigenVector
  5.   ;; mat 为方阵 ; eps 为精度、该精度不宜太小否收敛很慢, 宜为1e-6~1e-8
  6.   ;; by GSLS(SS) 2012-12-24
  7.   ;; (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)
  8.   ;; (setq a '((7. 3. -2.) (3. 4. 1.) (-2. -1. 3.))) --> (9.27492 (1.0 0.493399 -0.39736))
  9.   ;; (setq a '((4. 1. 0.) (1. 0. -1.) (1. 1. -4.))) --> (4.20303 (1.0 0.20303 0.146657))
  10.   ;; (setq a '((1. 1. 0.5) (1. 1. 0.25) (0.5 0.25 2.))) --> (2.53652587 (0.74822116 0.64966115 1.0))
  11.   ;; (PowerEigen a 1e-9)
  12.   (setq v0 (mapcar (function (lambda (x) 1.0)) (car mat))
  13. f0 1.0
  14. vk (mxv mat v0)
  15. f1 (apply (function max) (mapcar (function abs) vk))
  16. vk (pt* vk (1/ f1)))
  17.   (while
  18.     (not
  19.       (and
  20. (<= (abs (- f1 f0)) eps)
  21. (<=
  22.    ((lambda (x)
  23.       (sqrt (apply (function +) (mapcar (function *) x x))))
  24.      (pt- vk v0))
  25.    eps)))
  26.      (setq v0 vk
  27.     f0 f1
  28.     vk (mxv mat v0)
  29.     f1 (apply (function max) (mapcar (function abs) vk))
  30.     vk (pt* vk (1/ f1))))
  31.   (cons f1 vk))
  32. ;;
  33. (defun Power->EV (mat eps / m  v  v1 eig res)
  34.   ;; Get all EigenValues
  35.   ;; by GSLS(SS) 2013-1-3
  36.   ;; (PowerEigenValues a)
  37.   (while (cadr mat)
  38.     (setq eig (powereigen mat eps)
  39.    res (cons eig res))
  40.     (if res
  41.       (setq v (car mat)
  42.      v1 (cdar res)
  43.      v1 (pt* v1 (1/ (car v1)))
  44.      m (vtxv v1 v)
  45.      mat ([-] mat m)
  46.      mat (mapcar (function cdr) (cdr mat))
  47.       ))     
  48.     )
  49.   (reverse (mapcar (function car) (cons (car mat) res)))
  50.   )

反幂法用来求最小特征值和特征向量,其迭代过程与幂法雷同。


QR迭代法20世纪60年代产生的一种求中小规模的矩阵的特征值与特征向量的方法,先求得准确的矩阵模,将该模与叠代中间矩阵的主迹(脚标相同的元素积)做比较,如果在容差范围内,则停止迭代,否则进行下一步迭代。QR法的主要步骤:
R——矩阵A的上三角矩阵,A=QR
AHouseholder(镜像变换,由著名的数值分析专家Householder1958年为讨论矩阵特征值问题而提出来的,在Grandyang的百度空间有详细的公式和算法说明http://hi.baidu.com/grandyang/item/e024400bf3c879ce915718e2)变换,得到Q1R1
新矩阵A1=R1Q1,如果A1的主迹等于|A|停止迭代,否则继续
-----(镜像变换)---->Q2R2
A2=R2Q2,如果A2的主迹等于|A|停止迭代,否则继续
-----(镜像变换)---->Q3R3
A3=R3Q3 ....
...
Ak=RkQk
...
An=RnQn
QR法经过N次迭代总能得到收敛,收敛速度依赖于相邻特征值的比值,为了加速收敛可以引入原点平移法,此时迭代公式改为Ak-sk*I=RkQk+sk*I , sk为位移量。求得的所有特征值后,采用幂法的第4步求特征向量。
QR法可以把求解精度设得很高,如1e-10,所需增加的计算量并不是太多(一个5x5方阵,计算精度由1e-8提高到1e-10其迭代次数也仅是由24次提高到28次),当然由于VLISP自身的计算精度并不是很高(如乘*、除/、幂expt等函数),精度也不适合设得太高,建议采用1e-10。值得一提的是QR法同时还可以用来求超秩矩阵(超定线性方程组的系数)的特征值。

  1. ;;;----------------------------------Householder变换------------------------------;;;
  2. (defun householder  (mat / imat i m n d a b u hk q)
  3.   ;; Use Householder transformation to Get QR matrix for QR-decomposition
  4.   ;; (householder mat)
  5.   ;; (setq mat '((2. 3.) (1. 2.) (5. 7.)))
  6.   ;; (setq mat '((1. 2. -3.) (1. 3. 10.) (1. -1. 6.)))
  7.   (setq i    0
  8. n    (length (car mat))
  9. m    (length mat)
  10. imat ([i] m)
  11. q    imat)
  12.   (while (< i n)
  13.     (setq d (mapcar (function (lambda (x) (nth i x))) mat))
  14.     (repeat i
  15.       (setq d (cdr d)))
  16.     (setq a (* (sign (car d))
  17.         (sqrt (apply (function +)
  18.        (mapcar (function (lambda (x)
  19.       (* x x)))
  20.         d))))
  21.    b (* a (+ a (car d)))
  22.    d (cons (+ (car d) a)
  23.     (cdr d)))
  24.     (repeat i
  25.       (setq d (cons 0. d)))
  26.     (setq u (vtxv d d))
  27.     (if (not(zerop b))      
  28.       (setq u ([k*] (1/ b) u)))
  29.     (setq hk ([-] imat u)
  30.    q  ([*] q hk)
  31.    mat ([*] hk mat)
  32.    i (1+ i)))
  33.   (list q mat))
  34. ;;;------------------------------------QR分解迭代---------------------------------;;;
  35. ;;
  36. (defun QR->EV  (mat eps / d i b t1)
  37.   ;; QR-decomposition
  38.   ;; by GSLS(SS) 2013-1-2
  39.   (setq t1 (getvar "millisecs")
  40. d  (gc:determ mat)
  41. i  0);_Set the maximum number of iterations to prevent stack overflow
  42.   (while (and (not (equal d ([tr] mat) eps)) (< i 1000))
  43.     (setq b   (householder mat)
  44.    mat ([*] (cadr b) (car b))
  45.    i   (1+ i)))
  46.   (princ (strcat "\nQR-decomposition of iterations is : "
  47.    (rtos i 2)
  48.    " times ."))
  49.   (princ (strcat "\nTotel time of QR-decomposition : "
  50.    (rtos (- (getvar "millisecs") t1) 2 0)
  51.    " m.s."))
  52.   mat)

Givens变换(这是20世纪90年代的新算法,又称平面旋转变换) 是一种重要的正交交换, 它具有计算量小, 数值稳定的优点;欲将一个向量许多分量化为0采用Householder变换;欲将一个向量某个分量化为0采用Givens变换。Givens变换向量解和迭代法求解特征值和特征向量,主要步骤是:
1.先将矩阵A做右旋转相似变换,依次将从左到右的各列下三角元素变换为0,在这个过程中,下三角元素会逐步趋于0,因此需要多次迭代,将矩阵转换为上三角矩阵S则各个对角线上的元素即为所有特征值。
2.之后可以用幂法的第4步来求解特征向量。
  1. ;;--------------------------------Givens变换法-----------------------------------;;;
  2. (defun Givens  (mat Q / i j aij aii ajj aji ang a b m n0 n r Q1)
  3.   ;;(Givens a)
  4.   ;; Givens 变换求上三角阵
  5.   (setq i  (length mat)
  6. m0 i
  7. n0 (length (car mat))
  8. j -1)
  9.   (while (< j (1- n0))
  10.     (setq i m0
  11.    j (1+ j))
  12.     (while (> (setq i (1- i)) j)
  13.       (setq aij (nth j (nth i mat)))
  14.       (if (not (equal aij 0.0 1e-14))
  15. (progn
  16.    (setq aii (nth i (nth i mat))
  17.   ajj (nth j (nth j mat))
  18.   aji (nth i (nth j mat)))
  19.    (if (= aii ajj)
  20.      (if (> aij 0)
  21.        (setq ang _pi4)
  22.        (if (< aij 0)
  23.   (setq ang (- _pi4))))
  24.      (setq ang  (/ (atan (/  aij (- ajj aii) 0.5)) 2.))
  25.      )
  26.    (setq ang (max (min ang _pi4) (- _pi4)))
  27.    (setq m m0
  28.   Q1 nil)
  29.    (while (<= 0 (setq m (1- m)))
  30.      (setq n n0
  31.     r nil)
  32.      (while (<= 0 (setq n (1- n)))
  33.        (setq r (cons (cond ((and (= m i) (= n j))
  34.        (- (sin ang)))
  35.       ((and (= m j) (= n i))
  36.        (sin ang))
  37.       ((and (= m i) (= n i))
  38.        (cos ang))
  39.       ((and (= m j) (= n j))
  40.        (cos ang))
  41.       ((= m n)
  42.        1.)
  43.       (0.0))
  44.        r)))
  45.      (setq Q1 (cons r Q1))
  46.      )
  47.    (setq Q ([*] Q1 Q))
  48.    (setq mat ([*] ([*] Q1 mat) ([trp] Q1)))
  49.    )
  50. )      
  51.       )
  52.     )
  53.   (list mat Q))
  54. ;;-------------------
  55. (defun Givens->EV  (mat eps / d i b t1 q )
  56.   ;; Givens method main function
  57.   ;; by GSLS(SS) 2013-1-4
  58.   ;; (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)))
  59.   ;; (Givens->EV a 1e-12) --> (0.0660769 4.18334 18.2577 94.9829 7.50996)
  60.   (setq t1 (getvar "millisecs")
  61. d  (gc:determ mat)
  62. q  ([i] (length (car mat)))
  63. i  0 ;_Set the maximum number of iterations to prevent stack overflow
  64. )
  65.   (while (and (not (equal d ([tr] mat) eps)) (< i 1000))
  66.     (setq b (Givens mat Q)
  67.    mat (car b)
  68.    q (cadr b))
  69.     (setq i (1+ i)))
  70.   (print mat)
  71.   (princ (strcat "\nJacobi iterations is : "
  72.    (rtos i 2)
  73.    " times ."))
  74.   (princ (strcat "\nTotel time of Jacobi solving : "
  75.    (rtos (- (getvar "millisecs") t1) 2 0)
  76.    " m.s."))
  77.   (setq i -1)  
  78.   (mapcar (function (lambda (r1 )        
  79.         (nth (setq i (1+ i)) r1)))
  80.    mat)
  81.   )

仅求全部特征值而言,同精度下比较,Givens变换的迭代次数仅为QR分解迭代的1/5左右,总计算量约为QR法的一半,效率高了很多。
当然,Givens变换Jacobi方法结合起来,就可以同时求出特征值和特征向量;但是这里编程起来比较麻烦,需要逐次迭代将从大到小得非对角线上元素采用Givens变换消0;并同时记录变换矩阵的乘积Q=Q*Qi,上三角矩阵S主对角线上的特征值对应的变换矩阵积Q的各列为特征向量;这种方法也可以算是同步求解法之一,据我所知,该法是目前应用最为广泛的方法。


同步求解法首先得列并阵,行式[ (λiI-A)^t | I ]和列式[ λiI-A / I ]都可以;对于行式采用初等行变换将左阵转换为上三角矩阵,主对角线上的元素即为特征值,特征所对应右阵的行即为特征向量;对于列式采用初等列变换,将上阵转换为下三角方阵,主对角线上的元素即为特征值,其对应下阵的各列即为特征向量。该法最为精准,但是对于大型矩阵的求解,要书写计算式就麻烦了。同步求解法,还有其他的方法来实现,如Jacobi叠代法、行列互换法

下面是特征值与特征向量求解的主要函数,请下载附件文件做测试。
  1. ;;;-------------------------------------------------------------------------------;;;
  2. ;;;                         Get EigenValues and EigenVectors                      ;;;
  3. ;;;-------------------------------------------------------------------------------;;;
  4. (defun [E]  (mat / em i ai aii el vl res)
  5.   ;; get Eigenvalues and Eigenvectors
  6.   ;;
  7.   ;; First Use QR-decomposition to get all EigenValues ,
  8.   ;; Then Use Gaussian-elimination to get all EigenVectors .
  9.   ;;
  10.   ;; mat ---- a nonsingular matrix
  11.   ;;
  12.   ;; returns a list of all Eigenvalues and Eigenvectors , like '((Value1 Vector1) (Value2 Vector2) ... (Valuen Vectorn))
  13.   ;;
  14.   ;; e.g.
  15.   ;; (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))
  16.   ;;
  17.   ;; (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))
  18.   ;;
  19.   ;; (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))
  20.   ;;  ******if use Power method , this one will too slow .******
  21.   ;;
  22.   ;; (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))
  23.   ;;
  24.   ;; (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)))
  25.   ;;          ;_([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)
  26.   ;;        (4.18334 4.56872 -4.5105 3.56226 -1.68936 1.0) (0.0660769 6.41687 3.11375 -4.83309 -0.559052 1.0))
  27.   ;;
  28.   ;; (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))
  29.   ;;  ******if use Power method , this one will too slow .******
  30.   ;;
  31.   ;; (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))
  32.   ;;  ******if use Power method , this one will too slow .******
  33.   ;;        
  34.   ;;        get:     (setq vl ([E] a) v (nth 1 vl))
  35.   ;;                 
  36.   ;;        checkEg: (equal (gc:determ ([-] ([λ] (length (car a)) (car v)) a)) 0.0 1e-8)
  37.   ;;                 
  38.   ;;        checkEv  (equal (mxv a (cdr v)) (mxv ([λ] (length a) (car v)) (cdr v)) 1e-8)
  39.   ;;
  40.   ;; by GSLS(SS) 2013-1-3
  41.   ;;
  42.   ;; Any suggestions are welcome , and please post them to chlh_jd@126.com or 56454300@qq.com
  43.   ;;
  44.   ;; Thanks Robert for method suggestions .
  45.   ;; Thanks HighFlyingBird's help .
  46.   ;; Thanks ElpanovEvgeniy's Gaussian-elimination function , it's really a pretty codes .  ;;----------- QR method  ------------
  47.   (setq Em (QR->EV mat 1e-12));_1e--10 The tolerance between the matrix modules with the product of it's Trace ,
  48.                               ;_You can set a larger value like 1e-8 to reduce computations .
  49.     (setq i 0)
  50.     (while (and (setq ai (nth i Em)) (setq aii (nth i ai)))
  51.     (setq el (cons aii el)
  52.    i (1+ i)))                          
  53.   ;;-----------------------------------
  54.   
  55.   ;;----------- Power method ----------
  56.   ;_(setq el (Power->EV mat 1e-9))
  57.   ;;-----------------------------------
  58.   
  59.   (foreach a  el  ;_Get all EigenVectors
  60.     (setq vl (EigenValue->EigenVectors a mat))
  61.     (if vl
  62.       (foreach c  vl
  63. (setq res (cons (cons a c) res)))
  64.       (setq res (cons (cons a nil) res))))
  65.   
  66.   ;;----------- QR method--------------
  67.   (mapcar (function(lambda (x)
  68.        (mapcar (function (lambda (y)
  69.       (if (equal y 0.0 1e-9)
  70.         0.0
  71.         y))) x)))
  72.         res)
  73.   ;;-----------------------------------  ;;----------- Power method ----------
  74.   ;|(reverse (mapcar (function(lambda (x)
  75.        (mapcar (function (lambda (y)
  76.       (if (equal y 0.0 1e-9)
  77.         0.0
  78.         y))) x)))
  79.         res))|;
  80.   ;;-----------------------------------
  81.   
  82.   )
  83. ;;
  84. (defun EigenValue->EigenVectors  (λ mat / n imat ires)
  85. ;_(setq mat a λ -2. ires nil)
  86. ;_(if (equal λ 4. 1e-4)(princ))
  87. ;_(princ (rtos λ 2 8))
  88.   (setq n    (length mat)
  89. imat ([i] n)
  90. mat  ([trp] ([-] ([λ] n λ) mat))
  91. mat  (mapcar (function (lambda (x y)
  92.      (append x y)))
  93.        mat
  94.        imat);_Creat matrix [(λI-A)^t | I]
  95. mat  (eea:gauss mat);_Gaussian elimination
  96. imat nil)
  97.   (setq
  98.     mat (mapcar (function (lambda (x / y)
  99.        (setq x (reverse x))
  100.        (repeat n
  101.          (setq y (cons (car x) y)
  102.         x (cdr x)))
  103.        (setq imat (cons (reverse y) imat))
  104.        x))
  105.   mat)
  106.     )
  107.   (setq mat (reverse mat))
  108.   ;; Check ZeroVector of Left Matrix .
  109.   (while (<= (apply (function +)
  110.       (mapcar (function *) (car mat) (car mat)))
  111.       5e-7)
  112.     (setq mat  (cdr mat)
  113.    ires (cons (car imat) ires)
  114.    imat (cdr imat)))
  115.   ;; determine Linealy Independent
  116.   (while (and (cadr mat) (not (vslip (car mat) (cadr mat))))
  117.     (setq
  118.       ires (cons
  119.       (vrefzero (car imat) (cadr imat) (car mat) (cadr mat))
  120.       ires)
  121.       imat (cdr imat)
  122.       mat  (cdr mat)))
  123.   ;; Result deal ...
  124.   (mapcar
  125.     (function reverse)
  126.     (mapcar (function (lambda (b)
  127.    (if (equal b 0.0 5e-7)
  128.      0.0
  129.      b)))
  130.      ires))
  131.   )
  132. ;;;-------------------------------------高斯消元法--------------------------------;;;
  133. ;;;
  134. (defun eea:gauss  (lst)
  135.   ;; Gauss eliminationby
  136.   ;; by ElpanovEvgeniy
  137.   ;; http://www.theswamp.org/index.php?action=post;quote=163485;topic=13505.13505.msg163485#msg163485
  138.   ;; Implementation Gaussian elimination
  139.   ;; For text:
  140.   ;;  1x+2y+3z=2
  141.   ;;  10x+1y+8z=17
  142.   ;;  7z+2y=5
  143.   ;; (eea:gauss '((1.0 2.0 3.0) (10.0 1.0 8.0) (0.0 2.0 7.0)))
  144.   ;; =>
  145.   ;; ((1.0 2.0 3.0 2.0) (0.0 -19.0 -22.0 -3.0) (0.0 0.0 4.68421 4.68421))
  146.   ;;(gauss lst)
  147.   (if (car lst)
  148.     (if ;|(zerop (caar lst))|;  (equal (caar lst) 0.0 1e-8) ;_here changed for [E]
  149.       (if (vl-every ;|(function zerop)|;
  150.      (function (lambda (x) (equal x 0.0 1e-8))) ;_here changed for [E]
  151.      (mapcar (function car) lst))
  152. (if (cdr lst)
  153.    (cons
  154.      (car lst)
  155.      (mapcar
  156.        (function (lambda (x) (cons 0. x)))
  157.        (eea:gauss (mapcar (function cdr) (cdr lst)))
  158.        ) ;_  mapcar
  159.      ) ;_  cons
  160.    lst
  161.    ) ;_  if
  162. (eea:gauss
  163.    (cons
  164.      (mapcar
  165.        (function +)
  166.        (car lst)
  167.        (car (vl-remove-if
  168.        (function (lambda (x) (zerop (car x))))
  169.        (cdr lst)))
  170.        ) ;_  mapcar
  171.      (cdr lst)
  172.      ) ;_  cons
  173.    ) ;_  gauss
  174. ) ;_  if
  175.       (cons
  176. (car lst)
  177. (mapcar
  178.    (function (lambda (x) (cons 0. x)))
  179.    (eea:gauss
  180.      (mapcar
  181.        (function
  182.   (lambda (x / i)
  183.     (if (and (equal (car x) 0.0 1e-6)
  184.       (equal (caar lst) 0.0 1e-6)) ;_here add '0 / 0 = 1.'
  185.       (setq i 1.)
  186.       (setq i (/ (car x) (caar lst))))
  187.     ;|(setq i (/ (car x) (caar lst)))|;
  188.     (mapcar
  189.       (function -)
  190.       (cdr x)
  191.       (mapcar (function (lambda (x1) (* x1 i)))
  192.        (cdar lst))
  193.       );_  mapcar
  194.     );_  lambda
  195.   );_  function
  196.        (cdr lst)
  197.        );_  mapcar
  198.      );_  test
  199.    );_  mapcar
  200. );_  cons
  201.       );_  if
  202.     );_  if
  203.   );_  defun

非常感谢Robert和HighFlyingBird斑竹的热心帮助,同时也致谢ElpanovEvgeniy的高斯消元函数、Gile的矩阵高斯求模函数。
希望大家多多指导
  1. (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 '((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))

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

x

评分

参与人数 3明经币 +6 金钱 +60 收起 理由
tigcat + 1 很给力!
qjchen + 2 + 30 很给力!
highflybir + 3 + 30 赞一个!

查看全部评分

发表于 2022-6-10 17:29:35 | 显示全部楼层
landsat99 发表于 2022-6-8 18:46
eigen类的lsp实现,难得见到,赞。数值计算在commonlisp或scheme这类lisp中有什么优势特点吗?不了解lisp, ...

很久之前写过一些 矩阵求逆什么的 求方程解的,估计改改也能算特征值和特征向量
http://www.theswamp.org/index.php?topic=13505.msg163149#msg163149



但是可能是水平问题,就写得很辛苦,而autolisp的效率也比较低。所以个人觉得假如计算简单的有限元,几百根杆件以内的,是没有什么问题,非线性的可能就比较麻烦了。可能用objectarx会好些,3d3s就是这么做的。

不知道纯lisp会否好些,作为一个和fortran一样古老的语言活到现在,肯定有它厉害的地方。

帖子中的Evgeniy就是一个大高手,好像用common Lisp写过一个系统。


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

x
发表于 2022-6-13 12:00:55 | 显示全部楼层
qjchen 发表于 2022-6-10 17:29
很久之前写过一些 矩阵求逆什么的 求方程解的,估计改改也能算特征值和特征向量
http://www.theswamp.or ...

资深大神!! 深耕lsp十几年,大赞一个

这个矩阵运算的UI方式 太独特新颖了,创意很受启发。。。  相当精彩,果断收藏!
lsp果然神一般的存在于CAD
发表于 2022-6-8 18:46:03 来自手机 | 显示全部楼层
eigen类的lsp实现,难得见到,赞。数值计算在commonlisp或scheme这类lisp中有什么优势特点吗?不了解lisp,,请大神科普下
 楼主| 发表于 2013-1-4 19:56:54 | 显示全部楼层
本帖最后由 chlh_jd 于 2013-1-5 16:23 编辑

修改了下Jacobi函数,可以同时得到矩阵的特征值与特征向量,效率比QR法高了非常多;
使用的是依次对下三角元素消0叠代,每次迭代都必须进行下三角检查0检查,如果与0不在误差范围内则必须进行矩阵计算;
一般迭代收敛到1e-10需要n次迭代,n为方阵行数;每次迭代的最大矩阵乘法运算为1.5(n-1)^2次,但总的矩阵乘法次数有望在
[1.5+1.5*0.5+1.5*0.5^2+...+1.5*0.5^(n-1)](n-1)^2--->n=∞时为3(n-1)^2次以内;算法的复杂程度低于QR法(后续);且可以改进Givens一次变换的元素数量,设置一个阀值,仅对于元素大于阀值的所有下三角元素进行矩阵运算,依次减小阀值,直到阀值足够小。

  1. (defun Givens  (mat Q / i j aij aii ajj aji ang a b m n0 n r Q1)
  2.   ;;(Givens a)
  3.   ;; Givens 变换求上三角阵
  4.   ;; by GSLS(SS) 2013-1-4
  5.   (setq i  (length mat)
  6. m0 i
  7. n0 (length (car mat))
  8. j  -1)
  9.   (while (< j (1- n0))
  10.     (setq i m0
  11.    j (1+ j))
  12.     (while (> (setq i (1- i)) j)
  13.       (setq aij (nth j (nth i mat)))
  14.       (if (not (equal aij 0.0 1e-14))
  15. (progn
  16.    (setq aii (nth i (nth i mat))
  17.   ajj (nth j (nth j mat))
  18.   aji (nth i (nth j mat)))
  19.    (if (= aii ajj)
  20.      (if (> aij 0)
  21.        (setq ang _pi4)
  22.        (if (< aij 0)
  23.   (setq ang (- _pi4))))
  24.      (setq ang (/ (atan (/ aij (- ajj aii) 0.5)) 2.))
  25.      )
  26.    (setq ang (max (min ang _pi4) (- _pi4)))
  27.    (setq m  m0
  28.   Q1 nil)
  29.    (while (<= 0 (setq m (1- m)))
  30.      (setq n n0
  31.     r nil)
  32.      (while (<= 0 (setq n (1- n)))
  33.        (setq r (cons (cond ((and (= m i) (= n j))
  34.        (- (sin ang)))
  35.       ((and (= m j) (= n i))
  36.        (sin ang))
  37.       ((and (= m i) (= n i))
  38.        (cos ang))
  39.       ((and (= m j) (= n j))
  40.        (cos ang))
  41.       ((= m n)
  42.        1.)
  43.       (0.0))
  44.        r)))
  45.      (setq Q1 (cons r Q1))
  46.      )
  47.    (setq Q ([*] Q1 Q))
  48.    (setq mat ([*] ([*] Q1 mat) ([trp] Q1)))
  49.    )
  50. )
  51.       )
  52.     )
  53.   (list mat Q))
  54. ;;-------------------------------------------------------------------------------;;;
  55. (defun Jacobi  (mat eps / d i b t1 q)
  56.   ;; Jacobi method
  57.   ;; by GSLS(SS) 2013-1-4
  58.   ;;
  59.   ;; (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))
  60.   ;;
  61.   ;; (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))) -->
  62.   ;;      -->((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)
  63.   ;;           (94.9829 0.100149 0.192449 0.29428 0.723881 0.585089) (7.50996 -0.0487036 0.595184 0.450894 -0.586047 0.310849))
  64.   ;;
  65.   ;; (Jacobi a 1e-12)
  66.   (setq t1 (getvar "millisecs")
  67. d  (gc:determ mat)
  68. q  ([i] (length (car mat)))
  69. i  0 ;_Set the maximum number of iterations to prevent stack overflow
  70. )
  71.   (while (and (not (equal d ([tr] mat) eps)) (< i 1000))
  72.     (setq b   (Givens mat Q)
  73.    mat (car b)
  74.    q   (cadr b))
  75.     (setq i (1+ i)))
  76.   (princ (strcat "\nJacobi iterations is : "
  77.    (rtos i 2)
  78.    " times ."))
  79.   (princ (strcat "\nTotel time of Jacobi solving : "
  80.    (rtos (- (getvar "millisecs") t1) 2 0)
  81.    " m.s."))
  82.   (setq i -1)
  83.   (mapcar (function (lambda (a b)
  84.         (cons a b)))
  85.    (mapcar (function (lambda (r1)
  86.          (nth (setq i (1+ i)) r1)))
  87.     mat)
  88.    q)
  89.   )
 楼主| 发表于 2013-1-5 16:22:06 | 显示全部楼层
本帖最后由 chlh_jd 于 2013-1-5 16:49 编辑

改进QR法函数代码,可以同时求出特征值与特征向量;与Jacobi法相比,要达到准确的特征向量解,需要把求模精度提高1倍,否则求到的特征向量会反号,尽管如此得到的特征向量的精度任然小于Jacobi法,这是因为Jacobi法采用的Givens平面旋转变换比Householder正交变换的数值计算来得稳定。Householder每运算一遍要作3n次矩阵乘法,即QR法每迭代一次需要做3n次矩阵乘法,对于N维方阵,QR法要迭代5~6N次才能得到稳定向量解,则QR法需要共需要15~18n^2次矩阵乘法;与Jacobi法相比,
15~18n^2/[1.5n(n-1)^2] =10~12n倍,如果Jacobi法临近下限3(n-2)^2次矩阵迭代,则比值变为5~6n^2。可见QR法的效率明显低于Jacobi法。
当然QR法也可以仅把方阵变成海森伯格(Hessenberg)上三角矩阵,这样来减少计算量:Ak=QkRk,Ak+1=RkQk。
  1. ;;;----------------------------------Householder变换------------------------------;;;
  2. (defun householder  (q mat / i m n d a b u hk)
  3.   ;; Use Householder transformation to Get QR matrix for QR-decomposition
  4.   ;; (householder mat)
  5.   ;; (setq mat '((2. 3.) (1. 2.) (5. 7.)))
  6.   ;; (setq mat '((1. 2. -3.) (1. 3. 10.) (1. -1. 6.)))
  7.   (setq i 0
  8. n (length (car mat))
  9. )
  10.   (while (< i n)
  11.     (setq d (mapcar (function (lambda (x) (nth i x))) mat))
  12.     (repeat i
  13.       (setq d (cdr d)))
  14.     (setq a (* (sign (car d))
  15.         (sqrt (apply (function +)
  16.        (mapcar (function (lambda (x)
  17.       (* x x)))
  18.         d))))
  19.    b (* a (+ a (car d)))
  20.    d (cons (+ (car d) a)
  21.     (cdr d)))
  22.     (repeat i
  23.       (setq d (cons 0. d)))
  24.     (setq u (vtxv d d))
  25.     (if (not (zerop b))
  26.       (setq u ([k*] (1/ b) u)))
  27.     (setq hk  ([-] imat u)
  28.    q   ([*] q hk)
  29.    mat ([*] ([*] hk mat) ([trp] hk))
  30.    i   (1+ i)))
  31.   (list q mat))
  32. ;;;------------------------------------QR分解迭代---------------------------------;;;
  33. ;;
  34. (defun QR->EV  (mat eps / imat q d i b t1)
  35.   ;; QR-decomposition
  36.   ;; by GSLS(SS) 2013-1-2
  37.   ;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.)))
  38.   ;; (QR->EV a 5e-11)
  39.   (setq t1   (getvar "millisecs")
  40. d    (gc:determ mat)
  41. imat ([i] (length (car mat)))
  42. q    imat
  43. i    0) ;_Set the maximum number of iterations to prevent stack overflow
  44.   (while (and (not (equal d ([tr] mat) eps)) (< i 1000))
  45.     (setq b   (householder q mat)
  46.    mat (cadr b)
  47.    q   (car b)
  48.    i   (1+ i)))
  49.   (princ (strcat "\nQR-decomposition of iterations is : "
  50.    (rtos i 2)
  51.    " times ."))
  52.   (princ (strcat "\nTotel time of QR-decomposition : "
  53.    (rtos (- (getvar "millisecs") t1) 2 0)
  54.    " m.s."))
  55.   (setq i -1)
  56.   (mapcar (function (lambda (a b)
  57.         (cons a b)))
  58.    (mapcar (function (lambda (r1)
  59.          (nth (setq i (1+ i)) r1)))
  60.     mat)
  61.    ([trp] q))
  62.   )
发表于 2013-1-6 19:25:39 | 显示全部楼层
矩阵方面很好的资料,学习了。
发表于 2013-1-8 08:01:58 | 显示全部楼层
东西很好,就是太深了。有没有简单实用的东西呢?
 楼主| 发表于 2013-1-9 01:35:17 | 显示全部楼层
矩阵特征值与特征向量的应用在现代科学里面已经非常广泛了:图像问题,趋势问题,矩阵计算机...
朱小平老师有一本书叫《实对称矩阵的拟特征值理论与应用》(内容摘要
这里有介绍2个应用例子:天气Markov链的稳定状态、人口流动问题;
这里介绍了3个应用例子:动物繁殖的规律问题、商品的市场占有率问题、常染色体遗传问题

但是在ACAD上具体应用,我还没发现;我最近一直在想,二维边界映射问题是否可以通过特征向量矩阵来求解,如果可以那么任意边界1范围内的图形就可以变换到任意边界2内了。
发表于 2013-5-22 12:31:34 | 显示全部楼层
太深了,完全看不懂
发表于 2013-6-22 09:32:37 | 显示全部楼层
矩阵很强大,,,
发表于 2018-3-23 09:27:06 | 显示全部楼层
向老师学习,感谢分享
发表于 2019-7-26 13:37:56 | 显示全部楼层

感谢楼主分享,收藏了
您需要登录后才可以回帖 登录 | 注册

本版积分规则

小黑屋|手机版|CAD论坛|CAD教程|CAD下载|联系我们|关于明经|明经通道 ( 粤ICP备05003914号 )  
©2000-2023 明经通道 版权所有 本站代码,在未取得本站及作者授权的情况下,不得用于商业用途

GMT+8, 2024-11-6 07:45 , Processed in 0.254983 second(s), 30 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表