::-- hoxide [2005-10-02 02:51:44]

1. LISP

2. 资源

3. 代码合集

3.1. 数值分析

3.1.1. 基础工具集 szcom.lisp

#!lisp
;; szcom.lisp

(provide 'szcom)

(defpackage "szcom" 
  (:use "COMMON-LISP" "EXT")
  (:export "RMAT" "PMAT" "PDMAT" "PSMAT" ) 
  (:nicknames "SZCOM")
)

(in-package "szcom")

(defun RMAT (&optional input-stream) 
  (prog
   ((m (read input-stream))
    (n (read input-stream))
    tn
    )
   (if (not (and (integerp m)
                 (integerp n)))
       (error "RMAT: m or n must be intergers.~%") 
     )
   (setq A (make-array (list m n) :initial-element 0))
   (loop for i from 0 below m
         do
         (loop for j from 0 below n
               do
               (setq tn (read input-stream))
               (cond ((rationalp tn)
                      (setf (aref A i j) tn))
                     ((null tn) (error "RMAT: unexpected steam END.~%"))
                     (T (error "RMAT: the member of Array must be numbers~%"))
                     )
               )
         )
   (return A)
   )
  )
         

(defun PMAT (A)
  (prog
   ((m (car (array-dimensions A)))
    (n (cadr (array-dimensions A))))
   (format t "~%")
   (loop for i from 0 to (- m 1)
         do 
         (loop for j from 0 to (- n 1)
               do
               (format t "~12,6G" (aref A i j))
               )
         (format t "~%" )
         )
   )
  )

(defun PSMat (A)
  (prog
   ((n (decf (car (array-dimensions A)))))
   (loop for i from 0 to n
         do 
         (loop for j from 0 to n
               do
               (format t "~12,G~^"
                       (let ((ind (+ (- j i) 1)))
                         (if (and (<= 0 ind) (< ind 3))
                             (aref A i ind)
                           0
                           )
                         )
                       )
               )
         (format t "~%")
         )
   )
  )

(defun PDMAT (A)
  (prog
   ((n (car (array-dimensions A))))
   (format t "~%")
   (loop for i from 0 to (- n 1)
         do 
         (loop for j from 0 to (- n 1)
               do
               (format t "~12,6G~^" (if (= i  j) (aref A i) 0)
               )
               )
         (format t "~%" )
         )
   )
  )

3.1.2. LU 分解法解线性方程组

#!lisp
;; LU.lisp
(load "szcom")
(import '(szcom:PMAT szcom:RMAT))

(defun LU-S (A LU)
  (prog
   ((n  (car (array-dimensions A)) ))
   (loop for k from 0 to n
         do
         (loop for j from k  below n
               do
               (setf (aref LU k j) 
                     (- (aref A k j) 
                        (loop for q from 0 below k
                              sum (* (aref LU k q) (aref LU q j))
                              )
                        )
                     )
               )
         (loop for i from (+ k 1) below n
               do
               (setf (aref LU i k)
                     (/ 
                      (- (aref A i k) 
                         (loop for q from 0 below k
                               sum (* (aref LU i q) (aref LU q k))
                               )
                         )
                      (aref LU k k)
                      )
                     )
               )
         )
   )
  )

(defun LU-C (LU B X)
  (prog
    ((n  (car (array-dimensions LU))))
    (loop for k from 0 below n
          do
          (setf (aref X k 0)
                ( - (aref B k 0)
                    (loop for i from 0 below k
                          sum (* (aref X i 0) (aref LU k i)))
                    )
                )
          )
    (loop for k from (- n 1) downto 0
          do
          (setf (aref X k 0)
                (/
                 (- (aref X k 0)
                    (loop for j from (+ k 1) below n
                          sum (* (aref X j 0) (aref LU k j))
                          )
                    )
                 (aref LU k k)
                 )
                )
          )
    )
  )

(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PMat A)
(pprint A)
(format t "~% B: ~%")
(PMat B)
(pprint B)

(format t "~%Result of LU: ~%")
;(setq LU (make-array (array-dimensions A)))
(setq LU A)
(LU-S A LU)
(PMat LU)
(format t "~%In Pairs: ~%" )
(pprint LU)

(format t "~%Result of LUS: ~%")
;(setq X (make-array (array-dimensions B)))
(setq X B)
(LU-C LU B X)
(PMat X)
(pprint X)

执行结果:

bash-3.00$ cat ex2.dat
4 4
1  4 -2 3
2  2  0 4
3  0 -1 2
1  2  2 -3

4 1
6
2
1
8

bash-3.00$ ./ex2.lisp < ex2.dat

 A :

 1.00000     4.00000    -2.00000     3.00000
 2.00000     2.00000    0.000000     4.00000
 3.00000    0.000000    -1.00000     2.00000
 1.00000     2.00000     2.00000    -3.00000

#2A((1 4 -2 3) (2 2 0 4) (3 0 -1 2) (1 2 2 -3))
 B:

 6.00000
 2.00000
 1.00000
 8.00000

#2A((6) (2) (1) (8))
Result of LU:

 1.00000     4.00000    -2.00000     3.00000
 2.00000    -6.00000     4.00000    -2.00000
 3.00000     2.00000    -3.00000    -3.00000
 1.00000    0.333333    -.888889    -8.00000

In Pairs:

#2A((1 4 -2 3) (2 -6 4 -2) (3 2 -3 -3) (1 1/3 -8/9 -8))
Result of LUS:

 1.00000
 2.00000
0.000000
-1.00000

#2A((1) (2) (0) (-1))

3.1.3. 三对角阵的追赶法

#!lisp

;; TLU.lisp

(load "szcom")
(import '(szcom:PMAT szcom:PSMAT szcom:RMAT))


(defun shooting (A LU)
  (prog
   ((n  (car (array-dimensions A))))
   (loop for k from 0 below n
         do
         (setf (aref LU k 0) (aref A k 0))
         (setf (aref LU k 1)
               ( - (aref A k 1) 
                   (if (> k 0)
                       (* (aref LU k 0)
                          (aref LU (- k 1) 2)
                          )
                     0)
                   )
               )
         (if (< k n)
             (setf (aref LU k 2)
                   (/ (aref A k 2)
                      (aref LU k 1)
                      )
                   )
           )
         )
   )
  )

(defun shoot (LU B X)
  (prog
   ((n  (car (array-dimensions LU))))
   (loop for k from 0 below n
         do
         (setf (aref X k 0)
               (/
                (- (aref B k 0)
                   (if (> k 0)  
                       (* (aref LU k 0) (aref X (- k 1) 0))
                     0 
                     )
                   )
                (aref LU k 1)
                )
               )
         )
   (pprint X)
   (loop for k from (- n 1) downto 0
         do
         (decf (aref X k 0)
           (if (< k (- n 1)) 
               (* (aref LU k 2) (aref X (+ k 1) 0))
             0)
           )
         )
   )
  )


(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PSMat A)
(pprint A)
(format t "~% B: ~%")
(PMat B)
(pprint B)


(format t "~%Result of shooting: ~%")
(setq LU (make-array (array-dimensions A)))
(shooting A LU)
(PSMat LU)
(format t "~%In Pairs: ~%" )
(pprint LU)

(format t "~%Result of shooting result: ~%")
(setq X (make-array (array-dimensions B)))
(shoot LU B X)
(PMat X)
(pprint X)

例子:

bash-3.00$ cat ex3.dat
5 3
0  13   12
11  10   9
8   7    6
5   4    3
2   1    0
5 1
3
0
-2
6
8

bash-3.00$ ./ex3.lisp < ex3.dat

 A :
     13.         12.         0.0         0.0         0.0
     11.         10.          9.         0.0         0.0
     0.0          8.          7.          6.         0.0
     0.0         0.0          5.          4.          3.
     0.0         0.0         0.0          2.          1.

#2A((0 13 12) (11 10 9) (8 7 6) (5 4 3) (2 1 0))
 B:

 3.00000
0.000000
-2.00000
 6.00000
 8.00000

#2A((3) (0) (-2) (6) (8))
Result of shooting:
     13.    .9230769         0.0         0.0         0.0
     11.    -.15384616       -58.5         0.0         0.0
     0.0          8.        475.    1.263157900E-2     0.0
     0.0         0.0          5.    3.9368422    .7620321
     0.0         0.0         0.0          2.    -.5240642

In Pairs:

#2A((0 13 12/13)
    (11 -2/13 -117/2)
    (8 475 6/475)
    (5 374/95 285/374)
    (2 -98/187 0))
Result of shooting result:

#2A((3/13) (33/2) (-134/475) (32/17) (-396/49))
 5.71837
-5.94490
-.383673
 8.04082
-8.08163

#2A((1401/245) (-2913/490) (-94/245) (394/49) (-396/49))

3.1.4. 对称阵的LLT 分解法

#!lisp
;; LLT.lisp

(load "szcom")
(import '(szcom:RMAT szcom:PMAT))

(defun LLT (A)
  (prog
   ((n (car (array-dimensions A))))
   (loop for k from 0 below n
         do
         (setf (aref A k k)
               (sqrt
                (- (aref A k k) 
                   (loop for j from 0 to (- k 1)
                         sum (* (aref A k j)
                                (aref A j k)
                                )
                         )
                   )
                )
               )
         (loop for i from (+ k 1) below n
               do
               (setf (aref A i k)
                     (/ (- (aref A i k)
                           (loop for j from 0 to (- k 1)
                                 sum (* (aref A i j)
                                        (aref A k j)
                                        )
                                 )
                           )
                        (aref A k k)
                        )
                     )
               (setf (aref A k i ) (aref A i k))
               )
         )
   )
  )

(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PMat A)
(pprint A)
(pprint B)


(format t "~%Result of LLT: ~%")
(LLT A)
(PMat A)
(format t "~%In Pairs: ~%" )
(pprint A)



例子

bash-3.00$ cat LLT.dat
4 4
10     7      8    7
 7     5      6    5
 8     6     10    9
 7     5      9   10
4 1
 1
 2
 3
 4bash-3.00$ ./LLT.lisp < LLT.dat

 A :

 10.0000     7.00000     8.00000     7.00000
 7.00000     5.00000     6.00000     5.00000
 8.00000     6.00000     10.0000     9.00000
 7.00000     5.00000     9.00000     10.0000

#2A((10 7 8 7) (7 5 6 5) (8 6 10 9) (7 5 9 10))
#2A((1) (2) (3) (4))
Result of LLT:

 3.16228     2.21359     2.52982     2.21359
 2.21359    0.316228     1.26491    0.316228
 2.52982     1.26491     1.41421     2.12132
 2.21359    0.316228     2.12132    0.707108

In Pairs:

#2A((3.1622777 2.2135944 2.529822 2.2135944)
    (2.2135944 0.3162276 1.2649105 0.3162276)
    (2.529822 1.2649105 1.414214 2.1213198)
    (2.2135944 0.3162276 2.1213198 0.70710814))

3.1.5. 对称阵的LDLT分解法

#!lisp
;; LDLT.lisp

(load "szcom")
(import '(szcom:RMAT szcom:PMAT szcom:PDMAT))

(defun LDLT (A)
  (prog
   ((n (car (array-dimensions A))))
   (setq L (make-array (list n n) :initial-element 0)
         D (make-array (list n) :initial-element 0))
   (loop for k from 0 below n
         do
         (setf (aref D k)
               (- (aref A k k)
                  (loop for j from 0 to (- k 1)
                        sum (* (aref L k j)
                               (aref L k j)
                               (aref D j)
                               )
                        )
                  )
               )
         (setf (aref L k k) 1)
         (loop for i from (+ k 1) below n
               do
               (setf (aref L i k)
                     (/ (- (aref A i k)
                           (loop for j from 0 to (- k 1)
                                sum (* (aref L i j)
                                       (aref L k j)
                                       (aref D j)
                                       )
                                ))
                        (aref D k)
                        )
                     )
               )
         )
   (return (cons L  D))
   )
  )

(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PMat A)
(pprint A)
(pprint B)

(format t "~%Result of LDLT: ~%")
(setq res (LDLT A))
(PMat (car res))
(PDMAT (cdr res))
(format t "~%In Pairs: ~%" )
(pprint res)

例子:

bash-3.00$ cat LLT.dat
4 4
10     7      8    7
 7     5      6    5
 8     6     10    9
 7     5      9   10
4 1
 1
 2
 3
 4

bash-3.00$ ./LDLT.lisp < LLT.dat

 A :

 10.0000     7.00000     8.00000     7.00000
 7.00000     5.00000     6.00000     5.00000
 8.00000     6.00000     10.0000     9.00000
 7.00000     5.00000     9.00000     10.0000

#2A((10 7 8 7) (7 5 6 5) (8 6 10 9) (7 5 9 10))
#2A((1) (2) (3) (4))
Result of LDLT:

 1.00000    0.000000    0.000000    0.000000
0.700000     1.00000    0.000000    0.000000
0.800000     4.00000     1.00000    0.000000
0.700000     1.00000     1.50000     1.00000

 10.0000    0.000000    0.000000    0.000000
0.000000    0.100000    0.000000    0.000000
0.000000    0.000000     2.00000    0.000000
0.000000    0.000000    0.000000    0.500000

In Pairs:

(#2A((1 0 0 0) (7/10 1 0 0) (4/5 4 1 0) (7/10 1 3/2 1)) . #(10 1/10 2 1/

3.1.6. SOR 迭代法解线性方程组

#!lisp
;; SOR.lisp

(load "szcom")
(import '(szcom:PMAT szcom:RMAT))

(defun SOR (A B X &key ((:weight w) 1.0) ((:ep ep) 1.e-8 ))
  (prog ((n (car (array-dimensions B)))
         (k 0) (e ep) tx)
        (loop until (< e ep)
              do
              (setf e 0)
              (loop for k from 0 below n
                    do
                    (setf tx 
                          (/ 
                           (- (aref B k 0)
                              (loop for j from 0 below n
                                    sum
                                    (if (= j k) 0
                                      (* (aref A k j)
                                         (aref X j 0) )
                                      )))
                           (aref A k k)
                           )
                          )
                    (setf tx (+ (* w tx) (* (- 1 w) (aref X k 0)) ))
                    (setf e (max e (abs (- (aref X k 0) tx)))) 
                    (setf (aref X k 0) tx)
                    )
              (incf k)
              ;(pprint X)
              )
        (return k)
        )
  )

(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PMat A)
(pprint A)
(format t "~% B: ~%")
(PMat B)
(pprint B)


(format t "~%Result of SOR result: ~%")
(setq X (make-array (array-dimensions B) :initial-element 0))
(setq k (SOR A B X :weight 1.3 :ep 1e-5))
(format t "~% Loop ~D time(s)" k)
(PMat X)
(pprint X)

例子:

bash-3.00$ cat SOR.dat
4 4
-4  1  1  1
 1 -4  1  1
 1  1 -4  1
 1  1  1 -4
4 1
1
1
1
1

bash-3.00$ ./SOR.lisp < SOR.dat

 A :

-4.00000     1.00000     1.00000     1.00000
 1.00000    -4.00000     1.00000     1.00000
 1.00000     1.00000    -4.00000     1.00000
 1.00000     1.00000     1.00000    -4.00000

#2A((-4 1 1 1) (1 -4 1 1) (1 1 -4 1) (1 1 1 -4))
 B:

 1.00000
 1.00000
 1.00000
 1.00000

#2A((1) (1) (1) (1))
Result of SOR result:

 Loop 12 time(s)
-1.00000
-.999999
-1.00000
-1.00000

#2A((-1.0000014) (-0.99999917) (-1.0) (-1.0000004))

Lisp (last edited 2005-12-16 11:33:37 by ZoomQuiet)