Sistemi linearnih jednačina

16. zadatak

• Napisati M-fajl $sistem.m$ sa funkcijom $x = sistem(A, b)$ koja metodom proste iteracije rešava sistem jednačina $Ax = b$. Broj iteracija fiksirati na 50.

  $sistem.m$

  x
  function x = sistem(A,b)
  ​
  n = length(A); % Kada se length() primeni na matricu, vraća broj vrsti 
  A1 = zeros(n);
  b1 = zeros(n,1);
  ​
  % Potrebno je transformisati sistem Ax = b u oblik x = A1*x + b1
  for i = 1:n
      A1(i,:) = -A(i,:) / A(i,i); % i-tu vrstu delimo sa -A(i,i)
      A1(i,i) = 0;                % na dijagonali su 0
      b1(i) = b(i) / A(i,i);      % i-ti el. vektora b delimo sa A(i,i)
  end
  disp('A1 =')
  disp(A1)
  disp('b1 =')
  disp(b1)
  ​
  % Početna aproksimacija rešenja je vektor b1
  x = b1;
  ​
  % Za dobijanje nove aproksimacije koristi se rekurentna formula: x_n = A1*x_n-1 + b1
  for i = 1:50
      x= A1 * x + b1;
  end
  

   

• Napisati M-fajl $matrica.m$ sa funkcijom $[A, b, x] = matrica(broj, d)$ koja vraća kolonu $b$ duzine $d$ čiji su svi elementi jedinice, kvadratnu matricu $A_{d×d}$ koja iznad dijagonale ima jedinice, po dijagonali ima $10$ · $broj$, dok na prvoj poddijagonali ima vrednost $broj − 1$, na drugoj poddijagonali ima vrednost $broj − 2$, itd. kao i vektor $x$ koji je rešenje sistema $Ax = b$ (koristiti funkciju $sistem()$ za nalaženje vektora $x$ ).

  $matrica.m$

  x
  function [A, b, x] = matrica(broj, d)
  ​
  b = ones(d, 1); % vektor kolona dužine d sa svim jedinicama
  ad = broj*10*ones(d,1);
  ​
  % Formiramo matricu koja na glavnoj dijagonali ima vektor ad
  A = diag(ad);
  ​
  % Poziv diag(ad, 1) formiraće matricu koja na prvoj naddijagonali ima vektor ad 
  % Poziv diag(ad, -1) formiraće matricu koja na prvoj poddijagonali ima vektor ad 
  ​
  for i = 1:d-1
      ad = ones(d-i, 1);
      % popunjavamo matricu ispod (II sabirak) i iznad (I sabirak) glavne dijagonale
      A = A + diag(ad*(broj-i),-1*i) + diag(ad,i); 
  end
  ​
  disp('A =')
  disp(A)
  x = sistem(A, b);
  disp('Resenje sistema je: ')
  disp(x)
  %Tačno rešenje sistema se može dobiti naredbom x = A\b
  ​
  

17. zadatak

• Formirati M-fajl $dominantna.m$ sa funkcijom $d = dominantna(A)$ koja proverava da li je zadata matrica A dijagonalno dominantna. Funkcija vraća vrednost $1$ ako je matrica dijagonalno dominantna, inače vraća $0$.

  $dominantna.m$

  x
  function d = dominantna(A)
  ​
  n = length(A);
  d = zeros(n, 1); % el. d(i) nosiće informaciju o ispunjenosti uslova u i-toj vrsti
  ​
  for i = 1:n
      % Proveravamo da li je suma apsolutnih vrednosti vandijagonalnih el. i-te vrste 
      % manja od apsolutne vrednosti el. na dijagonali
      if sum(abs(A(i, :))) - abs(A(i, i)) <  abs(A(i, i)) 
          d(i) = 1; 
      else
          d(i) = 0;
      end
  end
  d = all(d==1);
  % ili d = all(d);
  

   

• Formirati M-fajl $sistem.m$ sa funkcijom $[iter, x] = sistem(A, b, tol)$ koja nalazi rešenje sistema $Ax = b$ Gaus-Zajdelovom metodom pod uslovom da je matrica $A$ dijagonalno dominantna. Inače ispisati poruku: "Matrica nije dijagonalno dominantna". Iterativni postupak se prekida kada za dve uzastopne iteracije važi $|x_k − x_{k−1}| ≤ tol$. Program vraća rešenje $x$ i broj iteracija $iter$.

  $sistem.m$

  x
  function [iter, x] = sistem2(A, b, tol)
  ​
  % Prekidamo izvršavanje programa ukoliko matrica nije dijagonalno dominantna
  if dominantna(A)==0
      error('Matrica nije dijagonalno dominantna');  
  end
  ​
  n = length(A);
  A1 = zeros(n);
  b1 = zeros(n, 1);
  ​
  ​
  for i = 1:n
      A1(i, :) = -A(i, :) / A(i, i);
      A1(i, i) = 0;
      b1(i) = b(i) / A(i, i);
  end
  ​
  q = norm(A1, inf); % uniformna norma matrice  
  tol = tol * (1-q) / q;
  ​
  x0 = b1; % početna aproksimacija rešenja
  x = x0;
  iter = 1;
  ​
  while 1 
       % za nalaženje nove apoksimacije koordinate x_i tačnog rešenja
       % koristimo vrednosti koordianta x_j, j=1,...i-1 dobijenih iz iste iteracije
       % zbog toga moramo koristiti for petlju i u njenom i-tom prolasku 
       % računati vrednost koordinate x_i (razlika u odnosu na metodu proste iteracije) 
       for i = 1:n
           x(i) = A1(i, :) * x + b1(i); 
       end
       iter = iter + 1;
       if all(abs(x - x0) <= tol) % kriterijum zaustavljanja
         break 
       end
       x0 = x;
  end
  ​
  % Provera:
  % disp('Tacno resenje:');
  % disp(A\b);
  

