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 Ax= b o xA= b?
È 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:
| Denota la soluzione dell'equazione matriciale xA = b, ottenuta con |
| Denota la soluzione dell'equazione matriciale Ax = b, ottenuta con |
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 Ax= b descrive tutte le soluzioni possibili. È possibile trovare la soluzione generale come segue:
Risolvendo il sistema omogeneo corrispondente Ax = 0. A questo scopo utilizzare il comando
null
, digitandonull(A)
. Questo restituisce una base per lo spazio della soluzione verso Ax = 0. Qualsiasi soluzione è una combinazione lineare di vettori di base.Trovando una soluzione particolare al sistema non omogeneo Ax = b.
È poi possibile scrivere qualsiasi soluzione per Ax = b 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
.
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 e 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 è
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')
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 m
xn
.
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:
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.
Funzione | Descrizione |
---|---|
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:
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.
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.
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
mldivide
| mrdivide
| pinv
| decomposition
| lsqminnorm