[matlab] Matrice Tridiagonale Fattorizzazione Lu

[atomico82]

ciao a tutti... sono in preda in una crisi psico fisica rivolta a Calcolo Numerico

Mi è stato dato come esercizio di implementare un codice x la risoluzione di un sistema reale AX=b usando la fattorizzazione LU.

La matrice A è una matrice tridiagonale.. purtroppo non trovo appunti in rete che mi risolvano i miei dubbi.

Qualcuno sa come fare? l'ha fatto? o mi può dare consigli? grazie

[lovaz]

Da quanto ne so il fatto che A sia tridiagonale non cambia il procedimento di
decomposizione

[atomico82]

Originariamente inviato da lovaz

Da quanto ne so il fatto che A sia tridiagonale non cambia il procedimento di
decomposizione

in teoria no.. ma mi è stato chiesto di fare un codice + "leggero" memorizzando in 3 vettori le 3 diagonali e di lavorare con quelle.

[stella_650]

Mi è stato dato come esercizio di implementare un codice x la risoluzione di un sistema reale AX=b usando la fattorizzazione LU.La matrice A è una matrice tridiagonale.. purtroppo non trovo appunti in rete che mi risolvano i miei dubbi.

ti serve una fattorizzazione con pivoting oppure no?
ti occorre capirla anche da un punto di vista teorico oppure è una questione esclusiva sul codice?
Chiedo queste cose perchè in base alla tua richiesta, ci sono delle considerazioni diverse da fare...

[atomico82]

mi serve solo in codice matlab.. che poi a capirla ci penso io. Scrivo i ltesto esatto dell'esercizio così è + facile capire.

Progettare ed implementare un codice per la risoluzione di un sistema reale di equazioni lineari Ax=b, dove A è una matrice tridiagonale di dimensione nxn e x e b sono vettori di lunghezza n.
Memorizzare le tre diagonali di A in tre vettori. Risolvere il sistema usando gli
algoritmi di fattorizzazione LU e di forwarde backward substitution, adattati
alla forma tridiagonaledella matrice dei coefficienti.

Grazie

p.s. x quanto riguarda il pivoting immagino di si anche se non sta espressamente scritto.

[atomico82]

allora... penso di aver trovato quello che mi interessa.. purtroppo però è in linguaggio fortram e non riesco a capire bene cosa fa..

Fattorizzazione LU per matrici tridiagonali :

Subroutine lut(d, e, f, n, ifail)

c
c  SCOPO: fattorizzazione di una matrice quadrata tridiagonale A
c         nella forma A=L*U, dove L è una matrice bidiagonale in-
c         feriore a diagonale unitaria e U bidiagonale superiore.
c
c  PARAMETRI DI INPUT:
c                      
c  - d,  vettore di reali,  diagonale principale 
c                           della matrice da fattorizzare
c  - e,  vettore di reali,  diagonale inferiore 
c                           della matrice da fattorizzare
c  - f,  vettore di reali,  diagonale superiore 
c                           della matrice da fattorizzare
c  - n,            intero,  ordine della matrice
c
c  PARAMETRI DI OUTPUT:
c
c  -     d,  vettore di reali,  diagonale principale
c                               della matrice fattorizzata
c  -     e,  vettore di reali,  diagonale inferiore
c                               della matrice da fattorizzata
c  -     f,  vettore di reali,  diagonale superiore
c                               della matrice da fattorizzata
c  - ifail,            intero,  indicatore di errore
c
c  ROUTINE RICHIAMATE:
c
c  - nessuna
c  
c      |a1   c1   0    …   0  |     d = (a1 a2 a3 … an)
c      |b1   a2   c2   …   0  |     e = (b1 b2 … bn-1)
c  A = |0    b2   a3   …   0  |     f = (c1 c2 … cn-1)
c      |…    …    …    …  cn-1|     
c      |0    0    0  bn-1 an  |
c
c      |1   0   .   .   0|       |d1  f1  0   .   .  | 
c      |e1  1   0   .   .|       |0   d2  f2  .   .  | 
c  L = |0   e2  1   .   .|   U = |.   0   d3  .   .  |   
c      |.   .   .   .   .|       |.   .   0   .  fn-1|
c      |0   .   .  en-1 1|       |0   .   .   .   dn |
c
c          |d1    f1       0        .       .        |
c          |e1d1  e1f1+d2  f2       .       .        |
c  A = LU =|0     e2d2     e2f2+f3  f3       .        |
c          |.     0        .      .       fn-1      |
c          |.     .        .      en-1fn-1  en-1fn-1+dn| 
c
c  f(i) = c(i)                                    diagonale superiore
c  d(1) = a(1) 
c  e(i) = b(i)/d(i) => e(i) = e(i)/d(i)           diagonale inferiore
c  d(i+1) = a(i+1)–e(i)f(i) = d(i+1)-e(i)f(i)     diagonale principale
c  i= 1,..,n-1
c

       integer n, ifail
       real    d(n), e(n-1), f(n-1)
      