18. zadatak

• Formirati M-fajl $LUdekompozicija.m$ sa funkcijom $x = LUdekompozicija(A, b)$ koja metodom LU dekompozicije vraća rešenje sistema $Ax=b$. Koristiti ugrađenu funkciju $lu()$.

  $LUdekompozicija.m$

  x
  function X = LUdekompozicija(A, b)
  ​
  % Rešavamo sistem Ax = b metodom LU dekompozicije
  % Na primer:
  % A = [24.21 2.42 3.85; 2.31 31.49 1.52; 3.49 4.85 28.72]; 
  % b = [1 0 0]';
  ​
  n = length(b);
  % Koristimo ugrađenu funkciju lu() koja vraća matrice L, U i P
  [L U P] = lu(A);
  ​
  disp('L =')
  disp(L)
  disp('U =')
  disp(U);
  ​
  % Sada se problem Ax = b svodi na rešavanje dva nova sistema
  % L*y = b
  % ----------------------------------------
  % L(1,1)*y(1)                                      = b(1)
  % L(2,1)*y(1) + L(2,2)*y(2)                        = b(2)
  % L(2,1)*y(1) + L(2,2)*y(2) + L(2,3)*y(3)          = b(2)
  % .....
  % L(n,1)*y(1) + L(n,2)*y(2)+    ...   + L(n,n)*y(n)= b(n)
  % ----------------------------------------
  % U*x = y
  % ----------------------------------------
  % U(1,1)*x(1) + U(1,2)*x(2) +   ...   + U(1,n)*x(n) = y(1)
  %               U(2,2)*x(2) +   ...   + U(2,n)*x(n) = y(2)
  %                         .....
  %                 U(n-1,n-1)*x(n-1) + U(n-1,n)*x(n) = y(n-1)
  %                                       U(n,n)*x(n) = y(n)
  % ----------------------------------------
  %%%%%%%%%%%%%% rešavamo sistem L*y = b %%%%%%%%%%%%%%
  y = zeros(n, 1);
  y(1) = b(1);
  ​
  for i = 2:n
      y(i) = (b(i) - (L(i, :) * y));
  end
  ​
  disp('y =')
  disp(y)
  %%%%%%%%%%%%%% rešavamo sistem U*x = y %%%%%%%%%%%%%%
  X = zeros(n, 1);
  X(n) = y(n) / U(n, n);
  ​
  for i = n-1:-1:1
      X(i) = (y(i) - U(i, :) * X) / U(i, i);
  end
  

  $LUdekomp.m$

  x
  function x=LUdekomp(A,b)
  ​
  % kod za LU dekompoziciju (nije ugrađena funkcija, jedinice po dijagonali kod U)
  ​
  n = length(A);
  L = zeros(n);
  U = eye(n);
  ​
  L(:, 1) = A(:, 1);
  U(1, 2:n) = A(1, 2:n) / L(1, 1);
  ​
  for i = 1:n
      for j = 1:n
          if j <= i
              L(i, j) = A(i, j) - (L(i, 1:(j-1)) * U(1:(j-1), j));
          else 
              U(i, j) = (A(i, j)-L(i, 1:i-1) * U(1:i-1, j)) / L(i, i); 
          end
      end
  end
  disp('L =');
  disp(L);
  disp('U =');
  disp(U);
  ​
  %%%%%%%%%%%%%% rešavamo sistem L*y = b %%%%%%%%%%%%%%
  y = zeros(n, 1);
  y(1) = b(1) / L(1, 1);
  for i = 2:n
      y(i) = (b(i) - (L(i, 1:i-1) * y(1:i-1))) / L(i, i);
  end
  disp('y=')
  disp(y)
  %%%%%%%%%%%%%% rešavamo sistem U*x = y %%%%%%%%%%%%%%
  x=zeros(n,1);
  x(n)=y(n);
  for j=n-1:-1:1
      x(j)=y(j)-(U(j,j+1:n)*x(j+1:n));
  end
  

• Formirati M-fajl $inverzna.m$ sa funkcijom $inverzna = inverz(A)$ koja nalazi matricu $A^{-1}$ korišćenjem funkcije iz fajla $LUdekompozicija.m$.

  $inverzna.m$

  x
  % primer: A = [7.05 0.11 0.13 0.15; 
  %               0.17 7.41 0.75 0.80; 
  %               0.97 1.02 8.11 0.11; 
  %               0.13 0.20 0.29 8.42];
  ​
  function inverzna = inverz(A)
  ​
  % dužina matrice A
  [n, m] = size(A);
  ​
  % E = jedinična matrica dimenzije nxn
  E = eye(n);   
  ​
  inverzna = zeros(n);
  ​
  for i = 1 : n
      % rešavamo n sistema: A*X(i) = E(i)
      % inverzna = [X(1) X(2) .. X(n)]
      % X(i) je i-ta kolona inverzne matrice 
      % E(i) je i-ta kolona jedinične matrice
      inverzna(:, i) = LUdekompozicija(A, E(:, i));
  end