方阵的逆矩阵函数及求解线性方程

snoopychen

方阵求逆是数值计算中非常重要的一个步骤，包括方程求解，特征值计算等都和
逆矩阵密切相关，数值计算方法有许多，比如可以采用LU分解等原理。
eachy版主编过一个，好像是用于ucs坐标系转换的3*3矩阵的，对于n*n阶的Lisp
代码一时在网上没有找到。FORTRAN和C倒是都有源码库。
本来用lisp来进行这种计算似乎不是很好，网上一些大牛也建议不要用lisp进行
矩阵求逆什么的。
但鉴于对lisp非常喜欢，把j.p. griffith的fortran代码改写成如下的lisp格式
各位帮忙测试有没有错误。
程序不是很完善，应该在那些提示singular matrix的不满秩矩阵中加一段判断
之后就直接退出，不再repeat n
在编写这个的过程中，遇到了不少表和数组有点不同的地方
1）表格式比数组格式当然要自由和方便，但是，对表中某个数值进行修改的时候
似乎特别麻烦，就算本程序中采用了std-setnth的函数，感觉它在修改某个数值的
时候必须先将此位置之前的所有拿到一个表A中，然后添加这个修改的数值，再把
后面那部分的数值添加到表A中，这个思路感觉挺耗时间的。虽然这个std函数号称用了二分法进行查找，但还是觉得不大好。而在fortran中，我们只要设A(i,j)=number，就可以方便的代替某个位置的数值了。
2）对表中某个数字进行数值运算的时候，比如把它乘一个变量的时候，必须先将其
调出，然后修改，然后再用std-setnth修改表中的项。

基于以上原因，请大侠们看看如何比较好的定义这样两个函数

A）a是一个表，func是一个可以用Lambda之类定义的函数，i,j分别是二维表的行列
这个函数用来将二维表a的i行j列元素进行func（比如sin什么的）函数修改
(xd-abcde (func a i j))
B）存不存在直接修改表中第i项，不要其他辅助表的方法

[pcode=lisp,true]
;;; N*N方阵求逆算法
;;;by qjchen at www.xdcad.net
;;;------------------------------------------------------------------------
;;;根据j.p. griffith  6/88的fortran程序改编?
;;;利用了类似gauss-jordan的消元法原理进行编写
;;;缺点：对于不能求逆的矩阵没有很好容错

(defun c:test ()
  (setq t1 (list 3.0 3.0 3.0 3.0 3.0))
  (setq t2 (list 3.0 4.0 3.0 3.0 3.0))
  (setq t3 (list 3.0 3.0 5.0 3.0 3.0))
  (setq t4 (list 3.0 3.0 3.0 6.0 3.0))
  (setq t5 (list 3.0 3.0 3.0 3.0 8.0))
  (setq tt (list t1 t2 t3 t4 t5))
  (setq tt1 (matrix_inversion tt))
(princ "\n              tt = ") (princ tt);Erase_DV
(getstring "\nPress Enter to continue...");Erase_DV
  
(princ "\n             tt1 = ") (princ tt1);Erase_DV
(getstring "\nPress Enter to continue...");Erase_DV
)

