Lisp & Array (9): eliminazione Gaussiana

Lisp »Tips 'n Tricks »Lisp & Array »1 »2 »3 »4 »5 »6 »7 »8 »9 »10

ELG implementa un classico metodo per la soluzione di un sistema lineare, cioè determinare i valori delle variabili che soddisfano simultaneamente un insieme di equazioni.
L'idea di base è l'eliminazione, cioè ricavare da una equazione una incognita e sostituirla nelle equazioni restanti. La sostituzione diminuisce la dimensione del problema. Iterando il procedimento si riduce il problema ad un problema ad una dimensione, ad una incognita. Quindi si trovano le altre incognite grazie ad una procedura di sostituzione all'indietro.
L'idea è attribuita a Gauss, ma pare fosse già nota ai cinesi di 2000 anni fa
ELG converte in AutoLISP il seguente codice in C:

eliminate()
{
   int i,j,k,max;
   float t;
   for (i=1; i<=n; i++)
   {
      max = i;
      for (j=i+1; j<=n; j++)
         if (fabs(A[j][i])>fabs(A[max][i])) max = j;
      for (k=i; k<=n+1; k++)
      {   
         t = A[i][k]; A[i][k] = A[max][k]; A[max][k] = t;
      }
      for (j=i+1; j<=n; j++)
         for (k=n+1; k>=i; k--)
            A[j][k] -= A[i][k]*A[j][i]/A[i][i];
   }
}
substitute()
{
   int j,k;
   float t;
   for (j=n; j>=1; j--)
   {
      t = 0.0;
      for (k=j+1; k<=n; k++) t += A[j][k]*X[k];
      X[j] = (A[j][n+1]-t)/A[j][j];
   }
} *

ELG

;|
   ELG.LSP (4 Luglio 2005)
   Copyright (C) 2005 Claudio Piccini.
   All rights reserved
   www.cg-cad.com

   Legge dal file ELG.TXT coefficienti e termini noti
   di un sistema di equazioni.

   Struttura di ELG.TXT:
   '(
   (0  0  0  0  0 0)
   (0  n  n  n  n 0)
   (0  n  n  n  n 0)
   (0  n  n  n  n 0)
   (0  0  0  0  0 0)
   )
   Vincoli:
      1. Ogni lista deve avere lo stesso numero di elementi.
      2. Numero Colonne=Numero Righe+1.
      3. La prima lista deve contenere solo 0
      4. L'ultima lista deve contenere solo 0
      5. Ogni lista deve iniziare con 0 e terminare con 0
      6. Il file deve chiamarsi ELG.TXT
  
   ELG determina i valori delle variabili che soddisfano un sistema
   di equazioni con il metodo della eliminazione Gaussiana (o cinese).

   Un esempio di sistema di 3 equazioni    
   x+3y-4z=8
   x+y-2z=2
   -x-2y+5z=-1
   e relativo file ELG.TXT
   '(
   (0  0  0  0  0 0)
   (0  1  3 -4  8 0)
   (0  1  1 -2  2 0)
   (0 -1 -2  5 -1 0)
   (0  0  0  0  0 0)
   )
|;

(defun LeggeFile ( / nDir nf 
                     i j e 
                     Lst rigo X 
 )
 (setq nDir (getvar "dwgprefix"))
 (setq nf (strcat nDir "elg.txt"))
 (setq Lst (load nf))
 (setq NR (length Lst))
 (setq NC (length (car Lst)))
 (setq i 0)
 (while (< i NR)
  (setq j 0)
  (while (< j NC)
   (setq e (nth j (nth i Lst)))
   (setq rigo (append rigo (list e)))
   (setq j (1+ j))
  )
  (setq X (append X (list rigo))) 
  (setq rigo nil)
  (setq i (1+ i))
 )
 (setq X X)
)

