Sistemi linearnih jednačina

U nastavku su razmatrane implementacije svih iterativnih metoda za rešavanje sistema linearnih jednačina koje su obrađene na časovima vežbi, kao i Gausove metode koja dolazi iz grupe direktnih metoda. Sisteme ćemo posmatrati u matričnom obliku:

Iterativne metode

Iterativne metode određuju niz aproksimacija rešenja koji konvergira tačnom rešenju. Dakle, imamo niz sve kvalitetnijih rešenja, a kriterijumi zaustavljanja su najčešće dostignut maksimalan broj iteracija ili zadovoljena unapred odabrana tačnost rešenja. Iterativni niz dobijamo tako što najpre transformišemo dati sistem u ekvivalentan oblik

odakle dobijamo niz

Taj niz ne mora uvek da konvergira, a potreban i dovoljan uslov konvergencije je da važi, gde je spektralni radijus matrice i predstavlja najveću po modulu sopstvenu vrednost matrice .

Apriorna ocena greške:

Aposteriorna ocena greške:

Radićemo u sledećoj vektorskoj i pridruženoj matričnoj normi:

U zavisnosti od odabira matrice i vektora razlikujemo više iterativnih metoda:

Metoda proste iteracije:

Ovako se može zvati bilo koja iterativna metoda koja konvergira rešenju, ali mi ćemo tako zvati metodu kod koje je

gde je dijagonalna matrica koja odgovara matrici .

% function [x,err] = SimpleIteration(A,b,k)

% [n, m] = size(A);

% if n~=m

% error('Matrica nije kvadratna!');

% end

% I = eye(n);

% D = diag(diag(A));

% B = (D+I)\(D+I-A); %MATLAB preporucuje ovo umesto inv(D+I)*(D+I-A)

% c = (D+I)\b;

% if max(abs(eig(B)))>1

% error('Metoda ne konvergira! Spektralni radius > 1');

% end

% x = c;

% fprintf('Pocetna aproksimacija:')

