RISOLUZIONE DI UN SISTEMA LINEARE TRAMITE FATTORIZZAZIONE LR

Contenuto

Script per la risoluzione di un sistema lineare denso

%----------------------------------------------------
%  Calcoli di Processo dell'Ingegneria Chimica
%  Davide Manca
%  28-Oct-2003 22:44:32
%   8-Nov-2004 22:50:48
%----------------------------------------------------

%  Inizializzazione dell'ambiente
clear all
clc

Inizializzazione matrice quadrata piena

a = [5. , 3. , 1. ; 2. , -10. , 1. ; 1. , 1. , 7.];
b = [38. , -36. , 64.];

Chiamata alla funzione che fattorizza e risolve il sistema lineare

[xSol iEnd] = DvdSolveA(a,b);

Visualizzazione dei risultati

if iEnd
   n = length(b);
   for i = 1: n
      disp(['x(',num2str(i),') = ',num2str(xSol(i))]);
   end
end

x(1) = 3
x(2) = 5
x(3) = 8

Fattorizzazione LR tramite trasformazione di Gauss di una matrice piena

%----------------------------------------------------
%  Funzione DvdSolveA(A,b)
%  Calcoli di Processo dell'Ingegneria Chimica
%  Davide Manca
%  28-Oct-2003 22:36:44
%   8-Nov-2004 22:50:48
%----------------------------------------------------
%
%  La funzione DvdSolveA determina la soluzione di un sistema quadrato
%  denso tramite fattorizzazione LR e quindi risoluzione del sistema
%  Input: matrice quadrata A e vettore termini noti b
%  Output: vettore soluzione xSol e indice di successo iSuccesso [false o true]
%  N.B.: i dati di input: A e b NON vengono modificati

function [xSol iSuccesso] = DvdSolveA(A,b)

   %  La prima cosa da fare è settare a FALSE l'esito della funzione
   iSuccesso = false;
   nRowsCols = size(A);
   m = length(nRowsCols);

   %  Analisi di consistenza delle dimensioni della matrice A e vettore b
   if m ~= 2
      disp('La matrice non ha dimensione 2. Impossibile continuare !!!')
      return
   end

   if nRowsCols(1) ~= nRowsCols(2)
      disp('La matrice non è quadrata. Impossibile continuare !!!')
      return
   end

   if nRowsCols(1) ~= length(b)
      disp('Dimensioni non consistenti tra matrice coefficienti e vettore termini noti. Impossibile continuare !!!')
      return
   end

   %  Inizio algoritmo risolutivo...
   xSol = b;            %  Per non modificare il vettore dei termini noti...
   a = A;               %  Per non modificare la matrice dei coefficienti...
   n = nRowsCols(1);    %  Per comodità di scrittura...

   %  Fattorizzazione di un sistema quadrato pieno con algoritmo semplificato (trasformazione di Gauss)
   %  Si ottiene la matrice: L * R  = A quindi occorre risolvere il sistema: L (R x) = b
   for k = 1: n-1
       for i = k+1: n
          if a(k,k) == 0
             disp('Errore !!! Impossibile continuare. Pivot NULLO !!!');
             return
          end
          if a(k,k) == 0
             disp(['Pivot nullo. Elemento a(',num2str(k),',',num2str(k),') == 0. Impossibile continuare!!!'])
             return
          end
          a(i,k) = a(i,k) / a(k,k);
          for j = k+1: n
             a(i,j) = a(i,j) - a(i,k)*a(k,j);
          end
       end
   end

   %  Preparazione delle matrici L ed R per soluzione forward e backward
   %  Caricamento matrici R ed L
   R = zeros(n, n);
   for i = 1: n
       for j = i: n
          R(i,j) = a(i,j);
       end
   end

   L = eye(n);    %  La pongo inizialmente uguale ad una matrice identità
   for i = 1: n
       for j = 1: i-1
          L(i,j) = a(i,j);
       end
   end

   %  Risoluzione sistemi triangolari sinitro e destro
   %  Si noti l'utilizzo strategico dell'indice "iSuccesso"...
   [xSol iSuccesso] = DvdSolveL(L,xSol);        %  Risoluzione sistema triangolare sinistro
   if iSuccesso
      [xSol iSuccesso] = DvdSolveR(R,xSol);     %  Risoluzione sistema triangolare destro
   end