(defun vArray ( ii jj valore / i j e rigo X )
 (setq i 0)
 (while (< i NR)
  (setq j 0)
  (while (< j NC)
   (setq e (nth j (nth i A)))
   (if (and (= j jj)(= i ii))
    (setq rigo (append rigo (list valore)))
    (setq rigo (append rigo (list e))) 
   )
   (setq j (1+ j))
  )
  (setq X (append X (list rigo))) 
  (setq rigo nil)
  (setq i (1+ i))
 )
 (setq X X)
)

(defun vLista ( ii valore / i e rigo XX )
 (setq i 0)
 (while (< i NR)
  (setq e (nth i X))
  (if (= i ii)
   (setq XX (append XX (list valore)))
   (setq XX (append XX (list e))) 
  )
  (setq i (1+ i))
 )
 (setq XX XX)
)

(defun eliminate ( N / i j k 
                       mas e1 e2 
                       val tt 
                       JK IK JI II
 )
 (setq i 1)
 (while (<= i N)
  (setq mas i)
  (setq j (1+ i))
  (while (<= j N)
   (setq e1 (nth i (nth j A)))
   (setq e2 (nth i (nth mas A)))
   (if (> (abs e1)(abs e2))(setq mas j))
   (setq j (1+ j))
  )
  (setq k i)
  (while (<= k (1+ N))
   (setq tt (nth k (nth i A)))
   (setq val (nth k (nth mas A)))
   (setq A (vArray i k val))
   (setq A (vArray mas k tt))
   (setq k (1+ k))
  )
  (setq j (1+ i))
  (while (<= j N)
   (setq k (1+ N))
   (while (>= k i)
    (setq JK (nth k (nth j A)))
    (setq IK (nth k (nth i A)))
    (setq JI (nth i (nth j A)))
    (setq II (nth i (nth i A)))
    (setq val (/ (* IK JI)(* 1.0 II)))
    (setq JK (- JK val))
    (setq A (vArray j k JK))
    (setq k (1- k))
   )
   (setq j (1+ j))
  )
  (setq i (1+ i))
 )
)

(defun substitute ( N / j k tt val
                        JK XK JN JJ
 )
 (setq j N)
 (while (>= j 1)
  (setq tt 0.0)
  (setq k (1+ j))
  (while (<= k N)
   (setq JK (nth k (nth j A)))
   (setq XK (nth k X))
   (setq tt (+ tt (* JK XK)))  
   (setq k (1+ k))
  )
  (setq JN (nth N (nth j A)))
  (setq JJ (nth j (nth j A)))
  (if (/= JJ 0)
   (setq val (/ (- JN tt) (* 1.0 JJ))) 
   (setq val 0)
  )
  (setq X (vLista j val))
  (setq j (1- j))
 )
)

(defun c:elg ( / A X 
                 NR NC 
                 i e str 
 )
 (setvar "cmdecho" 0)
 (setq A (LeggeFile))
 ; inizializza la lista X dei valori
 (setq i 0)
 (while (< i NR)
  (setq X (append X (list 0)))
  (setq i (1+ i))
 )
 (eliminate (1- NR))
 (substitute (1- NR))
 ; legge la lista X dei valori
 (princ "\n")
 (setq i 1)
 (while (< i (1- NR))
  (setq e (nth i X))
  (setq str (strcat "x" (itoa i) "="))
  (princ str)
  (princ e)
  (princ " ")
  (setq i (1+ i))
 )
 (setvar "cmdecho" 1)
 (princ)
)
;;;eof

Test del Lisp

0.005x₁+x₂=0.5
x₁+x₂=1
ELG.TXT:
'(
(0 0 0 0 0)
(0 0.005 1 0.5 0)
(0 1 1 1 0)
(0 0 0 0 0)
)
Command: elg
x1=0.502513 x2=0.497487

x₁+x₂+x₃=1
x₁+x₂+2x₃=2
x₁+2x₂+2x₃=1
ELG.TXT:
'(
(0 0 0 0 0 0)
(0 1 1 1 1 0)
(0 1 1 2 2 0)
(0 1 2 2 1 0)
(0 0 0 0 0 0)
)
Command: elg
x1=1.0 x2=-1.0 x3=1.0

