Main Content

La traduzione di questa pagina non è aggiornata. Fai clic qui per vedere l'ultima versione in inglese.

Sistemi di equazioni lineari

Considerazioni per i calcoli

Uno dei problemi più importanti nel calcolo tecnico è la risoluzione dei sistemi di equazioni lineari simultanee.

Nella notazione matriciale, questo problema generale assume la seguente forma: Date due matrici A e b, esiste una matrice x unica, per cui Axb o xAb?

È utile considerare un esempio 1x1. Ad esempio, l'equazione

7x = 21

presenta una soluzione unica?

La risposta ovviamente è sì. L'equazione presenta la soluzione unica x = 3. La soluzione viene ottenuta facilmente tramite la divisione:

x = 21/7 = 3.

La soluzione non viene ottenuta ordinariamente calcolando l'inverso di 7, cioè 7–1= 0,142857..., e quindi moltiplicando 7–1 per 21. Sarebbe più laborioso e, rappresentando 7–1 con un numero finito di cifre, meno accurato. Considerazioni simili si applicano anche alle serie di equazioni lineari con più di un'incognita; MATLAB® le risolve senza calcolare l'inverso della matrice.

Benché non si tratti di notazione matematica standard, MATLAB utilizza la familiare terminologia delle divisioni nel caso scalare per descrivere la soluzione di un sistema generale di equazioni simultanee. I due simboli della divisione, la barra, /, e la barra rovesciata, \, corrispondono alle due funzioni di MATLAB mrdivide e mldivide. Questi operatori sono utilizzati per le situazioni in cui la matrice ignota appare sulla sinistra o sulla destra della matrice dei coefficienti:

x = b/A

Denota la soluzione dell'equazione matriciale xA = b, ottenuta con mrdivide.

x = A\b

Denota la soluzione dell'equazione matriciale Ax = b, ottenuta con mldivide.

Si pensi alla "divisione" dei due lati dell'equazione Ax = b o xA = b per A. La matrice di coefficienti A si trova sempre nel "denominatore".

Le condizioni di compatibilità delle dimensioni per x = A\b richiedono che le due matrici A e b presentino lo stesso numero di righe. La soluzione x quindi presenta lo stesso numero di colonne di b e la dimensione delle sue righe è uguale alle dimensioni delle colonne di A. Per x = b/A i ruoli di righe e colonne risultano scambiati.

In pratica, le equazioni lineari della forma Ax = b sono più frequenti rispetto a quelle della forma xA = b. La barra rovesciata è quindi molto più utilizzata rispetto all'altro tipo di barra. La parte restante di questa sezione si concentra sull'operatore barra rovesciata; le proprietà corrispondenti dell'altro tipo di barra possono essere ricavate dall'identità:

(b/A)' = (A'\b').

La matrice dei coefficienti A non deve essere quadrata. Se A presenta le dimensioni mxn, vi sono tre casi:

m = n

Sistema quadrato. Cerca una soluzione esatta.

m > n

Sistema sovradeterminato, con più equazioni che incognite. Trova una soluzione basata sui minimi quadrati.

m < n

Sistema sottodeterminato, con meno equazioni che incognite. Trova una soluzione di base con almeno m componenti non zero.

L’algoritmo mldivide

L'operatore mldivide utilizza risolutori diversi per gestire diversi tipi di matrici dei coefficienti. Viene effettuata la diagnosi automatica dei vari casi tramite l'esame della matrice dei coefficienti. Per maggiori informazioni vedere la sezione "Algoritmi" della pagina di riferimento su mldivide.

Soluzione generale

La soluzione generale di un sistema di equazioni lineari Axb descrive tutte le soluzioni possibili. È possibile trovare la soluzione generale come segue:

  1. Risolvendo il sistema omogeneo corrispondente Ax = 0. A questo scopo utilizzare il comando null, digitando null(A). Questo restituisce una base per lo spazio della soluzione verso Ax = 0. Qualsiasi soluzione è una combinazione lineare di vettori di base.

  2. Trovando una soluzione particolare al sistema non omogeneo Ax = b.

È poi possibile scrivere qualsiasi soluzione per Axb come somma della soluzione particolare verso Ax =b, dal punto 2, oltre una combinazione lineare dei vettori di base dal punto 1.

La parte restante di questa sezione descrive come utilizzare MATLAB per trovare una particolare soluzione verso Ax = b, come nel punto 2.

