L'algoritmo di eliminazione di Gauss

Fortran di Giuseppe Ciaburro guida dal 04-06-2002

Popolarità:

economia news e media viaggi informatica internet salute e benessere int rattenimento e spettacolo sport tempo libero istruzio ne e formazione arte cultura scienza

Le categorie della guida

I link preferiti della guida

Il Sondaggio

Quale versione del fortran utilizzi?

Guarda i risultati

Algoritmi

A cura di Giuseppe Ciaburro

Pubblicato il 24/08/2001

» Invia tramite EMAIL

» Versione per la STAMPA

» Le vostre opinioni

Algoritmo di eliminazione di Gauss per la risoluzione di in sistema di equazioni lineare.

Il metodo di eleiminazione di Gauss si basa sulla trasformazione di un generico sistema non singolare in uno triangolare equivalente, avente cioé la stessa soluzione.
Se ad una equazione del sistema si sostituisce una combinazione lineare dell'equazione con un'altra dello stesso sistema, il nuovo risulta equivalente al precedente.
Il sistema sotto esame è del tipo e la matrice A è tridiagonale, e ciò semplifica il metodo di Gauss: dato il seguente sistema

$\begin{displaymath} \begin{array}{*{10}{c}} a_{11}x_1 &+& a_{12}x_2 & & & & & & ... ... & & & & a_{nn-1}x_{n-1} &+& a_{nn}x_n & =b_n \\ \end{array} \end{displaymath}$

eliminiamo l'incognita

dalla seconda equazione sottraendo da quest'ultima la prima equazione moltiplicata per

$\begin{displaymath} m_1=\frac{a_{21}}{a_{11}} \end{displaymath}$

ottenendo

$\begin{displaymath} \begin{array}{*{10}{c}} a_{11}x_1 &+& a_{12}x_2 & & & & & & ... ... & & a'_{22}x_2 &+& a_{23}x_3 & & & & & =b'_2 \\ \end{array} \end{displaymath}$

Il coefficiente $a_{23}$ rimane invariato perché $a_{13}=0$ e

.
Il sistema viene triangolarizzato in

passi: nel nostro caso nel termine

il numeratore è costante e pari a $-\gamma$ . La diagonale superiore rimane invariata e quindi per memorizzarla basta uno scalare float. La diagonale principale e la successione dei moltiplicatori

(evidenziati in neretto) sono memorizzate in una matrice $n \times 2$

$\begin{displaymath} \begin{array}{*{10}{c}} {\bf a_{11}} & - \gamma & & &\quad &... ...\ & &{\bf a'_{11}} & - \gamma & & {\bf m_2} \\ \end{array} \end{displaymath}$

secondo le seguenti linee di codice

$\displaystyle {\tt v(i,0)}$	$\textstyle {\tt =}$	$\displaystyle {\tt v(i,0) - V_{ds} * m_i }$
$\displaystyle {\tt v(i,1)}$	$\textstyle {\tt =}$	$\displaystyle {\tt m_i }$

con $i=1,2,\cdots n-1$ e

$\begin{displaymath} V_{ds} = -\gamma \quad {\rm e} \quad m_i= \frac{-\gamma}{v(i-1,0)} \end{displaymath}$

Con la seguente istruzione si opera la trasformazione del termine noto

$\begin{displaymath} {\tt b(i)= b(i) - b(i-1) * m_i } \end{displaymath}$

mentre il calcolo del vettore incognito viene effettuato grazie ad un ciclo iterativo, fatta eccezione per $x_{n-1}$

$\displaystyle {\tt x(n-1)}$	$\textstyle =$	$\displaystyle {\tt b(n-1)/v(n-1,0)}$
$\displaystyle {\tt x(k) }$	$\textstyle =$	$\displaystyle {\tt ( b(k)-V_{ds}*x(k+1))/v(k,0}$

con $k=n-2,n-3, \cdots 0$ .

Di seguito si riporta il listato del programma base per la risoluzione di un sistema tridiagonale.
Per semplicità il sistema è di ordine 3.


program pgauess
!Programma in fortran 90 che effettua
!eliminazione di Gauss

      integer, parameter :: n=4, ia=4
      real, dimension (ia,n):: a 
      real, dimension (n)::s,b,x
      integer, dimension(n)::l
interface
      subroutine gauss(n,a,ia,l,s)
      integer, intent(in) :: n,ia
      real, dimension (:,:), intent(inout)::a
      real, dimension (n)::s
      integer, dimension(n) , intent(out)::l
      end subroutine gauss
      subroutine solve(n,a,ia,l,b,x)
      integer, intent(in)::n,ia
      real, dimension(:,:), intent(in)::a
      real, dimension(:), intent(in) :: b
      integer, dimension(n), intent(in)::l
      real, dimension(n), intent(out) :: x
      end subroutine solve
end interface

      a(1,1)=3.0; a(1,2)=-13.0; a(1,3)=9.0; a(1,4)=3.0
      a(2,1)=-6.0; a(2,2)=4.0; a(2,3)=1.0; a(2,4)=-18.0
      a(3,1)=6.0; a(3,2)=-2.0; a(3,3)=2.0; a(3,4)=4.0
      a(4,1)=12.0; a(4,2)=-8.0; a(4,3)=6.0; a(4,4)=10.0
      print*, "the coefficient matrix a"
      b(1)=-19.0; b(2)=-34.0; b(3)=16.0; b(4)=26.0
      do i=1,n
      print*, (a(i,j),j=1,n)
      end do
      print*, "the following are the constants b(i)"
      print*, (b(i), i=1,n)
      call  gauss(n,a,n,l,s)
      print*, "the matrix a after gauss"
      do i=1,n
         print*, (a (i,j), j=1,n)
      end do
      call solve(n,a,ia,l,b,x)
      print*, "the solution x(i) ", (x(i), i=1,n)
end program pgauess
  
subroutine gauss(n,a,ia,l,s)    
      integer, intent(in) :: n, ia
      real, dimension (:,:), intent(inout):: a
      real, dimension (n):: s
      integer, dimension(n), intent(out)::l
      do  i = 1,n
        l(i) = i  
        smax = 0.0
        do j = 1,n
          smax = amax1(smax,abs(a(i,j)))
        end do
        s(i) = smax 
    end do 
     do  k = 1,n-1
        rmax = -1.0
        do i = k,n
          r = abs(a(l(i),k))/s(l(i))  
          if(r <= rmax)  exit    
          j = i   
          rmax = r
         end do 
  
        lk = l(j) 
        l(j) = l(k) 
        l(k) = lk  
        do i = k+1,n      
          xmult = a(l(i),k)/a(lk,k)   
          do j = k+1,n    
            a(l(i),j) = a(l(i),j) - xmult*a(lk,j) 
          end do 
          a(l(i),k) = xmult 
        end do 
     end do   
end subroutine gauss 
  
subroutine solve(n,a,ia,l,b,x)  
      integer, intent(in) :: n, ia
      real, dimension (:,:), intent(in):: a
      real, dimension (:), intent(in) :: b
      real, dimension (n), intent(out):: x
      integer, dimension(n), intent(in)::l
      integer :: k, i, j
      do k = 1,n-1
         do i = k+1,n      
          b(l(i)) = b(l(i)) - a(l(i),k)*b(l(k)) 
         end do 
     end do 
      x(n) = b(l(n))/a(l(n),n)
     do i = n-1,1,-1     
        sum = b(l(i))       
        do j = i+1,n      
          sum = sum - a(l(i),j)*x(j)  
        end do 
        x(i) = sum/a(l(i),i)
     end do 
end subroutine solve

Feed