有没有更快的求逆矩阵的方法

dnbcgrass · 发表于 2013-6-10 18:09:13

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？立即注册

×

本帖最后由 dnbcgrass 于 2013-6-10 18:10 编辑

下面是生成n阶非奇异矩阵(n×n阶方阵)函数

;======================== Cszdfz ==================================================================
;; 用法:
;; (Cszdfz 矩阵的阶数XN)
;; 返回值:
;; (n+2)/(2n+2) -1/2 0 0 ... 0 1/(2n+2)
;; -1/2 1 -1/2 0 ... 0 0
;; 0 -1/2 1 -1/2 ... 0 0
;; ..............................................
;; 0 ... -1/2 1 -1/2
;; 1/(2n+2) 0 ..............................-1/2 (n+2)/(2n+2)
(defun Cszdfz(XN / B i j zb dx)
(setq B nil i 1)
(repeat xn
(setq j 1 zb nil)
(repeat xn
(if (= i j)
(progn
(setq dx 1.0)
(if (or (= i 1) (= i xn))
(setq dx (/ (+ xn 2.0) (+ (* 2.0 xn) 2.0)))
)
)
(if (or (= (1- i) j) (= (1+ i) j))
(setq dx -0.5)
(progn
(setq dx 0.0)
(if (or (and (= i 1) (= j xn))
(and (= i xn) (= j 1))
)
(setq dx (/ 1.0 (+ (* 2.0 xn) 2.0)))
)
)
)
)
(setq zb (cons dx zb)
j (1+ j)
)
)
(setq B (cons (reverse zb) B)
i (1+ i)
)
)
(setq B (reverse B))
)

;========================   Cszdfz   ================================================================== 
;; 用法: 
;;        (Cszdfz 矩阵的阶数XN) 
;; 返回值: 
;;    (n+2)/(2n+2)  -1/2   0      0     ...     0     1/(2n+2) 
;;    -1/2          1      -1/2   0     ...     0     0 
;;    0             -1/2   1      -1/2  ...     0     0 
;;    .............................................. 
;;    0  ...                              -1/2  1   -1/2 
;;    1/(2n+2) 0 ..............................-1/2  (n+2)/(2n+2) 
(defun Cszdfz(XN  / B i j zb dx) 
     (setq B nil i 1) 
     (repeat xn 
         (setq j 1 zb nil) 
         (repeat xn 
             (if (= i j) 
                 (progn 
                     (setq dx 1.0) 
                     (if (or (= i 1) (= i xn)) 
                         (setq dx (/ (+ xn 2.0) (+ (* 2.0 xn) 2.0))) 
                     ) 
                 ) 
                 (if (or (= (1- i) j) (= (1+ i) j)) 
                     (setq dx -0.5) 
                     (progn 
                         (setq dx 0.0) 
                         (if (or (and (= i 1) (= j xn)) 
                                 (and (= i xn) (= j 1)) 
                             ) 
                             (setq dx (/ 1.0 (+ (* 2.0 xn) 2.0))) 
                         ) 
                     ) 
                 ) 
             ) 
             (setq zb (cons dx zb) 
                   j  (1+ j) 
             ) 
        ) 
        (setq B (cons (reverse zb) B) 
              i (1+ i) 
        ) 
     ) 
     (setq B (reverse B)) 
)

下面是求n阶非奇异矩阵的逆矩阵的函数，好像是由Highflybird提供，出自何处记不得清了，
用上面生成的矩阵，采用该函数(Inverse)求得的逆矩阵为：
n n-1  n-2 ... 1 2
n-1  n n-1 ... 3 2
n-2  n-1  n ... 4 3
..........................
2 3 4 ... n n-1
1 2 3 ... n-1  n
但是，当n大于500以上时，速度很慢，大家有没有更快的方法？