c  Controllo della dimensione 
       if (n .le. 0) then
          ifail=3
          goto 999
       endif          
                          
c  Fattorizzazione
       do 10 k=1,n-1

c  Controllo della singolarità di un minore principale
c  non viene controllato se d(n) .eq. 0, ma non occorre
c  in quanto d(n) non è utilizzato per calcolare e(n-1),
c  se anche fosse nullo non si effettua la divisione per zero
          if (d(k) .eq. 0) then
             ifail=6
             goto 999
          endif                                    
    
c  Calcolo di un elemento della diagonale inferiore          
          e(k)=e(k)/d(k)

c  Calcolo di un elemento della diagonale principale
          d(k+1)=d(k+1)-e(k)*f(k)

   10  continue   
                 
       ifail=0
       
  999  continue

       End

Soluzione di un sistema tridiagonale (fattorizzato LU) con back e forward substitution :

Subroutine slut(d, e, f, b, n, ifail)

c
c  SCOPO: risoluzione di un sistema lineare tridiagonale A*x=b me-
c         diante la risoluzione del sistema equivalente (L*U)*x=b,
c         dove L*U è la fattorizzazione di A.
c
c  PARAMETRI DI INPUT:
c
c  - d,  vettore di reali,  diagonale principale di U
c  - e,  vettore di reali,  diagonale inferiore di L
c  - f,  vettore di reali,  diagonale superiore di U 
c  - b,  vettore di reali,  colonna dei termini noti
c  - n,            intero,  ordine del sistema
c
c  PARAMETRI DI OUTPUT:
c
c  -     b,  vettore di reali,  soluzione del sistema
c  - ifail,            intero,  indicatore di errore
c
c  ROUTINE RICHIAMATE:
c  
c  - nessuna.
c
c      |1   0   .   .   0|   |b1|   |b1|
c      |e1  1   0   .   .|   |b2|   |b2|
c  L = |0   e2  1   .   .| * |b2| = |b3|
c      |.   .   .   .   .|   |..|   |..|
c      |0   .   .  en-1 1|   |bn|   |bn|
c
c      |d1  f1  0   .  .  |    |b1|   |b1|
c      |0   d2  f2  .  .  |    |b2|   |b2|
c  U = |.   0   d3  .  .  |  * |b3| = |b3|
c      |.   .   0   . fn-1|    |..|   |..|
c      |0   .   .   .  dn |    |bn|   |bn|
c

       integer n, ifail
       real    d(n), e(2:n), f(n-1), b(n)
       
c     Controllo della dimensione
       if (n .le. 0) then
          ifail=3
          goto 999
       endif
       
c     Controllo della singolarità di U
       do 10 i=1,n
          if (d(i) .eq. 0) then
             ifail=4
             goto 999
          endif
   10  continue
   
c     Risoluzione di L*y=b (con y=U*x)
c     forward substirution v. colonna
       do 20 i=2,n
          b(i)=b(i)-e(i)*b(i-1)
   20  continue
       
c     Risoluzione di U*x=y 
c     backward substitution v. colonna
       b(n)=b(n)/d(n)
       do 30 i=n-1,1,-1
             b(i)=(b(i)-f(i)*b(i+1))/d(i)
   30  continue

       ifail=0
  
  999  continue
      
      End

chi me lo traduce in matlab?