(defun matrix_inversion (a / nmax ipiv indxr indxc n n i j k l ll a-temp
                           temp big irow icol pivinv
                        )
  (setq nmax 50)
  (setq ipiv (create nmax))
  (setq indxr (create nmax))
  (setq indxc (create nmax))
  (setq n (length a))
  (setq i 1)
  (repeat n
    (setq big 0)
    (setq j 1)
    (repeat n
      (if (/= ([n] ipiv j) 1)
        (progn
          (setq k 1)
          (repeat n
            (if (= ([n] ipiv k) 0)
              (progn
                (if (>= (abs ([n,m] a j k)) big)
                  (progn
                    (setq big ([n,m] a j k))
                    (setq irow j)
                    (setq icol k)
                  )
                  (progn
                    (if (> ([n] ipiv k) 1)
                      (princ "singular matrix")
                    )
                  )
                )
              )
            )
            (setq k (1+ k))
          )
        )
      )
      (setq j (1+ j))
    )
    (setq temp (1+ ([n] ipiv icol)))
    (setq ipiv (std-setnth temp (1- icol) ipiv))
    (if (/= irow icol)
      (progn
        (setq l 1)
        (repeat n
          (setq dum ([n,m] a irow l))
          (setq a-temp ([n,m] a icol l))
          (setq a (qj-setnmth a-temp (1- irow) (1- l) a))
          (setq a (qj-setnmth dum (1- icol) (1- l) a))
          (setq l (1+ l))
        )
      )
    )
    (setq indxr (std-setnth irow (1- i) indxr))
    (setq indxc (std-setnth icol (1- i) indxc))
    (if (= ([n,m] a icol icol) 0)
      (princ "singular matrix")
    )
    (setq pivinv (/ 1.0 ([n,m] a icol icol)))
    (setq a (qj-setnmth 1 (1- icol) (1- icol) a))
    (setq l 1)
    (repeat n
      (setq a-temp (* ([n,m] a icol l) pivinv))
      (setq a (qj-setnmth a-temp (1- icol) (1- l) a))
      (setq l (1+ l))
    )
    (setq ll 1)
    (repeat n
      (if (/= ll icol)
        (progn
          (setq dum ([n,m] a ll icol))
          (setq a (qj-setnmth 0.0 (1- ll) (1- icol) a))
          (setq l 1)
          (repeat n
            (setq a-temp (- ([n,m] a ll l) (* dum ([n,m] a icol l))))
            (setq a (qj-setnmth a-temp (1- ll) (1- l) a))
            (setq l (1+ l))
          )
        )
      )
      (setq ll (1+ ll))
    )
    (setq i (1+ i))
  )
  (setq l n)
  (repeat n
    (if (/= ([n] indxr l) ([n] indxc l))
      (progn
        (setq k 1)
        (repeat n
          (setq dum ([n,m] a k ([n] indxr l)))
          (setq a-temp ([n,m] a k ([n] indxc l)))
          (setq a (qj-setnmth a-temp (1- k) (1- ([n] indxr l)) a))
          (setq a (qj-setnmth dum (1- k) (1- ([n] indxc l)) a))
          (setq k (1+ k))
        )
      )
    )
    (setq l (1- l))
  )
  a
)


;;创建一个有m个0.0元素的表
(defun create (m / lista)               ; (setq m 5)
  (setq lista nil)
  (repeat m
    (setq lista (append
                  lista
                  (list 0.0)
                )
    )
  )
  lista
)

;创建一个m*m个0.0元素的表
(defun createm (m / lista listb)
  (setq m 5)
  (setq lista nil)
  (repeat m
    (setq listb nil)
    (repeat m
      (setq listb (append
                    listb
                    (list 0.0)
                  )
      )
    )
    (setq lista (append
                  lista
                  (list listb)
                )
    )
  )
  lista
)

;;; std代替表中第n个元素的函数
(defun std-%setnth (new i lst / fst len)
  (cond
    ((minusp i)
      lst
    )
    ((> i (setq len (length lst)))
      lst
    )
    ((> i (/ len 2))
      (reverse (std-%setnth new (1- (- len i)) (reverse lst)))
    )
    (t
      (append
        (progn
          (setq fst nil)               ; ; possible vl lsa compiler bug
          (repeat (rem i 4)
            (setq fst (cons (car lst) fst)
                  lst (cdr lst)
            )
          )
          (repeat (/ i 4)
            (setq fst (cons (cadddr lst) (cons (caddr lst) (cons
                                                                 (cadr lst)
                                                                 (cons
                                                                       (car lst)
                                                                       fst
                                                                 )
                                                           )
                                         )
                      )
                  lst (cddddr lst)
            )
          )
          (reverse fst)
        )
        (if (listp new)
          new
          (list new)
        )                               ; v0.4001
        (cdr lst)
      )
    )
  )
)
(defun std-setnth (new i lst)
  (std-%setnth (list new) i lst)
)


;代替二维表中第i行第j列元素的函数（i和j从1开始）
(defun qj-setnmth (new i j lst / listb lista)
  (setq listb lst)
  (setq lista (nth i lst))
  (setq lista (std-setnth new j lista))
  (setq listb (std-setnth lista i listb))
  listb
)

;;; 获取某个数组表第几项第几项的数值
(defun [n,m] (a n m / i)               ; n是行，m是列
  (setq i (nth (1- m) (nth (1- n) a)))
  i
)

;;; 获取某个单列数组第几项的数值
(defun [n] (a n / i)                       ; n是行，m是列
  (setq i (nth (1- n) a))
  i
)

[/pcode]

狂刀

先不管写的如何，你的精神是值得称道的