;============= Inverse ========================================================================
;; Matrix Inverse - gile & Lee Mac
;; Uses Gauss-Jordan Elimination to return the inverse of a non-singular nxn matrix.
;; Args: m - nxn matrix
(defun Inverse ( m / c f p r d)
(defun ff ( p m )
(mapcar '(lambda ( x ) (mapcar '(lambda ( a b ) (- a (* (car x) b))) (cdr x) p)) m)
)
(setq m (mapcar 'append m (imat (length m))));原矩阵单位矩阵
(while m
(setq c (mapcar '(lambda ( x ) (abs (car x))) m))
(repeat (vl-position (apply 'max c) c)
(setq m (append (cdr m) (list (car m))))
)
(if (equal 0.0 (caar m) 1e-14)
(setq m nil
r nil
)
(setq p (mapcar '(lambda ( x ) (/ (float x) (caar m))) (cdar m))
m (ff p (cdr m))
r (cons p (ff p r))
)
)
)
(reverse r)
)
;;============= imat ===========================================================================
;; Identity Matrix - Lee Mac
;; Args: n - matrix dimension
(defun imat ( n / i j l m )
(repeat (setq i n)
(repeat (setq j n)
(setq l (cons (if (= i j) 1.0 0.0) l)
j (1- j)
)
)
(setq m (cons l m)
l nil
i (1- i)
)
)
m
)

;=============   Inverse  ======================================================================== 
;; Matrix Inverse  -  gile & Lee Mac 
;; Uses Gauss-Jordan Elimination to return the inverse of a non-singular nxn matrix. 
;; Args: m - nxn matrix 
(defun Inverse ( m / c f p r d) 
    (defun ff ( p m ) 
        (mapcar '(lambda ( x ) (mapcar '(lambda ( a b ) (- a (* (car x) b))) (cdr x) p)) m) 
    ) 
   
    (setq m (mapcar 'append m (imat (length m))));原矩阵 单位矩阵 
    (while m 
        (setq c (mapcar '(lambda ( x ) (abs (car x))) m)) 
        (repeat (vl-position (apply 'max c) c) 
            (setq m (append (cdr m) (list (car m)))) 
                  
        ) 
        (if (equal 0.0 (caar m) 1e-14) 
            (setq m nil 
                  r nil 
            ) 
            (setq p (mapcar '(lambda ( x ) (/ (float x) (caar m))) (cdar m)) 
                  m (ff p (cdr m)) 
                  r (cons p (ff p r)) 
            ) 
        ) 
    ) 
    (reverse r) 
) 
 
;;=============   imat  =========================================================================== 
;; Identity Matrix  -  Lee Mac 
;; Args: n - matrix dimension 
(defun imat ( n / i j l m ) 
    (repeat (setq i n) 
        (repeat (setq j n) 
            (setq l (cons (if (= i j) 1.0 0.0) l) 
                  j (1- j) 
            ) 
        ) 
        (setq m (cons l m) 
              l nil 
              i (1- i) 
        ) 
    ) 
    m 
)

851eb6ae7e2c4dd · 发表于 2013-6-10 19:54:51

在网上搜一搜，选个好的算法，只不过把其他语言转换成lisp好了。

XDSoft · 发表于 2013-6-10 20:03:16

LEE-MAC的Identity Matrix ，楼主测试了吗，效率如何？

XDSoft · 发表于 2013-6-10 20:19:23

XDSoft 发表于 2013-6-10 20:03
LEE-MAC的Identity Matrix ，楼主测试了吗，效率如何？

提供篇资料，只要是数学方法，应该不会太慢的，看看组合下LISP的函数，尽量使用基本函数。

下载地址：晓东文库

eachy · 发表于 2013-6-10 20:24:15

在AutoCAD中讨论广义逆矩阵没多大意义，作为知识了解还是必要的

Lispboy · 发表于 2013-6-10 20:30:16

楼主，这种N级的测试效率的，里面一定不能有APPEND，你的代码里面有APPEND,你看看用CONS和LISP组合能替代不，效率慢应该是APPEND的原因。

Free-Lancer · 发表于 2013-6-10 22:44:02

在AutoCAD下只要讨论4x4就够了来自: Android客户端

账号		自动登录	找回密码
密码			立即注册

[研讨] 有没有更快的求逆矩阵的方法

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分

点评