::-- 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))