Risoluzione forward di un sistema triangolare sinistro

%----------------------------------------------------
%  Funzione DvdSolveL(L,b)
%  Calcoli di Processo dell'Ingegneria Chimica
%  Davide Manca
%  28-Oct-2003 22:36:44
%   8-Nov-2004 22:50:48
%----------------------------------------------------
%
%  La funzione DvdSolveL determina la soluzione di un sistema triangolare
%  sinistro utilizzando l'algoritmo per RIGA
%  Input: matrice triangolare sinistra L e vettore termini noti b
%  Output: vettore soluzione xSol e indice di successo iSuccesso [false o true]
%  N.B.: i dati di input: L e b NON vengono modificati

function [xSol iSuccesso] = DvdSolveL(L,b)

   iSuccesso = false;
   nRowsCols = size(L);
   m = length(nRowsCols);

   if m ~= 2
      disp('La matrice non ha dimensione 2. Impossibile continuare !!!')
      return
   end

   if nRowsCols(1) ~= nRowsCols(2)
      disp('La matrice non è quadrata. Impossibile continuare !!!')
      return
   end

   if nRowsCols(1) ~= length(b)
      disp('Dimensioni non consistenti tra matrice coefficienti e vettore termini noti. Impossibile continuare !!!')
      return
   end

   xSol = b;            %  Per non modificare il vettore dei termini noti...
   n = nRowsCols(1);    %  Per comodità di scrittura...

   %  Risoluzione sistema triangolare sinistro: versione per RIGA
   for i = 1: n
       for k = 1: i-1
          xSol(i) = xSol(i) - L(i,k) * xSol(k);
       end
       if L(i,i) == 0
          disp(['La matrice è singolare. Elemento L(',num2str(i),',',num2str(i),') == 0. Impossibile continuare!!!'])
          return
       end
       xSol(i) = xSol(i) / L(i,i);
   end

   iSuccesso = true;

Risoluzione backward di un sistema triangolare destro

%----------------------------------------------------
%  Funzione DvdSolveR(R,b)
%  Calcoli di Processo dell'Ingegneria Chimica
%  Davide Manca
%  28-Oct-2003 22:08:12
%   8-Nov-2004 22:50:48
%----------------------------------------------------
%
%  La funzione DvdSolveR determina la soluzione di un sistema triangolare
%  destro utilizzando l'algoritmo per RIGA
%  Input: matrice triangolare destra r e vettore termini noti b
%  Output: vettore soluzione xSol e indice di successo iSuccesso [false o true]
%  N.B.: i dati di input: R e b NON vengono modificati

function [xSol iSuccesso] = DvdSolveR(R,b)

   iSuccesso = false;
   nRowsCols = size(R);
   m = length(nRowsCols);

   if m ~= 2
      disp('La matrice non ha dimensione 2. Impossibile continuare !!!')
      return
   end

   if nRowsCols(1) ~= nRowsCols(2)
      disp('La matrice non è quadrata. Impossibile continuare !!!')
      return
   end

   if nRowsCols(1) ~= length(b)
      disp('Dimensioni non consistenti tra matrice coefficienti e vettore termini noti. Impossibile continuare !!!')
      return
   end

   xSol = b;            %  Per non modificare il vettore dei termini noti...
   n = nRowsCols(1);    %  Per comodità di scrittura...

   %  Risoluzione sistema triangolare destro: versione per RIGA
   for i = n: -1: 1
       for k = i+1: n
          xSol(i) = xSol(i) - R(i,k) * xSol(k);
       end
       if R(i,i) == 0
          disp(['La matrice è singolare. Elemento R(',num2str(i),',',num2str(i),') == 0. Impossibile continuare!!!'])
          return
       end
       xSol(i) = xSol(i) / R(i,i);
   end

   iSuccesso = true;