Metode numerice în MATLAB

Table of Contents
Rezolvarea sistemelor algebrice liniare Sisteme pătratice Sisteme supradeterminate Sisteme subdeterminate Factorizarea QR Vectori și valori proprii eig și condeig Probleme proprii generalizate SVD Polinoame şi potrivirea datelor (data fitting) Evaluarea și determinarea rădăcinilor Convoluția Potrivirea datelor Aplicație: recensământul SUA

Rezolvarea sistemelor algebrice liniare

Fie m numărul de ecuaţii şi n numărul de necunoscute. Instrumentul fundamental de rezolvare a sistemelor de ecuaţii liniare este operatorul \.
El tratează trei tipuri de sisteme de ecuaţii liniare, pătratice (), supradeterminate () şi subdeterminate (). Mai general, operatorul \ poate fi utilizat pentru a rezolva , unde B este o matrice cu p coloane; în acest caz MATLAB rezolvă sistemele AX(:; j) = B(:; j) pentru j = 1 : p. Sistemele de forma se pot rezolva cu X = B/A.
Operatorul \ se bazează pe algoritmi diferiţi în funcţie de matricea coeficienţilor. Diversele cazuri, care sunt diagnosticate automat prin examinarea matricei sistemului includ:

Sisteme pătratice

Dacă A este o matrice pătratică nesingulară de ordinul n, atunci A\b este soluţia sistemului , calculată prin factorizare LU cu pivotare parţială. În timpul rezolvării, MATLAB calculează rcond(A) şi tipăreşte un mesaj de avertisment dacă rezultatul este mai mic decât eps:
x = hilb(15)\ones(15,1)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.460912e-19.
x = 15×1
10.461 -1498.8 51797 -7.469e+05 5.4714e+06 -2.1654e+07 4.2474e+07 -1.158e+07 -1.1811e+08 2.0545e+08

Sisteme supradeterminate

Dacă , în general sistemul nu are nici o soluţie. Expresia MATLAB A\b calculează soluţia sistemului în sensul celor mai mici pătrate, adică minimizează norma euclidiană a reziduului (adică norm(A*x-b)) peste toţi vectorii x. Dacă A are rang maxim m, atunci soluţia este unică. Dacă A are rangul k mai mic decât m (în acest caz spunem că A este deficientă de rang), A\b calculează o soluţie de bază cu cel mult m elemente nenule (k este determinat şi x este calculat utilizând factorizarea QR cu pivotare pe coloană). În ultimul caz MATLAB dă un mesaj de avertisment.
Soluţia se mai poate calcula şi cu x_min=pinv(A)*b, unde pinv(A) este pseudo-inversa lui A. Dacă A este deficientă de rang, x_min este soluţia unică de normă euclidiană minimală.
Pseudo-inversa Moore-Penrose a lui A, notată cu generalizează noţiunea de inversă pentru matrice dreptunghiulare şi deficiente de rang. Pseudo-inversa a lui A este matricea unică care satisface condiţiile
Vom ilustra cu următoarele exemple:
Y=pinv(ones(3))
Y = 3×3
0.11111 0.11111 0.11111 0.11111 0.11111 0.11111 0.11111 0.11111 0.11111
A=[0 0 0 0; 0 1 0 0; 0 0 2 0]; pinv(A)
ans = 4×3
0 0 0 0 1 0 0 0 0.5 0 0 0

Sisteme subdeterminate