% disp(x');

% fprintf('Rezultati kroz iteracije:')

% for i = 1:k

% xOld = x;

% x = B*xOld+c;

% disp(x');

% end

% Bnorm = max(sum(abs(B'))); %uniformna norma bez ugradjene funkicje norm()

% xNorm = max(abs(x-xOld));

% err = Bnorm/(1-Bnorm)*xNorm;

% end

A = [3 0.5 -2; -0.2 2 -0.9; -2 -2.5 5];

b = [-1.5 -0.9 -0.5]';

k = 15;

[x,err] = SimpleIteration(A,b,k)

Pocetna aproksimacija:
   -0.3750   -0.3000   -0.0833
Rezultati kroz iteracije:
   -0.4729   -0.4500   -0.3472
   -0.6106   -0.5857   -0.4863
   -0.6976   -0.6818   -0.6120
   -0.7702   -0.7574   -0.7020
   -0.8238   -0.8144   -0.7726
   -0.8655   -0.8582   -0.8260
   -0.8971   -0.8916   -0.8671
   -0.9214   -0.9171   -0.8984
   -0.9399   -0.9366   -0.9223
   -0.9541   -0.9516   -0.9406
   -0.9649   -0.9630   -0.9546
   -0.9731   -0.9717   -0.9653
   -0.9795   -0.9784   -0.9735
   -0.9843   -0.9835   -0.9797
   -0.9880   -0.9874   -0.9845
x =
-0.9880 -0.9874 -0.9845
err = 0.0525

Jakobijeva metoda:

Ova metoda se dobija za sledeći izbor:

Dovoljan uslov da Jakobijeva metoda konvergira je da matrica sistema bude dijagonalno dominantna, odnosno da apsolutna vrednost svakog elemetna na glavnoj dijagonali bude veća od zbira apsolutnih vrednosti ostalih elemenata iz zajedničke vrste.

% function x = Jacobi(A,b,tol)

% n = length(A) %pretpostavimo da je A pravilno uneta

% %Provera da li je matrica dijagonalno dominantna:

% for i = 1:n

% if 2*abs(A(i,i))<=sum(abs(A(i,:)))

% error('Matrica sistema nija dijagonalno dominantna');

% end

% end

% D = diag(diag(A));

% B = D\(D-A);

% c = D\b;

% Bnorm = norm(B,'inf');

% cNorm = norm(c,'inf');

% k = ceil(log((tol*(1-Bnorm)/cNorm))/log(Bnorm))-1

% err = (1-Bnorm)*tol/Bnorm;

% x = c;

% fprintf('Pocetna aproksimacija:')

% disp(x');

% fprintf('Rezultati kroz iteracije:')

% for i = 1:k

% xOld = x;

% x = B*xOld+c;

% disp(x');

% if (norm(abs(x-xOld),'inf')<err)

% break;

% end

% end

% end

A = [10 -1 2; 3 -12 2; -1 3 8];

b = [13 17 4]';

x = Jacobi(A,b,1e-2)

n = 3
k = 8
Pocetna aproksimacija:
    1.3000   -1.4167    0.5000
Rezultati kroz iteracije:
    1.0583   -1.0083    1.1938
    0.9604   -0.9531    1.0104
    1.0026   -1.0082    0.9775
    1.0037   -1.0031    1.0034
    0.9990   -0.9985    1.0016
x =
0.9990 -0.9985 1.0016

Metoda Gaus-Zajdela:

Ova metoda predstavlja modifikaciju Jakobijeve metode:

, gde su i donje i gornje trougaona matrica koje odgovaraju matrici .

% function x = GaussSeidel(A,b,k)

% n = length(A);

% I = eye(n);

% D = diag(diag(A));

% B = D\(D-A);

% c = D\b;

% B1 = tril(B,-1);

% B2 = triu(B,1);

% Binv = inv(I-B1);

% x = c;

% fprintf('Pocetna aproksimacija:')

% disp(x');

% fprintf('Rezultati kroz iteracije:')

% for i = 1:k

% xOld = x;

% x = Binv*(B2*xOld+c);

% disp(x');

% end

% end

A = [3 0.5 -2; -0.2 2 -0.9; -2 -2.5 5];

b = [-1.5 -0.9 -0.5]';

k = 10;

x = GaussSeidel(A,b,n)

Pocetna aproksimacija:
   -0.5000   -0.4500   -0.1000
Rezultati kroz iteracije:
   -0.4917   -0.5442   -0.5687
   -0.7885   -0.7848   -0.8078
   -0.9077   -0.9043   -0.9152
   -0.9594   -0.9578   -0.9627
   -0.9821   -0.9814   -0.9836
   -0.9921   -0.9918   -0.9928
   -0.9965   -0.9964   -0.9968
   -0.9985   -0.9984   -0.9986
   -0.9993   -0.9993   -0.9994
   -0.9997   -0.9997   -0.9997
x =
-0.9997 -0.9997 -0.9997

Direktne metode

Direktne metode predstavljaju klasu metoda koje nakon konačnog skupa koraka daju rešenje sistema (odnosno njegovu aproksimaciju) bez nekog naknadnog poboljšanja.

Podsetimo se Gausove metode eliminacije. Neka nam je dat sistem koji ima jedinstveno rešenje:

Pretpostavimo da je . Tada ćemo markirati nepoznatu , koju onda možemo da eliminišemo iz narednih jednačina tako što drugoj, trećoj,...,-toj jednačini dodamo prvu jednačinu pomnoženu sa , redom. Ako je matrica sistema, tada ćemo dodavanje prve vrste pomnožene pomenutim koeficijentima ostalim vrstama dobiti sledećim množenjem matrica:

gde je i .

U slučaju da je , tada možemo zameniti prvu jednačinu sa nekom od narednih jednačina kod koje je koeficijent uz prvu nepoznatu različit od 0. Nije teško uveriti se da se zamena dve vrste u matrici, -te i -te, može dobiti na sledeći način:

(Malo pojašnjenje: Matrica P ima jedinice na pozicijama i i na glavnoj dijagonali, sem na pozicijama i . Svi ostali elementi u matrici su 0.)

Nastavljajući ovaj postupak, imaćemo niz markiranih promenljih u svakoj od narednih jednačina koje ćemo eliminisati iz jednačina koje dolaze nakon tekuće, sve dok ne stignemo do sistema u trougaonom obliku:

Rešenja sistema sada dobijamo ,,hodom unazad'':

Metoda koja će biti implementirana u nastavku je vrlo slična opisanom Gausovom algoritmu, uz sledeću modifikaciju: Pretpostavimo da smo stigli do vrste u našem algoritmu eliminacije. Tada ćemo tražiti najveću od vrednosti (dakle, razmatramo elemente u -toj koloni počev od ), i po potrebi zameniti mesta tekuće -te vrste i vrste u kojoj se nalazi traženi maksimum.

% function x = Gauss(A,b)

%

% [n, m] = size(A);

% if n~=m

% error('Matrica nije kvadratna!');

% end

%

% for i=1:n

% [a, ind] = max(abs(A(i:end,i))); %interesuje nas samo indeks maksimalnog

% ind=ind+i-1; %prilagodjavamo dimenziji matrice jer je ind indeks maksimalnog

% %gledajuci vektor dobijen od elemenata od i-te do poslednje vrste

% P = eye(n);

% L = eye(n);

% if (i~=ind) %ako ima potrebe za zamenom vrsta

% P(i,i)=0;

% P(ind,ind)=0;

% P(i,ind)=1;

% P(ind,i)=1;

% A = P*A;

% b = P*b;

% end

% L(i+1:end,i) = -(A(i+1:end,i)/A(i,i)); %Formiramo matricu transformacija odredjenu i-tom kolonom

% A = L*A;

% b = L*b;

% end

% %Odredjujemo resenje hodom unazad:

% x = zeros(n,1);

% x(n) = b(n)/A(n,n);

% for i = n-1:-1:1

% suma = 0;

% for j = i+1:n

% suma = suma + x(j)*A(i,j);

% end

% x(i) = (b(i)-suma)/A(i,i);

% end

% end

A = [10 -1 2; 3 -12 2; -1 3 8];

b = [13 17 4]';

x = Gauss(A,b)

x =
1 -1 1

Prethodnu metodu možemo da iskoristimo za određivanje inverza matrice. Naime, problem određivanja matrice , takve da je , odnosno

može da se posmatra kao problem rešavanja sistema linearnih jednačina:

% function X = GaussInverse(A)

%

% [n, m] = size(A);

% if n~=m

% error('Matrica nije kvadratna!');

% end

%

% E = eye(n);

% X = zeros(n);

% for i = 1:n

% X(:,i) = Gauss(A,E(:,i));

% end

% end

A = [-2 2 1; -1 0 -2; 4 -3 1];

X = GaussInverse(A)

X =
-6.0000 -5.0000 -4.0000 -7.0000 -6.0000 -5.0000 3.0000 2.0000 2.0000