x₁+3x₂-4x₃=8
x₁+x₂-2x₃=2
-x₁-2x₂+5x₃=-1
ELG.TXT:
'(
(0 0 0 0 0 0)
(0 1 3 -4 8 0)
(0 1 1 -2 2 0)
(0 -1 -2 5 -1 0)
(0 0 0 0 0 0)
)
Command: elg
x1=1.0 x2=5.0 x3=2.0

'(
(0 0 0 0 0 0 0)
(0 6 2 0 0 2.25 0)
(0 2 6 1 0 0.45 0)
(0 0 1 8 3 0.15 0)
(0 0 0 3 10 0.075 0)
(0 0 0 0 0 0 0)
)
x1=0.395411 x2=-0.0612342 x3=0.0265823 x4=-0.000474684

'(
(0 0 0 0 0)
(0 6.0 2.175 8.28 0)
(0 2.175 1.387 3.623 0)
(0 0 0 0 0)
)
x1=1.00361 x2=1.03832

Analisi del Lisp

eliminate()
{
   int i,j,k,max;
   float t;
   for (i=1; i<=n; i++)
   {
      max = i;
      for (j=i+1; j<=n; j++)
         if (fabs(A[j][i])>fabs(A[max][i])) max = j;
      for (k=i; k<=n+1; k++)
      {   
         t = A[i][k];
         A[i][k] = A[max][k];
         A[max][k] = t;
      }
      for (j=i+1; j<=n; j++)
         for (k=n+1; k>=i; k--)
            A[j][k] -= A[i][k]*A[j][i]/A[i][i];
   }
}

 [...]
  (setq k i)
  (while (<= k (1+ N))
   (setq tt (nth k (nth i A)))    ; t = A[i][k]
   (setq val (nth k (nth mas A))) ; val = A[max][k]
   (setq A (vArray i k val))      ; A[i][k] = val
   (setq A (vArray mas k tt))     ; A[max][k] = t
   (setq k (1+ k))
  )
  (setq j (1+ i))
  (while (<= j N)
   (setq k (1+ N))
   (while (>= k i)
    (setq JK (nth k (nth j A)))        ; A[j][k]
    (setq IK (nth k (nth i A)))        ; A[i][k]
    (setq JI (nth i (nth j A)))        ; A[j][i]
    (setq II (nth i (nth i A)))        ; A[i][i]
    (setq val (/ (* IK JI)(* 1.0 II))) ; A[i][k]*A[j][i]/A[i][i]
    (setq JK (- JK val))               ; A[j][k]-A[i][k]*A[j][i]/A[i][i]
    (setq A (vArray j k JK))           ; A[j][k] = JK
    (setq k (1- k))
   )
   (setq j (1+ j))
  )
  (setq i (1+ i))
  [...]

substitute()
{
   int j,k;
   float t;
   for (j=n; j>=1; j--)
   {
      t = 0.0;
      for (k=j+1; k<=n; k++) t += A[j][k]*X[k];
      X[j] = (A[j][n+1]-t)/A[j][j];
   }
}


(defun substitute ( N / j k tt val
                        JK XK JN JJ
 )
 (setq j N)
 (while (>= j 1)
  (setq tt 0.0)
  (setq k (1+ j))
  (while (<= k N)
   (setq JK (nth k (nth j A))) ; A[j][k]
   (setq XK (nth k X))         ; X[k]
   (setq tt (+ tt (* JK XK)))  ; t += A[j][k]*X[k]  
   (setq k (1+ k))
  )
  (setq JN (nth N (nth j A))) ; A[j][n+1]
  (setq JJ (nth j (nth j A))) ; A[j][j]
  (if (/= JJ 0)
   (setq val (/ (- JN tt) (* 1.0 JJ))) ; (A[j][n+1]-t)/A[j][j] 
   (setq val 0)
  )
  (setq X (vLista j val)) ; X[j] = val
  (setq j (1- j))
 )
)

Per approfondire:
(*) R. Sedgewick, Algoritmi in C++, 1993, Addison-Wesley

Lisp »Tips 'n Tricks

Ultimo Aggiornamento_Last Update: 4 Luglio 2005