snoopychen

对，狂刀兄说的好，程序的编写确实比较罗嗦
可能多用点mapcar和apply会好很多，把lisp这个灵气十足的语言说得拖拖沓沓实在
是不应该。
不过这个函数由于是直接改编的，所以只是进行了直译（没有仔细读懂，不过还是用
visual fortran,lisp,matlab和excel都作过测试，矩阵运算是对的）
（感觉应该可以编一个程序，自动对fortran程序进行改编，这样子，许多的函数库就
可以自动转过来的了）

最近手痒，又改造了一下，加了一点std函数，改造成可以计算非齐次线性方程组的
小程序，见笑了。

程序如下，说明如内，演示如后
[php]
;;;线性方程求解by qjchen at www.xdcad.net
;;;说明：对于针对1x+2y+3z+0a-5b+6d=5之类从a-z的变量线性方程进行求解
;;;要求方程数量和未知量的个数相等，就是说
;;;;要求，所有方程的变量先后顺序一样，所有系数必须存在，包括0
;;;;最后答案的精度取3位，字高取250，可自己修改，答案中的x1,x2***对应相对位置的变量
;;;;缺点：对无解情况没有给出提示，直接出错

;;;main program
(defun c:test (/ lstall1 lstall2 n sstext obj ent text lsttext lstnum
                 lstnum1 lstnum2 i temp temp2
)
  (setq lstall1 nil
        lstall2 nil
  )
  (setq sstext (ssget))
  (if sstext
    (progn
      (setq n (sslength sstext)
            i 0
      )
      (repeat n
        (setq obj (ssname sstext i))
        (setq ent (entget obj))
        (setq text (cdr (assoc 1 ent)))
        (setq lsttext (STD-STRING->STRLIST text
                                           "xyzabcdefghijklmnopqrstuvw="
                      )
        )
        (setq lstnum (mapcar
                       '(lambda (x)
                          (atof x)
                        )
                       lsttext
                     )
        )
        (setq lstnum1 (reverse (cdr (reverse lstnum))))
        (setq lstnum2 (list (last lstnum)))
        (setq lstall1 (append
                        lstall1
                        (list lstnum1)
                      )
        )
        (setq lstall2 (append
                        lstall2
                        (list lstnum2)
                      )
        )
        (setq i (1+ i))
      )
    )
  )
  (setq temp (matrix_inversion lstall1))
  (setq temp2 (matrix-mul temp lstall2))
  (setq str "")       
  (setq tempstr (mapcar
                  '(lambda (x)
                     (rtos (nth 0 x) 2 3)
                   )
                  temp2
                )
  )
  (setq i 1)
  (foreach x tempstr
    (setq str (strcat str "x" (itoa i) "=" x " "))
    (setq i (1+ i))
  )
  (setq pp (getpoint "方程解的放置位置:\n"))
  (command "text" pp 2.5 0 str)
)