Sistemi quadrati

La situazione più comune interessa una matrice quadrata dei coefficienti A e un singolo vettore colonna sul lato destro b.

Matrice non singolare dei coefficienti

Se la matrice A è non singolare, allora la soluzione, x = A\b, presenta le stesse dimensioni di b. Ad esempio:

A = pascal(3);
u = [3; 1; 4];
x = A\u

x =
      10
     -12
       5

È possibile confermare che A*x è esattamente pari a u.

Se A e b sono quadrate e delle stesse dimensioni, anche x= A\b presenta le stesse dimensioni:

b = magic(3);
X = A\b

X =
      19    -3    -1
     -17     4    13
       6     0    -6

È possibile confermare che A*x è esattamente uguale a b.

Entrambi gli esempi generano soluzioni esatte e composte da numeri interi. Questo accade perché si è scelta una matrice dei coefficienti pascal(3), che è una matrice a rank pieno (non singolare).

Matrice singolare dei coefficienti

Una matrice quadrata A è singolare se non presenta colonne linearmente indipendenti. Se A è singolare, la soluzione di Ax = b non esiste oppure non è unica. L'operatore barra rovesciata, A\b, genera un avviso se A è quasi singolare o se rileva una singolarità esatta.

Se A è singolare e Ax = b ha una soluzione, è possibile trovare una soluzione particolare non unica digitando

P = pinv(A)*b

pinv(A) è una pseudo-inversa di A. Se Ax = b non presenta una soluzione esatta, pinv(A) restituisce una soluzione basata sui minimi quadrati.

Ad esempio:

A = [ 1     3     7
     -1     4     4
      1    10    18 ]

è singolare, come è possibile verificare digitando

rank(A)

ans =

     2

Poiché A non è a rank pieno, presenta alcuni singoli valori uguali a zero.

Soluzioni esatte. Per b =[5;2;12], l'equazione Ax = b presenta una soluzione esatta, data da

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

Verificare che pinv(A)*b sia una soluzione esatta digitando

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

Soluzioni basate sui minimi quadrati. Se tuttavia b = [3;6;0], Ax = b non presenta una soluzione esatta. In questo caso pinv(A)*b restituisce una soluzione basata sui minimi quadrati. Digitando

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

non si ottiene nuovamente il vettore originale b.

È possibile determinare se Ax = b presenta una soluzione esatta trovando la forma a gradini ridotta della riga della matrice aumentata [A b]. Per eseguire questa operazione sull'esempio, digitare

rref([A b])
ans =
    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

Poiché la riga inferiore contiene tutti zeri, a eccezione dell'ultimo elemento, l'equazione non ha soluzione. In questo caso pinv(A) restituisce una soluzione basata sui minimi quadrati.

Sistemi sovradeterminati

Questo esempio dimostra quanto siano frequenti i sistemi sovradeterminati in diversi tipi di adattamento delle curve a dati sperimentali.

Una quantità y viene misurata in diversi valori di tempo t per produrre le osservazioni sotto riportate. Le seguenti dichiarazioni consentono di inserire i dati e di visualizzarli in una tabella.

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
B = table(t,y)
B=6×2 table
     t      y  
    ___    ____

      0    0.82
    0.3    0.72
    0.8    0.63
    1.1     0.6
    1.6    0.55
    2.3     0.5

Provare a modellare i dati con una funzione di decadimento esponenziale

y(t)=c1+c2e-t.

L'equazione precedente dichiara che il vettore y dovrebbe essere approssimato da una combinazione lineare di altri due vettori. Uno è un vettore costante contenente tutti uno, l'altro è il vettore con i componenti exp(-t). È possibile calcolare i coefficienti incogniti c1 e c2 con un adattamento basato sui minimi quadrati, che minimizza la somma dei quadrati delle deviazioni dei dati dal modello. Vi sono sei equazioni in due incognite, rappresentate da una matrice 6x2.

E = [ones(size(t)) exp(-t)]
E = 6×2

    1.0000    1.0000
    1.0000    0.7408
    1.0000    0.4493
    1.0000    0.3329
    1.0000    0.2019
    1.0000    0.1003

Utilizzare l'operatore barra rovesciata per ottenere la soluzione basata sui minimi quadrati.

c = E\y
c = 2×1

    0.4760
    0.3413