În cazul unui sistem subdeterminat putem fie să nu avem nici o soluţie fie să avem o infinitate. În ultimul caz, A\b produce o soluţie de bază cu cel mult k elemente nenule, unde k este rangul lui A. În general, această soluţie nu are norma euclidiană minimă; soluţia cu norma minimă se poate calcula cu pinv(A)*b. Dacă sistemul nu are nici o soluţie (este incompatibil), atunci A\b este o soluţie în sensul celor mai mici pătrate. Exemplul următor ilustrează diferenţa dintre soluţia obţinută cu \ şi pinv:
A = [1 1 1; 1 1 -1]; b=[3; 1];
x=A\b; y = pinv(A)*b;
[x y]
ans = 3×2
2 1 0 1 1 1
[norm(x) norm(y)]
ans = 1×2
2.2361 1.7321
MATLAB foloseşte factorizarea QR cu pivotare pe coloană. Fie exemplul
R=fix(10*rand(2,4)), b=fix(10*rand(2,1))
R = 2×4
7 5 0 7 3 0 5 9
b = 2×1
1 5
Sistemul are 2 ecuaţii şi 4 necunoscute. Deoarece matricea coeficienţilor conţine întregi mici, este recomandabil să afişăm soluţia în format raţional. Soluţia particulară se obţine cu:
format rat
p=R\b
p =
0 -26/45 0 5/9
O componentă nenulă este p(2), deoarece R(:,2) este coloana cu cea mai mare normă. Cealaltă este p(4), deoarece R(:,4) rămâne dominantă după eliminarea lui R(:,2).
Soluţia completă a unui sistem supradeterminat poate fi caracterizată prin adăugarea unui vector arbitrar din spaţiul nul al matricei sistemului, care poate fi găsit cu funcţia null cu o opţiune care cere o bază raţională.
Z=null(R,'r')
Z =
-5/3 -3 7/3 14/5 1 0 0 1
Se poate verifica că R*Z este zero şi orice vector de forma x=p+Z*q, unde q este un vector arbitrar, satisface R*x=b.
R*Z
ans =
0 0 0 0
x=p+Z*ones(2,1)
x =
-14/3 41/9 1 14/9
format short
r=R*x-b
r = 2×1
10-14 ×
-0.3553 0
Factorizarea LU şi Cholesky
Funcţia lucalculează o factorizare LUP cu pivotare parţială. Apelul [L,U,P]=lu(A) returnează factorii triunghiulari şi matricea de permutare. Forma [L,U]=lu(A) returnează , deci L este o matrice triunghiulară cu liniile permutate.
format short g
A = gallery('fiedler',3), [L,U]=lu(A)
A = 3×3
0 1 2 1 0 1 2 1 0
L = 3×3
0 1 0 0.5 -0.5 1 1 0 0
U = 3×3
2 1 0 0 1 2 0 0 2
Factorizarea LU este definită şi pentru matrice dreptunghiulare:
[L,U,P]=lu([1,2,3;3,0,4])
L = 2×2
1 0 0.33333 1
U = 2×3
3 0 4 0 2 1.6667
P = 2×2
0 1 1 0
Rezolvarea sistemului Ax=b cu x=A/b şi cu A pătratică este echivalentă cu secvenţa MATLAB:
[L,U] = lu(A); x = U\(L\b);
Determinantul şi inversa se calculează de asemenea prin factorizare LU:
det(A)=det(L)*det(U)=+-prod(diag(U))
inv(A)=inv(U)*inv(L)
Comanda chol(A), unde A este hermitiană şi pozitiv definită calculează R triunghiulară superior astfel încât . Factorizarea Cholesky permite înlocuirea sistemului A*x=b cu R’*R*x=b. Deoarece operatorul \ recunoaşte sisteme triunghiulare, sistemul se poate rezolva rapid cu x=R\(R’\R\b). Dăm un exemplu de factorizare Cholesky:
A=pascal(4)
A = 4×4
1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20
R=chol(A)
R = 4×4
1 1 1 1 0 1 2 3 0 0 1 3 0 0 0 1

Factorizarea QR