;;; N*N方阵求逆算法
;;; by qjchen at www.xdcad.net
;;; ------------------------------------------------------------------------
;;; 根据j.p. griffith  6/88的fortran程序改编?
;;; 利用了类似gauss-jordan的消元法原理进行编写
;;; 缺点：对于不能求逆
;;; N阶矩阵求逆的算法
(defun matrix_inversion (a / nmax ipiv indxr indxc n n i j k l ll a-temp
                           temp big irow icol pivinv
                        )
  (setq nmax 50)
  (setq ipiv (create nmax))
  (setq indxr (create nmax))
  (setq indxc (create nmax))
  (setq n (length a))
  (setq i 1)
  (repeat n
    (setq big 0)
    (setq j 1)
    (repeat n
      (if (/= ([n] ipiv j) 1)
        (progn
          (setq k 1)
          (repeat n
            (if (= ([n] ipiv k) 0)
              (progn
                (if (>= (abs ([n,m] a j k)) big)
                  (progn
                    (setq big ([n,m] a j k))
                    (setq irow j)
                    (setq icol k)
                  )
                  (progn
                    (if (> ([n] ipiv k) 1)
                      (princ "singular matrix")
                    )
                  )
                )
              )
            )
            (setq k (1+ k))
          )
        )
      )
      (setq j (1+ j))
    )

    (setq temp (1+ ([n] ipiv icol)))
    (setq ipiv (std-setnth temp (1- icol) ipiv))
    (if (/= irow icol)
      (progn
        (setq l 1)
        (repeat n
          (setq dum ([n,m] a irow l))

          (setq a-temp ([n,m] a icol l))
          (setq a (qj-setnmth a-temp (1- irow) (1- l) a))
          (setq a (qj-setnmth dum (1- icol) (1- l) a))
          (setq l (1+ l))
        )
      )
    )

    (setq indxr (std-setnth irow (1- i) indxr))
    (setq indxc (std-setnth icol (1- i) indxc))
    (if (= ([n,m] a icol icol) 0)
      (princ "singular matrix")
    )
    (setq pivinv (/ 1.0 ([n,m] a icol icol)))
    (setq a (qj-setnmth 1 (1- icol) (1- icol) a))
    (setq l 1)
    (repeat n
      (setq a-temp (* ([n,m] a icol l) pivinv))
      (setq a (qj-setnmth a-temp (1- icol) (1- l) a))
      (setq l (1+ l))
    )
    (setq ll 1)
    (repeat n
      (if (/= ll icol)
        (progn
          (setq dum ([n,m] a ll icol))
          (setq a (qj-setnmth 0.0 (1- ll) (1- icol) a))
          (setq l 1)
          (repeat n
            (setq a-temp (- ([n,m] a ll l) (* dum ([n,m] a icol l))))
            (setq a (qj-setnmth a-temp (1- ll) (1- l) a))
            (setq l (1+ l))
          )
        )
      )
      (setq ll (1+ ll))
    )
    (setq i (1+ i))
  )

  (setq l n)
  (repeat n
    (if (/= ([n] indxr l) ([n] indxc l))
      (progn
        (setq k 1)
        (repeat n
          (setq dum ([n,m] a k ([n] indxr l)))
          (setq a-temp ([n,m] a k ([n] indxc l)))
          (setq a (qj-setnmth a-temp (1- k) (1- ([n] indxr l)) a))
          (setq a (qj-setnmth dum (1- k) (1- ([n] indxc l)) a))
          (setq k (1+ k))
        )
      )
    )
    (setq l (1- l))
  )
  a
)


;;; 创建一个有m个0.0元素的表
(defun create (m / lista)               ; (setq m 5)
  (setq lista nil)
  (repeat m
    (setq lista (append
                  lista
                  (list 0.0)
                )
    )
  )
  lista
)

;;; 创建一个m*m个0.0元素的表
(defun createm (m / lista listb)
  (setq m 5)
  (setq lista nil)
  (repeat m
    (setq listb nil)
    (repeat m
      (setq listb (append
                    listb
                    (list 0.0)
                  )
      )
    )
    (setq lista (append
                  lista
                  (list listb)
                )
    )
  )
  lista
)