In altri termini, l'adattamento basato sui minimi quadrati dei dati è

y(t)=0.4760+0.3413e-t.

Le dichiarazioni seguenti valutano il modello in corrispondenza di incrementi con spaziamento regolare in t, quindi tracciano il risultato insieme ai dati originali:

T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

E*c non è esattamente uguale a y, ma la differenza potrebbe risultare inferiore rispetto agli errori di misurazione nei dati originali.

Una matrice rettangolare A presenta una carenza di rank se non ha colonne linearmente indipendenti. Se A presenta una carenza di rank, la soluzione basata sui minimi quadrati di AX = B non è unica. A\B genera un errore se A presenta una carenza di rank e genera una soluzione basata sui minimi quadrati. È possibile utilizzare lsqminnorm per trovare la soluzione di X con norma minima tra tutte le soluzioni.

Sistemi sottodeterminati

Questo esempio mostra come la soluzione dei sistemi sottodeterminati non sia unica. I sistemi lineari sottodeterminati presentano più incognite che equazioni. L'operazione di divisione sulla sinistra della matrice in MATLAB trova una soluzione di base basata sui minimi quadrati, che ha al massimo m componenti non zero per una matrice dei coefficienti mxn.

Ecco un esempio casuale:

R = [6 8 7 3; 3 5 4 1]
rng(0);
b = randi(8,2,1)
R =

       6              8              7              3       
       3              5              4              1       


b =

       7       
       8      

Il sistema lineare Rp = b implica due equazioni in quattro incognite. Poiché la matrice dei coefficienti contiene numeri interi di valore ridotto, è appropriato usare il comando format per visualizzare la soluzione in formato razionale. Questa particolare soluzione si ottiene con

format rat
p = R\b
p =

       0       
      17/7     
       0       
     -29/7    

Uno dei componenti non zero è p(2), perché R(:,2) è la colonna di R con la norma più grande. L'altro componente non zero è p(4), perché R(:,4) risulta dominante dopo l'eliminazione di R(:,2).

La soluzione generale completa del sistema sottodeterminato può essere caratterizzata aggiungendo p a una combinazione lineare arbitraria dei vettori spazio nulli, individuabili tramite la funzione null con un'opzione che richiede una base razionale.

Z = null(R,'r')
Z =

      -1/2           -7/6     
      -1/2            1/2     
       1              0       
       0              1       

È possibile confermare che R*Z è zero e che R*x - b residuo è piccolo per qualsiasi vettore x, dove

x = p + Z*q

Poiché le colonne di Z sono i vettori spazio nulli, il prodotto Z*q è una combinazione lineare di questi vettori:

Zq=(x1x2)(uw)=ux1+wx2.

Come esempio illustrativo, scegliere un q arbitrario e costruire x.

q = [-2; 1];
x = p + Z*q;

Calcolare la norma del residuo.

format short
norm(R*x - b)
ans =

   2.6645e-15

Quando vi sono molte soluzioni disponibili, quella con la norma minima risulta di particolare interesse. È possibile usare lsqminnorm per calcolare la soluzione con norma minima basata sui minimi quadrati. Questa soluzione presenta il valore minimo possibile per norm(p).

p = lsqminnorm(R,b)
p =

    -207/137   
     365/137   
      79/137   
    -424/137  

Soluzione di sistemi con più membri sulla destra

Alcuni problemi riguardano la soluzione di sistemi lineari che presentano la stessa matrice dei coefficienti A, ma diversi membri sulla destra b. Quando sono disponibili contemporaneamente diversi valori di b, è possibile costruire b come matrice con diverse colonne e risolvere contemporaneamente tutti i sistemi di equazioni utilizzando un solo comando barra rovesciata: X = A\[b1 b2 b3 …].

A volte, tuttavia, i diversi valori di b non sono tutti disponibili contemporaneamente, il che significa che è necessario risolvere di seguito diversi sistemi di equazioni. Quando si risolve uno di questi sistemi di equazioni con la barra (/) o la barra rovesciata (\), l'operatore fattorizza la matrice dei coefficienti A e utilizza questa decomposizione in matrice per calcolare la soluzione. A ogni successiva risoluzione di un sistema di equazioni simili con diverso b, tuttavia, l'operatore calcola la stessa decomposizione di A, che è un calcolo ridondante.