Factorizarea QR a lui A este o descompunere de forma , unde Q este unitară, adică (ortogonală în cazul real ) și R este triunghiulară superior, cu . Factorizarea QR se aplică la rezolvarea sistemelor, chiar și a celor care nu sunt pătratice (este mai stabilă decât LUP) și la rezolvarea problemelor proprii.
A=[1,2,3; 3,2,1; 1, -4, 1];
b=[6;6;-2];
[Q,R]=qr(A)
Q = 3×3
-0.30151 0.34463 0.889 -0.90453 0.19146 -0.381 -0.30151 -0.91901 0.254
R = 3×3
-3.3166 -1.206 -2.1106 0 4.7482 0.30634 0 0 2.54
x=R\(Q'*b)
x = 3×1
1 1 1

Vectori și valori proprii

eig și condeig

MATLAB utilizează rutine LAPACK pentru a calcula valori şi vectori proprii. Valorile proprii ale unei matrice se calculează cu funcţia eig: e=eig(A) pune valorile proprii ale lui A în vectorul e. Forma [V,D]=eig(A), unde A este matrice pătratică de ordinul n, returnează în coloanele lui V n vectori proprii ai lui A şi în matricea diagonală D valorile proprii ale lui A. Are loc relaţia A*V=V*D. Nu orice matrice are n vectori proprii liniari independenţi, deci matricea V returnată de eig poate fi singulară (sau datorită erorilor de rotunjire nesingulară, dar foarte prost condiţionată).
Matricea din exemplul următor are o valoare proprie dublă 1, dar numai un vector propriu linear independent:
[V,D]=eig([2, -1; 1,0])
V = 2×2
0.70711 0.70711 0.70711 0.70711
D = 2×2
1 0 0 1
Vectorii proprii sunt scalaţi astfel ca norma lor euclidiană să fie egală cu unu (lucru posibil, căci dacă x este un vector propriu, atunci orice multiplu al său este de asemenea vector propriu).
Dacă A este hermitiană, atunci MATLAB returnează valorile proprii sortate crescător şi matricea vectorilor proprii unitară (în limita preciziei de lucru):
[V,D]=eig([2,-1;-1,1])
V = 2×2
-0.52573 -0.85065 -0.85065 0.52573
D = 2×2
0.38197 0 0 2.618
norm(V'*V-eye(2))
ans =
2.3417e-16
În exemplul următor vom calcula valorile proprii ale matricei lui Frank (nehermitiană):
F = gallery('frank',5)
F = 5×5
5 4 3 2 1 4 4 3 2 1 0 3 3 2 1 0 0 2 2 1 0 0 0 1 1
e = eig(F)'
e = 1×5
10.063 3.5566 1 0.099375 0.28117
Dacă λ este valoare proprie a matricei F, atunci este de asemenea valoare proprie a lui F:
1./e
ans = 1×5
0.099375 0.28117 1 10.063 3.5566
Motivul este acela că polinomul caracteristic este anti-palindromic, adică termenii egal depărtaţi de extrem sunt numere opuse:
poly(F)
ans = 1×6
1 -15 55 -55 15 -1
Funcţia condeig calculează numărul de condiţionare pentru valori proprii. Acesta se defineşte prin
Forma c=condeig(A) returnează un vector al numerelor de condiţionare ale valorilor proprii ale lui A. Forma [V,D,s] = condeig(A) este echivalentă cu
[V,D] = eig(A),
s = condeig(A).
Un număr de condiţionare mare indică o valoare proprie sensibilă la perturbaţii ale matricei. Exemplul următor afişează în prima linie valorile proprii ale matricei lui Frank de ordinul 6 şi în a doua linie numerele lor de condiţionare:
A = gallery('frank',6);
[V,D,s] = condeig(A);
[diag(D)'; s']
ans = 2×6
12.974 5.3832 1.8355 0.5448 0.07708 0.18576 1.3059 1.3561 2.0412 15.326 43.521 56.695

Probleme proprii generalizate

MATLAB poate rezolva şi probleme de valori proprii generalizate, adică probleme de forma: fiind date două matrice pătratice de ordinul n, A şi B, să se găsească scalarii λ şi vectorii x astfel încât . Valorile proprii generalizate se pot calcula cu e = eig(A,B), iar [V,D] = eig(A,B) returnează o matrice diagonală D a valorilor proprii şi o matrice pătratică de ordinul n a vectorilor proprii V astfel încât A*V=B*V*D. Teoria corespunzătoare este mai complicată decât cea a valorilor proprii standard: putem să nu avem nici o valoare proprie, putem avea un număr finit de valori proprii sau o infinitate, sau valori proprii infinit de mari. Dacă B este singulară, se pot obţine valori proprii NaN. Dăm un exemplu de rezolvare a unei probleme proprii generalizate:
A = gallery('triw',3), B = magic(3)
A = 3×3
1 -1 -1 0 1 -1 0 0 1
B = 3×3
8 1 6 3 5 7 4 9 2
[V,D]=eig(A,B); V, eigvals = diag(D)'
V = 3×3
-1 -1 0.35263 0.48441 -0.4574 0.3867 0.21993 -0.25164 -1
eigvals = 1×3
0.27507 0.029193 -0.34593

SVD

Se numeşte descompunere cu valori singulare (singular value decomposition – SVD) descompunerea
unde Σ este o matrice diagonală reală, iar U şi V sunt matrice unitare (ortogonale în cazul real).
Există două variante de SVD: una completă, care se aplică unei matrice dreptunghiulare şi care returnează matricele U de dimensiune , Σ de dimensiune şi V de dimensiune şi una economică sau redusă în care U are dimensiunea , Σ are dimensiunea şi V are dimensiunea .
SVD este un instrument util de analiză a aplicaţiilor dintr-un spaţiu vectorial cu valori în alt spaţiu, posibil de dimensiune diferită. Dacă A este pătratică, simetrică şi pozitiv definită SVD şi descompunerea cu valori proprii coincid. Spre deosebire de descompunerea cu valori proprii, SVD există întotdeauna.
Fie matricea:
A = [9, 4; 6, 8; 2 7];
Descompunerea sa SVD completă este
[U,S,V]=svd(A)
U = 3×3
-0.6105 0.71744 0.33552 -0.66459 -0.23361 -0.70975 -0.43082 -0.65629 0.61942
S = 3×2
14.936 0 0 5.1883 0 0
V = 2×2
-0.69254 0.72138 -0.72138 -0.69254
iar cea redusă
[U,S,V]=svd(A,0)
U = 3×2
-0.6105 0.71744 -0.66459 -0.23361 -0.43082 -0.65629
S = 2×2
14.936 0 0 5.1883
V = 2×2
-0.69254 0.72138 -0.72138 -0.69254
În ambele cazuri se poate verifica că U*S*V’ este egală cu A, în limita erorilor de rotunjire.
norm(U*S*V'-A)
ans =
2.5942e-15

Polinoame şi potrivirea datelor (data fitting)

MATLAB reprezintă un polinom
printr-un vector linie p=[p(1), p(2), ..., p(n+1)] al coeficienţilor, ordonaţi descrescător după puterile variabilei.
Să considerăm trei probleme legate de polinoame:

Evaluarea și determinarea rădăcinilor

Evaluarea se face cu ajutorul schemei lui Horner, implementată în MATLAB prin funcţia polyval. În comanda y=polyval(p,x) x poate fi o matrice, în acest caz evaluarea făcându-se element cu element (deci, în sens tablou). Evaluarea în sens matricial, adică obţinerea matricei
,
unde X este o matrice pătratică, se poate face cu comanda Y = polyvalm(p,X).
Rădăcinile (reale şi complexe) ale polinomului p se pot obţine cu z = roots(p). Funcţia poly realizează operaţia inversă, adică construieşte polinomul cunoscând rădăcinile. Ea acceptă ca argument şi o matrice pătratică A, caz în care p=poly(A) calculează polinomul caracteristic al lui A, adică .
Funcţia polyder calculează coeficienţii derivatei unui polinom, fără a-l evalua.
Ca exemplu, să considerăm polinomul . Rădăcinile lui le obţinem cu
p = [1 -1 -1]; z = roots(p)
z = 2×1
-0.61803 1.618
Verificăm, în limitele erorilor de rotunjire, că acestea sunt rădăcinile:
polyval(p,z)
ans = 2×1
-1.1102e-16 2.2204e-16
Observăm că p este polinomul caracteristic al unei anumite matrice
A = [0 1; 1 1]; cp = poly(A)
cp = 1×3
1 -1 -1
Teorema Cayley-Hamilton ne spune că orice matrice satisface polinomul său caracteristic. Aceasta se verifică în limita erorilor de rotunjire:
polyvalm(cp, A)
ans = 2×2
1.1102e-16 0 0 1.1102e-16

Convoluția

Înmulţirea şi împărţirea polinoamelor se realizează cu conv şi deconv. Sintaxa lui deconv este [q,r]=deconv(g,h), unde g este deîmpărţitul, h împărţitorul, q câtul şi r restul. În exemplul următor vom împărţi la , obţinând câtul şi restul 0. Polinomul iniţial se va obţine apoi cu conv.
g = [1 -2 -1 2]; h=[1 -2];
[q,r] = deconv(g,h)
q = 1×3
1 0 -1
r = 1×4
0 0 0 0
conv(h,q)+r
ans = 1×4
1 -2 -1 2

Potrivirea datelor

Să presupunem că avem m observaţii măsurate în valorile specificate :
Modelul nostru este o combinaţie de n funcţii de bază
Matricea de proiectare (design matrix) va fi matricea cu elementele
ale cărei elemente pot depinde de α. În notaţie matricială, modelul se poate exprima ca:
Reziduurile sunt diferenţele dintre valorile observate şi cele date de model
sau, în notaţie matricială,
Ne propunem să minimizăm o anumită normă a reziduurilor. Cele mai frecvente alegeri sunt
și
O explicaţie intuitivă, fizică, a celei de-a doua alegeri ar fi aceea că anumite observaţii sunt mai importante decât altele şi le vom asocia ponderi, . De exemplu, dacă la observaţia i eroarea este aproximativ , atunci putem alege . Deci, avem de a face cu o problemă discretă de aproximare în sensul celor mai mici pătrate. Problema este liniară dacă nu depinde de α şi neliniară în caz contrar.
Orice algoritm de rezolvare a unei probleme de aproximare în sensul celei mai mici pătrate fără ponderi poate fi utilizat la rezolvarea unei probleme cu ponderi prin scalarea observaţiilor şi a matricei de proiectare. În MATLAB aceasta se poate realiza prin
A=diag(w)*A
y=diag(w)*y
Dacă problema este liniară şi avem mai multe observaţii decât funcţii de bază, suntem conduşi la rezolvarea sistemului supradeterminat , pe care îl rezolvăm prin metoda celor mai mici pătrate.
Abordarea teoretică se bazează pe rezolvarea ecuaţiilor normale
Dacă funcţiile de bază sunt liniar independente şi deci nesingulară, soluţia este
sau
unde este pseudo-inversa lui A. Ea se poate calcula cu funcţia MATLAB pinv.
Fie sistemul arbitrar. Dacă A este o matrice cu şi are rangul n, atunci fiecare din următoarele trei instrucţiuni
x=A\b
x=pinv(A)*b
x=inv(A'*A)*A'*b
calculează aceeaşi soluţie în sensul celor mai mici pătrate, deşi operatorul \ o face cel mai repede.
Totuşi, dacă A nu are rangul complet, soluţia în sensul celor mai mici pătrate nu este unică. Există mai mulţi vectori care minimizează norma . Soluţia calculată cu x=A\b este o soluţie de bază; ea are cel mult r componente nenule, unde r este rangul lui A. Soluţia calculată cu x=pinv(A)*b este soluţia cu normă minimă (ea minimizează norm(x)). Încercarea de a calcula o soluţie cu x=inv(A’*A)*A’*b eşuează dacă A’*A este singulară. Iată un exemplu care ilus­trează diversele situaţii.
Matricea
A=[1,2,3; 4,5,6; 7,8,9; 10,11,12];
este deficientă de rang. Dacă
b=A(:,2);
atunci o soluţie evidentă a lui A*x=b este x=[0,1,0]'. Nici una dintre abordările de mai sus nu calculează pe x.
Operatorul \ ne dă
x=A\b
Warning: Rank deficient, rank = 2, tol = 1.459426e-14.
x = 3×1
0.5 0 0.5
Această soluţie are două componente nenule. Varianta cu pseudoinversă ne dă
x=pinv(A)*b
x = 3×1
0.33333 0.33333 0.33333
Se observă ca norm(y)=0.5774<norm(x)=0.7071.
A treia variantă eşuează complet:
z=inv(A'*A)*A'*b
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.712392e-18.
z = 3×1
3.5625 -0.625 4.5
Abordarea bazată pe ecuaţii normale are mai multe dezavantaje.
Ecuaţiile normale sunt întotdeauna mai prost condiţionate decât sistemul supradeterminat iniţial.
Numărul de condiţionare se ridică de fapt la pătrat:
Pentru o matrice dreptunghiulară X, numărul de condiţionare ar putea fi definit prin
În reprezentarea în virgulă flotantă, chiar dacă coloanele lui A sunt liniar independente, ar putea fi aproape singulară.
MATLAB evită ecuaţiile normale. Operatorul \ foloseşte intern factoriza­rea QR. Soluţia se poate exprima prin c=R\(Q'*y).
Dacă baza în care se face aproximarea este , se poate folosi funcţia polyfit. Comanda p=polyfit(x,y,n) calculează coeficienţii polinomului de aproximare discretă de grad n în sensul celor mai mici pătrate pentru datele x şi y.
Dacă , se returnează coeficienţii polinomului de interpolare.
O cantitate y este măsurată în diferite momente de timp, t, pentru a produce următoarele observaţii:
tab1mcmmp.png
Aceste date pot fi introduse MATLAB prin
t=[0,0.3,0.8,1.1,1.6,2.3]';
y=[0.82,0.72,0.63,0.60,0.55,0.50]';
Vom încerca să modelăm datele cu ajutorul unei funcţii de forma Coeficienţii necunoscuţi se vor calcula prin metoda celor mai mici pătrate. Avem 6 ecuaţii şi două necunoscute, reprezentate printr-o matrice
E=[ones(size(t)),exp(-t)]
E = 6×2
1 1 1 0.74082 1 0.44933 1 0.33287 1 0.2019 1 0.10026
c=E\y
c = 2×1
0.47595 0.34132
Urmează reprezentarea grafică pe puncte echidistante, completată cu datele originale:
T=[0:0.1:2.5]';
Y=[ones(size(T)),exp(-T)]*c;
plot(T,Y,'-',t,y,'o')
xlabel('t'); ylabel('y');
Se poate vedea că , dar diferenţa este minimă în sensul celor mai mici pătrate. Dacă matricea este deficientă de rang (adică nu are coloane liniar independente), atunci soluţia în sensul celor mai mici pătrate a sistemului nu este unică. În acest caz operatorul \ dă un mesaj de avertizare şi produce o soluţie de bază cu cel mai mic număr posibil de elemente nenule.

Aplicație: recensământul SUA

Al doilea exemplu are ca date de intrare rezultatele de recensământelor obţinute de U. S. Census pentru anii 1900–2010, din zece în zece ani, exprimate în milioane de oameni:
%CENSUS - example with census data
% polynomial fit
%population
y = [ 75.995 91.972 105.711 123.203 131.669 150.697 ...
179.323 203.212 226.505 249.633 281.422 308.786]';
t = (1900:10:2010)'; % census years
Dorim să modelam datele și să facem predicții pentru anii 1975 și 2015
x = (1890:1:2019)'; % evaluation years
w = [1975,2015]; % prediction years
Vom modela datele printr-un polinom de gradul 3
Fitting direct:
c=polyfit(t,y,3)
Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.
c = 1×4
-6.5861e-06 0.047733 -109 80029
Echivalent
A=[t.^3,t.^2,t, ones(size(t))];
c=A\y
Warning: Rank deficient, rank = 3, tol = 6.912993e-05.
c = 4×1
4.1278e-06 -0.015109 13.846 0
Rezultatele sunt inacceptabile din punct de vedere numeric.
Vom normaliza datele cu
mt=mean(t); st=std(t);
s=(t-mt)/st;
xs=(x-mt)/st;
coeficienții
cs=polyfit(s,y,3)
cs = 1×4
-0.30871 11.837 76.561 166.49
sau echivalent
A=[s.^3,s.^2,s, ones(size(s))];
c=A\y
c = 4×1
-0.30871 11.837 76.561 166.49
predicțiile:
zs=polyval(cs,xs);
est=polyval(cs,(w-mt)/st);
reprezentare grafică
plot(t,y,'o',x,zs,'-',w,est,'*')
for i=1:length(w)
text(w(i),est(i)-20,num2str(est(i)))
end
title('U.S. Population', 'FontSize', 14)
xlabel('year', 'FontSize', 12)
ylabel('Millions', 'FontSize', 12)