Comandi essenziali di MATLAB
Comandi base
MATLAB può essere usato come una calcolatrice:
>> 1+sin(3)*exp(-2)
ans =
1.0191Le operazioni standar +, -, *, / funzionano come uno si aspetta. Ci sono più o meno tutte le funzioni standard sin(x), cos(x), exp(x) e cosi via.
La cosa più importante è che ogni cosa in MATLAB è una matrice. Ad esempio
>> a = 1
a =
1
>> size(a)
ans =
1 1cioè la variabile a è assegnata allo scalare 1 ma è trattata come una matrice 1\times 1 come il comando size fa vedere. Il comando size infatti restituisce una matrice 1\times 2 dove il primo elemento è il numero di righe e il secondo numero è il numero di colonne.
Inizializzare una matrice in MATLAB è molto facile infatti
>> a = [ 1 2 3 ]
a =
1 2 3
>> size(a)
ans =
1 3
>> b = [ 1 2 3; 4 5 6 ]
b =
1 2 3
4 5 6
>> size(b)
ans =
2 3come si vede dall’esempio, una matrice è definita elencando gli elementi tra parentesi quadre. Gli elementi sono elencati per righe e separati da spazi (o virgole opzionali). Le righe sono separate da ;.
Nel caso di matrici con molte righe può essere utile spezzare il comando su più righe con il comando di continuazione ..., ad esempio
>> a = [ 1 2 3; ...
4 5 6; ...
7 8 9 ]
a =
1 2 3
4 5 6
7 8 9
>> size(a)
ans =
3 3Accedere agli elementi di una matrice
Accedere agli elementi di una matrice è piuttosto facile e intuitivo, ad esempio
>> a = [ 1 2 3; ...
4 5 6; ...
7 8 9 ]
a =
1 2 3
4 5 6
7 8 9
>> a(2,3)
ans =
6
>> a(2,1)
ans =
4
>> a(1)
ans =
1
>> a(2)
ans =
4
>> a(3)
ans =
7
>> a(4)
ans =
2
>> a(5)
ans =
5notate che ci sono due modi di accedere ad un elemento di una matrice,
con due indici
a(i,j)doveiè l’indice di riga ejè l’indice di colonna dell’elemento.con un indice
a(i)doveiè l’indice dell’elementoiesimo della matrice srotolata per colonne.
accedere con un solo indice è comodo nel caso di matrici 1\times n o n \times 1 perché vengono trattate di fatto come dei vettori (di fatto lo sono).
Accedere a sottoblocchi di una matrice
Accedere a sottoblocchi di matrici è altrettanto facile che accedere a degli elementi della stessa. Se al posto degli indici si mettono dei vettori riga allora è possibile estrarre blocchi di righe e colonne. Se al posto degli indici si mette un vettore allora si estraggono gli elementi elencati.
>> a = [ 1 2 3 1 1 1; ...
4 5 6 2 2 2; ...
7 8 9 3 2 1 ]
a =
1 2 3 1 1 1
4 5 6 2 2 2
7 8 9 3 2 1
>> a( [1,3], [1,2,5] )
ans =
1 2 1
7 8 2
>> a( [1,3,5,8,12] )
ans =
1 7 5 6 3Questo modo di operare è molto utile nella scrittura di algoritmi che operano su matrici a blocchi o su istruzioni condizionali che vedremi più avanti.
operazioni sulle matrici
Per le matrici gli operatori +, -, * funzionano come uno si aspetta:
>> a = [ 1 2 3; ...
4 5 6; ...
7 8 9 ]
a =
1 2 3
4 5 6
7 8 9
>> b = [ 1 0 1; ...
2 1 1; ...
1 0 2 ]
b =
1 0 1
2 1 1
1 0 2
>> a+b
ans =
2 2 4
6 6 7
8 8 11
>> a-b
ans =
0 2 2
2 4 5
6 8 7
>> a*b
ans =
8 2 9
20 5 21
32 8 33Un commento a parte va fatto per gli operatori / e \. Infatti a/b è formalmente equivalente a a*inv(b) dove inv(b) è l’operatore che costruisce la matrice inversa di b. In pratica a/b non costruisce l’inversa di b ma risolve il problema a=x*b. Invecie a\b è formalmente equivalente a inv(a)*b. In pratica a\b non costruisce l’inversa di b ma risolve il problema a*x=b.
>> a = [ 4 1 1; ...
1 4 1; ...
1 2 4 ]
a =
4 1 1
1 4 1
1 2 4
>> b = [ 3 0 1; ...
2 3 1; ...
1 1 3 ]
b =
3 0 1
2 3 1
1 1 3
>> norm( a/b - a*inv(b) )
ans =
1.2776e-16
norm( a\(\backslash\)b - inv(a)*b )
ans =
1.5348e-16Il comando norm(a) è la norma L^2 della matrice a che vale
\|\bm{A}\| = \sqrt{\sigma(\bm{A}^T\bm{A})}
dove \sigma(\bm{B}) è il raggio spettrale di \bm{B} cioè il massimo del modulo degli autovalori di \bm{B}.
Un’altra operazione importante importante sulle matrici è la trasposizione fatta con l’apice:
>> a = [ 4 2 2; ...
1 4 1; ...
1 2 4 ]
a =
4 2 2
1 4 1
1 2 4
>> a'
ans =
4 1 1
2 4 2
2 1 4operazioni elemento per elemento
Un’altra caratteristica molto comoda di MATLAB è la possibilità di operare elemento per elemento su una matrice in modo elementare. GLi operatori elemento per elemento hanno il prefisso punto .. Ad esempio
a = [ 4 2 2; ...
1 2 4 ];
b = [ 1 2 3; ...
2 3 3 ];
a .* b % prodotto elemento per elemento
ans =
4 4 6
2 6 12
a ./ b % divisione elemento per elemento
ans =
4.0000 1.0000 0.6667
0.5000 0.6667 1.3333
a .^ b % elevamento a potenza elemento per elemento
ans =
4 4 8
1 8 64Nell’esempio ^ è l’operatore elevamento a potenza, inoltre si nota che ponendo alla fine di una istruzione ; si ha la soppressione dell’output. Cioè normalmente il risultato di ogni istruzione viene stampato subito dopo, il ; a fine istruzione sopprime la stampa. Attenzione che gli operatori funzionano in modo diverso se prefissati da ., infatti:
a = [ 2 2; ...
1 2 ];
a.^3
ans =
8 8
1 8
a^3
ans =
20 28
14 20
a*a*a
ans =
20 28
14 20tutte le funzioni scalari di MATLAB funzionano elemento per elemento, ad esempio
a = [ 2.1 4.2 2.1; ...
1 2.001 -0.1];
sin(a)
ans =
0.8632 -0.8716 0.8632
0.8415 0.9089 -0.0998
cos(a)
ans =
-0.5048 -0.4903 -0.5048
0.5403 -0.4171 0.9950
exp(a)
ans =
8.1662 66.6863 8.1662
2.7183 7.3964 0.9048Attenzione che se invece volessimo fare l’esponenziale di matrice cioè
\exp(\bm{A}) = \bm{I} + \bm{A} + \frac{1}{2}\bm{A}^2 + \cdots
bisogna usare la funzione expm:
>> a = [ 2.1 4.2; ...
1 2.001 ];
>> exp(a) % esponenziale componente per componente
ans =
8.1662 66.6863
2.7183 7.3964
>> expm(a) % esponenziale matriciale
ans =
31.4019 60.8176
14.4804 29.9683Costruzione di successioni
Spesso programmando con MATLAB è conveniente avere degli operatori per costruire delle successioni di interi o reali a passo costante. In MATLAB l’operatore : (due punti) serve allo scopo
>> 1:4 % succesione da 1 a 4
ans =
1 2 3 4
>> 13:34 % succesione da 13 a 34
ans =
Columns 1 through 12
13 14 15 16 17 18 19 20 21 22 23 24
Columns 13 through 22
25 26 27 28 29 30 31 32 33 34
>> 1:4:13 % succesione da 1 a 13 con passo 4
ans =
1 5 9 13
>> 0.24:0.1:3
ans =
Columns 1 through 7
0.2400 0.3400 0.4400 0.5400 0.6400 0.7400 0.8400
Columns 8 through 14
0.9400 1.0400 1.1400 1.2400 1.3400 1.4400 1.5400
Columns 15 through 21
1.6400 1.7400 1.8400 1.9400 2.0400 2.1400 2.2400
Columns 22 through 28
2.3400 2.4400 2.5400 2.6400 2.7400 2.8400 2.9400queste successioni possono essere utilmente utilizzate per estrarre blocchi di matrici, ad esempio
>> a = [ 1 2 3 4 5;
6 7 8 9 10;
11 12 13 14 15 ];
>> a(1:2,3)
ans =
3
8
>> a(2:end,1:end-1)
ans =
6 7 8 9
11 12 13 14la parola chiave end contiene la dimensione massima di riga o di colonna della matrice su cui opera. E’ molto comoda nella operazione sulle matrici. Ad esempio se volessimo ruotare gli elementi di un vettore basta fare le seguenti operazioni
>> v = [ 1 2 3 4 5 ]
v =
1 2 3 4 5
>> w = [ v(end) v(1:end-1) ]
w =
5 1 2 3 4notate che se in un elemento di una matrice si mette un blocco matrice questo viene spiattellato
>> v = [ 1 [1 2 3 4] ]
v =
1 1 2 3 4cioè in MATLAB non si possono fare strutture ricorsive (non con le matrici perlomeno) cioè matrici i cui elementi a loro volta sono matrici. Questa caratteristica in realtà è molto comoda perché diventa molto facile costruire matrici complesse a partire da blocchi semplici, ad esempio
>> I3 = eye(3) % matrice identita' 3x3
I3 =
1 0 0
0 1 0
0 0 1
>> I4 = eye(4) % matrice identita' 4x4
I4 =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
>> M = [ I3 zeros( 3, 4); ones(4,3) I4 ]
M =
1 0 0 0 0 0 0
0 1 0 0 0 0 0
0 0 1 0 0 0 0
1 1 1 1 0 0 0
1 1 1 0 1 0 0
1 1 1 0 0 1 0
1 1 1 0 0 0 1La funzione eye(n) costruisce una matrice identità n\times n, mentre zeros(n,m) e ones(n,m) costruiscono matrici di 0 o di 1 n\times m.
Una scorciatoia per 1:end è usare semplicemente due punti com segue
>> a = [ 1 2 3 4;
5 6 7 8 ];
>> a(:,1)
ans =
1
5
>> a(1:end,1)
ans =
1
5
>> a(2,:)
ans =
5 6 7 8
>> a(2,1:end)
ans =
5 6 7 8Principali funzioni di matrice
Esistono un gran numero di funzioni per manipolare le matrici, ad esempio:
>> M = [ 1 2 3 4;
5 6 7 8;
1 1 1 1;
2 3 1 0 ];Calcola solo gli autovalori di M
>> e = eig(M)
e =
11.5977
-3.0282
-0.5695
-0.0000Calcola gli autovalori e gli autovetori di M
>> [V,e] = eig(M)
V =
0.3200 0.5808 0.8121 0.1667
0.8883 0.3638 -0.4641 0.1667
0.1421 -0.0542 -0.3190 -0.8333
0.2972 -0.7262 0.1527 0.5000
e =
11.5977 0 0 0
0 -3.0282 0 0
0 0 -0.5695 0
0 0 0 -0.0000
>> M * V(:,1) - e(1) * V(:,1) % controlla la prima coppia autovalore autovettore
ans =
1.0e-13 *
0
0.1066
0.0155
0.0178Decompositione LU
>> [L,U,P] = lu(M)
L =
1.0000 0 0 0
0.2000 1.0000 0 0
0.4000 0.7500 1.0000 0
0.2000 -0.2500 -0.0000 1.0000
U =
5.0000 6.0000 7.0000 8.0000
0 0.8000 1.6000 2.4000
0 0 -3.0000 -5.0000
0 0 0 0.0000
P =
0 1 0 0
1 0 0 0
0 0 0 1
0 0 1 0
>> norm(P*L*U - M) % controllo la decomposizione
ans =
0Fattorizzazioen QR
>> [Q,R] = qr(M)
Q =
-0.1796 0.7671 0.5690 -0.2357
-0.8980 -0.3002 0.2188 0.2357
-0.1796 -0.2668 -0.0875 -0.9428
-0.3592 0.5003 -0.7878 -0.0000
R =
-5.5678 -7.0046 -7.3638 -8.0822
0 0.9672 0.4336 0.4002
0 0 2.3635 3.9392
0 0 0 -0.0000
>> norm(Q*R - M) % controllo la fattorizzazione
ans =
1.7681e-15Decomposizione ai valori singolari
>> [U,S,V] = svd(M)
U =
-0.3636 0.4329 0.7904 0.2357
-0.9034 0.0224 -0.3575 -0.2357
-0.1349 -0.1026 -0.2870 0.9428
-0.1829 -0.8953 0.4062 0.0000
S =
14.5996 0 0 0
0 2.9146 0 0
0 0 0.5982 0
0 0 0 0.0000
V =
-0.3686 -0.4627 -0.7888 0.1667
-0.4679 -0.6137 0.6138 0.1667
-0.5296 0.1569 -0.0206 -0.8333
-0.6039 0.6203 0.0240 0.5000
>> norm(U*S*V' - M) % controllo la decompositione
ans =
5.0559e-15Norma 1 di matrice
>> norm(M,1)
ans =
13Norma inifinito di matrice
>> norm(M,'inf')
ans =
26Norma di Frobenius
norm(M,'fro')
ans =
14.8997Alcune utili funzioni di matrice
Ci sono alcune funzioni di matrice che sono molto utili nella costruzione di algoritmi.
Le funzioni max(a) e min(a) restituiscono un vettore riga nel quale ogni elemento è il massimo (o il minimo) dei valori nella colonna. Se la matrice è un vettore riga o colonna max(a) e min(a) restituisce il massimo o minimo del vettore. In modo analogo la funzione sum(a) restituisce il vettore riga con la somma dei valori sulle colonne di a. Se la matrice è un vettore riga o colonna sum(a) restituisce la somma delle componenti del vettore. In modo analogo funziona la funzione prod(a).
>> a = [ 1 2 3 -4;
0 1 0 -2;
2 3 0 -1 ];
>> max(a)
ans =
2 3 3 -1
>> min(a)
ans =
0 1 0 -4
>> sum(a)
ans =
3 6 3 -7
>> max(max(a))
ans =
3
>> min(min(a))
ans =
-4
>> sum(sum(a))
ans =
5
>> max(sum(abs(a))) % norma 1
ans =
7
>> norm(a,1)
ans =
7
>> max(sum(abs(a'))) % norma inifinito
ans =
10
>> norm(a,'inf')
ans =
10Se assegnamo i risultati della funzione max(a) o min(a) ed una coppia di variabili
[res,idx] = max(a)
[res,idx] = min(a)allora res contiene il vettore riga con gli elementi massimi (o minimi) della colonna e idx contiene gli indici di riga dell’elemento massimo (o minimo) ad esempio
>> a = [ 1 2 3 -4;
0 1 0 -2;
2 3 0 -1 ];
>> [res,idx] = max(a)
res =
2 3 3 -1
idx =
3 3 1 3
>> [res,idx] = min(a)
res =
0 1 0 -4
idx =
2 2 2 1La funzione length(a) è equivalente a max(size(a)).
Parte grafica
In MATLAB è molto facile fare il grafico di una funzione con il comando plot, ad esempio il comando
>> x = [0:0.1:10];
plot(x, sin(x))produce la figura
mentre il seguente
>> x = [0:0.1:10];
plot(x,sin(x),'r*')produce la figura
con il comando hold si possono combinare vari plot per ottenere grafici complessi. Ad esempio i comandi
>> plot(x, sin(x))
>> hold on
>> plot(x, sin(x),'r*')
>> xlabel('asse x')
>> ylabel('asse y')
>> title({'un titolo','su piu` righe'})
>> hold offproducono la figura
La parte grafica di MATLAB è molto complessa e estesa, per vedere in dettaglio le varie possibilità basta digitare doc plot nella finestra di comando.
Oltre alle matrici
Nella costruzione di algoritmi complessi a volte non bastano solo vettori e matrici ma occorrono strutture dati più complesse. MATLAB fornisce due possibilità interessanti, i cell array e le strutture. Un cell array può essere visto come un vettore eterogeneo che può contenere qualunque cosa. La sintassi è simile a quelle delle matrici dove vengono sostituite le parentesi tonde e quadre con le parentesi grafe. Ad esempio
>> a = { [1 2 3], {'una stringa'}, 2, [ 1 2; 3 4] } % esempio di cell array
a =
[1x3 double] 'una stringa' [2] [2x2 double]
>> disp('accedo ad a{1}'); a{1}
accedo ad a{1}
ans =
1 2 3
>> disp('accedo ad a{2}'); a{2}
accedo ad a{2}
ans =
una stringa
>> disp('accedo ad a{3}'); a{3}
accedo ad a{3}
ans =
2
>> disp('accedo ad a{4}'); a{4}
accedo ad a{4}
ans =
1 2
3 4Il seguente esempio fa vedere come si possono costruire strutture molto complesse:
% esempio di cell array
>> a = { [1 2 3], { 1 2 3 { 2 [ 1 2; 3 4] } } }
a =
[1x3 double] {1x4 cell}
>> disp('accedo ad a{1}'); a{1}
accedo ad a{1}
ans =
1 2 3
>> disp('accedo ad a{2}'); a{2}
accedo ad a{2}
ans =
[1] [2] [3] {1x2 cell}
>> disp('accedo ad a{2}{1}'); a{2}{1}
accedo ad a{2}{1}
ans =
1
>> disp('accedo ad a{2}{2}'); a{2}{2}
accedo ad a{2}{2}
ans =
2
>> disp('accedo ad a{2}{3}'}); a{2}{3}
accedo ad a{2}{3}
ans =
3
>> disp('accedo ad a{2}{4}'); a{2}{4}
accedo ad a{2}{4}
ans =
[2] [2x2 double]
>> disp('accedo ad a{2}{4}{1}'); a{2}{4}{1}
accedo ad a{2}{4}{1}
ans =
2
>> disp('accedo ad a{2}{4}{2}'); a{2}{4}{2}
accedo ad a{2}{4}{2}
ans =
1 2
3 4nell’esempio precedente si nota l’uso della funzione disp che stampa una stringa a terminale.
Per accedere agli elementi di un cell array si usano gli interi. Un struttura è simile al cell array nel senso che permette di costruite strutture dati complesse ma l’accesso agli elementi è tramite l’operatore punto . e i nomi dei campi corrispondenti. Ad esempio
% esempio di sruttura
>> a . pippo = [ 1 2 3 ]; % il campo pippo e' un vettore riga
>> a . pluto = [ 1; 3 ]; % il campo pluto e' un vettore colonna
>> a . paperino = { 'a' 'b' 3 }; % il campo paperino e' un cell array
>> a . subs.a = 1; % il campo subs e' una struttura a sua volta
>> a . subs.b = 2; % il campo subs e' una struttura a sua volta
>> a . subs.c = 3; % il campo subs e' una struttura a sua volta
>> disp(a)
pippo: [1 2 3]
pluto: [2x1 double]
paperino: {'a' 'b' [3]}
subs: [1x1 struct]
>> disp(a.subs)
a: 1
b: 2
c: 3Strutture e matrici si possono combinare. Ad esempio si possono fare vettori di strutture e strutture di vettori.
L’istruzione find
MATLAB supporta le istruzioni condizionali if–else–end cicli while–end e for–end. Con queste istruzioni si può programmare in maniera procedurale, ma non è la maniera giusta di programmare in MATLAB. Se si riesce è meglio cercare di costruire gli algoritmi in maniera vettoriale cercando di fare blocchi di istruzioni che operano su matrice e vettori e non sulle loro componenti. A questo proposito la funzione find permette di vettorizzare cicli condizionali. L’istruzione find si può usare essenzialmente in 2 modi:
Modo 1
Si usa find per farsi restituire gli indici di riga e colonna che soddisfano una condizione logica, ad esempio
>> A = [ 1 2 3 -4;
0 1 0 -2;
2 3 0 -1 ];
>> [I,J] = find( A > 0 )
I =
1
3
1
2
3
1
J =
1
1
2
2
2
3
>> [I,J] = find( A > 0 \& A < 2 )
I =
1
2
J =
1
2Modo 2
Si usa find per farsi restituire un indice corrispondente che soddisfano una condizione logica, ad esempio
>> A = [ 1 2 3 -4;
0 1 0 -2;
2 3 0 -1 ];
>> [IDX] = find( A > 0 )
IDX =
1
3
4
5
6
7
>> [IDX] = find( A > 0 \& A < 2 )
IDX =
1
5Possiamo usare la funzione find per fare cose che normalmente faremmo con dei cicli. Ad esempio se vogliamo cambiare tutte le componenti negative di una matrice con i loro quadrati:
>> A = [ 1 2 3 -4;
0 1 0 -2;
2 3 0 -1 ]
>> IDX = find( A < 0 )
IDX =
10
11
12
>> A(IDX) = A(IDX) .^ 2
A =
1 2 3 16
0 1 0 4
2 3 0 1Instruzioni condizionali
L’istruzione if–elseif–else–end:
>> a = 1;
>> if a > 0
disp('a>0');
else
disp('a<=0');
end
a>0
>> b = 0;
>> if b > 0
disp('b>0');
elseif b < 0
disp('b<0');
else
disp('b==0');
end
b==0Il ciclo for–end:
>> x = [];
>> for i=1:10
x = [ x i ];
end
>> disp(x);
1 2 3 4 5 6 7 8 9 10
>> x = [];
>> for i=10:-1:1
x = [ x i ];
end
>> disp(x);
10 9 8 7 6 5 4 3 2 1
>> A = zeros(5,5);
>> for i=1:5
for j=1:5
A(i,j) = (i-j)^2 + i*j;
end
end
>> disp(A);
1 3 7 13 21
3 4 7 12 19
7 7 9 13 19
13 12 13 16 21
21 19 19 21 25Il ciclo while–end:
>> a = 129;
>> b = 0;
>> while a > 1
a = floor(a/2);
b = b + 1;
end
>> disp( b );
7Ogni ciclo for o while può essere interrotto da una istruzione break in ogni punto del ciclo.
Confronti
Nelle istruzioni condizionali abbiamo già visto alcune istruzioni condizionali MATLAB supporta le seguenti
eq,==uguaglianzane,~=diseguaglianzalt,<"minoregt,>"maggiorele,<=minore o ugualege,>=maggiore o uguale
e gli operatori logici
and,&andor,|ornot,~notxoror esclusivoanyvera se qualche elemento del vettore è non zeroallvera se ogni elemento del vettore è non zero
Gli script
In MATLAB si possono definire funzioni e scriverle su file esterni che poi possono essere richiamati da riga di comando. Il file deve chiamarsi comando.m dove comando è il nome del comando che si vuole usare. All’interno di comando.m si troverà la definizione della funzione comando tipo
function res = comando( x, y, x )
...consideriamo ad esempio il file distanza.m
function res = distanza( x, y )
res = sqrt( sum((x-y).^2) );la funzione distanza(x,y) restituisce la distanza tra due vettori vettori, ad esempio:
>> d = distanza( [1 2 3], [ 0 1 -2]);
>> disp(d)
5.1962Alcune considerazioni
A differenza di altri linguaggi di programmazione il blocco
functionnon necessariamente è terminato daend.Il nome della funzione è necessariamente lo stesso del file
.m.È possibile definire altre funzioni all’interno dello stesso script. Queste funzioni sono pero visibili solo all’interno dello script.
Gli argomenti sono passati per valore (non è esattamente così ma rende l’idea) quindi se modificati all’interno della funzione le modifiche non sono viste dal programma chiamante.
Restituzione di più valori
Se servono più valori da restituire questi si mettono in una lista. Ad esempio dato lo script trenorme.m
% calcola varie norme di un vettore
%
% _______
% n / n
% ===== / =====
% \ / \
% max |x | , > |x | , / > x
% | i| / | i| / /
% ===== / =====
% i = 1 \/ i = 1
%
function [norma2,norma1,normainf] = trenorme( x )
norma2 = sqrt(sum(x^2));
norma1 = sum(abs(x));
normainf = max(abs(x));calcola le norme 1, 2 e \infty di un vettore. Il seguente comando produce:
>> [a,b,c] = trenorme([1 2 3])
a =
3.7417
b =
6
c =
3Numero variabile di argomenti
Una altra caratteristica molto utile di MATLAB è la possibilità di usare un numero variabile di argomenti in ingresso e uscita. Per avere più informazioni su questa possibilità usare
>> doc varargin
>> doc varargoutEsempio complesso
%
% Risolve il sistema lineare A*x = b
%
function x = linearsistemsolve( A, b )
[nr,nc] = size(A);
% controllo dimensioni
if nr ~= nc
error( sprintf( 'la matrice A (%d x %d) non è quadrata ', nr,nc) );
end
neq = length(b);
if nr ~= neq
error( sprintf( 'la matrice A (%d x %d) non è compatibile con la lunghezza del vettore b, length(b) = %d ', nr, nc, neq ) );
end
M = [ A b ]; % costruisco la matrice aumentata
for i=1:neq-1
% cerco l'elemento di massimo modulo colonna i-esima
[V,k] = max(abs(M(i:end,i)));
% rimetto l'indice k alla giusta posizione
k = k + i-1;
% scambio la riga i con la riga k
BF = M(i,:);
M(i,:) = M(k,:);
M(k,:) = BF;
% controllo singolarita'
if V == 0
error( 'Matrice singolare' );
end
% costruisco la matrice di Frobenius
L = eye(neq);
L(i+1:end,i) = -M(i+1:end,i)./M(i,i);
% azzero gli elementi della colonna i
M = L * M;
end
% alloco il vettore soluzione
x = zeros(neq,1);
% a questo punto la matrice e' in forma triangolare
% applichiamo le operazioni di ritorno
x(neq) = M(neq,neq) / M(neq,neq+1);
for i=neq-1:-1:1
x(i) = ( M(i,neq+1) - M(i,i+1:end-1) * x(i+1:end) ) / M(i,i);
endEsempio di uso
>> A = [ 1 2 3; 1 2 4; 2 3 4 ];
>> b = [ 1 2 3 ]
>> x = linearsistemsolve(A,b)
x =
4
-3
1Bibliografia
- Davis, T. A. (2005). MATLAB Primer (7th ed.). Chapman & Hall/CRC.
- Otto, S. R., & Denier, J. P. (2005). An Introduction to Programming and Numerical Methods in MATLAB. Springer.