La soluzione a questo problema consiste nel calcolare preventivamente la decomposizione di A, quindi riutilizzare i fattori per risolvere i diversi valori di b. In pratica, però, il calcolo preventivo della decomposizione in questo modo può risultare difficile, perché è necessario sapere quale decomposizione calcolare (LU, LDL, Cholesky e così via) e anche come moltiplicare i fattori per risolvere il problema. Ad esempio, con la decomposizione LU è necessario risolvere due sistemi lineari per risolvere il sistema originale Ax = b:

[L,U] = lu(A);
x = U \ (L \ b);

Il metodo consigliato per la risoluzione di sistemi lineari con più membri sulla destra consecutivi, invece, consiste nell'utilizzare gli oggetti decomposition. Questi oggetti consentono di sfruttare i vantaggi delle performance del calcolo preventivo della decomposizione in matrice, ma non richiedono di sapere come utilizzare i fattori delle matrici. È possibile sostituire la precedente decomposizione LU con:

dA = decomposition(A,'lu');
x = dA\b;

Se non si è certi del tipo di decomposizione da utilizzare, decomposition(A) sceglie il tipo corretto in base alle proprietà di A, in modo simile a quanto viene effettuato con l'operatore barra retroversa.

Ecco un semplice test dei possibili vantaggi di questo approccio in termini di performance. Il test risolve 100 volte lo stesso sistema lineare sparso utilizzando sia (\) che decomposition.

n = 1e3;
A = sprand(n,n,0.2) + speye(n);
b = ones(n,1);

% Backslash solution
tic
for k = 1:100
    x = A\b;
end
toc
Elapsed time is 9.006156 seconds.
% decomposition solution
tic
dA = decomposition(A);
for k = 1:100
    x = dA\b;
end
toc
Elapsed time is 0.374347 seconds.

Per questo problema, la soluzione decomposition è molto più rapida rispetto all'uso della sola barra retroversa, pur con una sintassi semplice.

Metodi iterativi

Se la matrice dei coefficienti A è grande e sparsa, solitamente i metodi di fattorizzazione non sono efficienti. I metodi iterativi generano una serie di soluzioni approssimative. MATLAB offre diversi metodi iterativi per la gestione delle matrici di input grandi e sparse.

FunzioneDescrizione
pcg

Metodo del gradiente coniugato precondizionato. Questo metodo è appropriato per una matrice A con coefficiente hermitiano definito positivo.

bicg

Metodo del gradiente biconiugato

bicgstab

Metodo stabilizzato del gradiente biconiugato

bicgstabl

Metodo BiCGStab(l)

cgs

Metodo quadrato gradiente coniugato

gmres

Metodo del residuo minimo generalizzato

lsqr

Metodo LSQR

minres

Metodo del residuo minimo. Questo metodo è appropriato per una matrice A con coefficiente hermitiano.

qmr

Metodo del residuo quasi-minimo

symmlq

Metodo LQ simmetrico

tfqmr

Metodo QMR senza trasposizione

Calcolo multithread

MATLAB supporta il calcolo multithread per diverse funzioni numeriche dell'algebra lineare e orientate agli elementi. Queste funzioni vengono eseguite automaticamente su più thread. Per una maggior rapidità di esecuzione di una funzione o espressione su più CPU, occorre soddisfare una serie di condizioni:

  1. La funzione deve eseguire operazioni facilmente ripartibili in sezioni eseguibili contemporaneamente. Queste sezioni devono poter essere eseguite con poche comunicazioni tra i processi. Devono richiedere poche operazioni sequenziali.

  2. Le dimensioni dei dati devono essere sufficientemente grandi, in modo che i vantaggi dell'esecuzione simultanea giustifichino il tempo richiesto per la partizione dei dati e la gestione di thread di esecuzione distinti. Ad esempio, la maggior parte delle funzioni risulta velocizzata solo quando gli array contengono almeno diverse migliaia di elementi.

  3. L'operazione non è vincolata alla memoria; il tempo di accesso alla memoria non deve costituire la maggior parte del tempo di elaborazione. Come regola generale, le funzioni complesse sono più velocizzabili rispetto alle funzioni semplici.

inv, lscov, linsolve e mldivide mostrano un notevole aumento della velocità sugli array a doppia precisione (nell'ordine di almeno 10.000 elementi) quando la modalità multithread è attiva.

Vedi anche

| | | |

Argomenti complementari

Siti web esterni