;;; stdlib函数---------------------------------------------------------------
(defun STD-STRING->LIST (s / lst)
  (if (stringp s)
    (while (/= s "")
      (setq lst (cons (ascii (substr s 1 1)) lst)
            s (substr s 2)
      )
    )
  )
  (reverse lst)
)
(defun STD-STRING->STRLIST (s delims)
  (std-strtok s delims)
)
(defun STD-STRTOK (s delims / len s1 i c lst)
  (setq delims (std-string->list delims)
        len (strlen s)
        s1 ""
        i (1+ len)
  )
  (while (> (setq i (1- i))
            0
         )
    (setq c (substr s i 1))
    (if (member (ascii c) delims)
      (if (/= s1 "")                       ; no null tokens
        (setq lst (cons s1 lst)
              s1 ""
        )
      )
      (setq s1 (strcat c s1))
    )
  )
  (if (/= s1 "")
    (cons s1 lst)                       ; no ("" "1" "2")!
    lst
  )
)
(defun STRINGP (s)
  (= (type s) 'STR)
)


;;; std代替表中第n个元素的函数
(defun std-%setnth (new i lst / fst len)
  (cond
    ((minusp i)
      lst
    )
    ((> i (setq len (length lst)))
      lst
    )
    ((> i (/ len 2))
      (reverse (std-%setnth new (1- (- len i)) (reverse lst)))
    )
    (t
      (append
        (progn
          (setq fst nil)               ; ; possible vl lsa compiler bug
          (repeat (rem i 4)
            (setq fst (cons (car lst) fst)
                  lst (cdr lst)
            )
          )
          (repeat (/ i 4)
            (setq fst (cons (cadddr lst) (cons (caddr lst) (cons
                                                                 (cadr lst)
                                                                 (cons
                                                                       (car lst)
                                                                       fst
                                                                 )
                                                           )
                                         )
                      )
                  lst (cddddr lst)
            )
          )
          (reverse fst)
        )
        (if (listp new)
          new
          (list new)
        )                               ; v0.4001
        (cdr lst)
      )
    )
  )
)
(defun std-setnth (new i lst)
  (std-%setnth (list new) i lst)
)


;;; 代替二维表中第i行第j列元素的函数（i和j从1开始）
(defun qj-setnmth (new i j lst / listb lista)
  (setq listb lst)
  (setq lista (nth i lst))
  (setq lista (std-setnth new j lista))
  (setq listb (std-setnth lista i listb))
  listb
)

;;; 获取某个数组表第几项第几项的数值
(defun [n,m] (a n m / i)               ; n是行，m是列
  (setq i (nth (1- m) (nth (1- n) a)))
  i
)

;;; 获取某个单列数组第几项的数值
(defun [n] (a n / i)                       ; n是行，m是列
  (setq i (nth (1- n) a))
  i
)


;;; ---------------------------------------------------------------

;;; 矩阵相乘
(defun matrix-mul (m1 m2 / m3 m4 m33 x y)
  (setq m3 nil
        m33 nil
  )

  (setq m4 (mtr m2))
  (foreach x m1
    (setq m33 nil)
    (foreach y m4
      (setq m (apply
                '+
                (mapcar
                  '*
                  x
                  y
                )
              )
      )
      (setq m33 (append
                  m33
                  (list m)
                )
      )


    )
    (setq m3 (append
               m3
               (list m33)
             )
    )


  )
  m3
)

(defun mtr (x)
  (apply
    'mapcar
    (cons 'list x)
  )
)


[/php]

aeo

[php]
/*************************************************************************
* 高斯--约当列主元素法求矩阵方程A的逆矩阵，其中A是N*N的矩阵
* 输入: n----方阵A的行数  
*       a----矩阵A
* 输出: det--A的行列式的值
*       a----A的逆矩阵
**************************************************************************/
double gaussian_jodan(int n,double a[N][N])
{  
  int i,j,k,mk;  
  int p[N];  /*记录主行元素在原矩阵中的位置*/  
  double det,m,f;   
  det = 1.0;  
  for(k=0;k<n;k++)  
  {    m=a[k][k];  /*选第K列主元素*/   
       mk=k;   
       for(i=k+1;i<n;i++)      
          if(fabs(m)<fabs(a[k]))      
          {        
             m=a[k];        
             mk=i;      
          }   
          if(fabs(m)<EPS) return(0);   
          if(mk!=k)   
          {      
             for(j=0;j<n;j++) /*将第K列主元素换行到主对角线上*/      
             {        
                f=a[k][j];        
                a[k][j]=a[mk][j];        
                a[mk][j]=f;      
             }      
             p[k]=mk;      
             det = -det;   
          }   
          else      
          p[k]=k;   
          det=det*m;   
          for(j=0;j<n;j++)  /*计算主行元素*/      
          if(j!=k)        
            a[k][j]=a[k][j]/a[k][k];   
          a[k][k]=1.0/a[k][k];   
          for(i=0;i<n;i++) /*消元*/   
          {      
            if(i!=k)      
            {        
              for(j=0;j<n;j++)         
                if(j!=k)            
                   a[j]=a[j]-a[k]*a[k][j];        
              a[k]=-a[k]*a[k][k];      
            }   
           }
  }  
  for(k=n-2;k>=0;k--) /*按主行在原矩阵中的位置交换列*/  
  {   
    if(p[k]!=k)      
    for(i=0;i<n;i++)      
    {        
       f=a[k];        
       a[k]=a[p[k]];        
       a[p[k]]=f;      
    }  
  }  
  return(det);
}   
[/php]

snoopychen

对，aeo版主说的对
本来求解线性方程可以直接用高斯消元法的
我拿逆矩阵来编程本来就有点浪费时间
不过原意是想编逆矩阵的，后来才搞着玩编求方程，所以兜了一个圈
过几天空了，就学着把aeo版主的c代码也编成lisp试试
不过我现在的编法硬是把表改成数组的形式来表达，非常弱，但是又不大知道怎么更好
直接用表的形式来实现这些过程，请各位教一教

楚界

跟vb的一样