Základy algoritmizace a programování

Petr Sváček, Luděk Beneš, Olga Majlingová
Ústav technické matematiky

Přednáška č. 12

MATLAB: Soustavy lin. rovnic, prime metody. Iteracni metody.

Stručně resumé přednášky (zapsané v MATLABU)

  1. Gaussova eliminace (krok po kroku)
    Soustava rovnic
    Definujeme matici A a pravou stranu b. 
    
A=[5 2 0 3; 2 6 1 1; 1 0 3 2; 1 -2 2 5]
b =[ 14; 1; 9; 16 ]
    Upravy budeme delat v rozsirene matici AA=(A,b) 
    
AA=[A,b];
n=size(A,1);
    Dopredny chod Gaussovy eliminace (uprava na horni trojuhelnikovy tvar, bez pivotace)
    Vezmeme prvni radek a pomoci nej eliminujeme prvky v prvnim sloupci a 2. az 4. radku:
    
mult=AA(2,1)/AA(1,1);
AA(2,:)=AA(2,:)-mult*AA(1,:)

mult=AA(3,1)/AA(1,1);
AA(3,:)=AA(3,:)-mult*AA(1,:)

mult=AA(4,1)/AA(1,1);
AA(4,:)=AA(4,:)-mult*AA(1,:)
    Zapis pomoci indexu i a j (zacneme znova) 
    
AA=[A ,b];

% Eliminujeme v 1. sloupci a j-tem radku
i=1; 
j=2 
mult=AA(j,i)/AA(i,i);
AA(j,:)=AA(j,:)-mult*AA(i,:)

j=3 
mult=AA(j,i)/AA(i,i);
AA(j,:)=AA(j,:)-mult*AA(i,:)

j=4 
mult=AA(j,i)/AA(i,i);
AA(j,:)=AA(j,:)-mult*AA(i,:)
    Pokracujeme s 2. sloupcem 
    
% Eliminujeme ve 2. sloupci a j-tem radku
i=2; 

j=3 
mult=AA(j,i)/AA(i,i);
AA(j,:)=AA(j,:)-mult*AA(i,:)
j=4 
mult=AA(j,i)/AA(i,i);
AA(j,:)=AA(j,:)-mult*AA(i,:)
    Pokracujeme s 3. sloupcem 
    
% Eliminujeme ve 3. sloupci a j-tem radku
i=3; 
j=4
mult=AA(j,i)/AA(i,i);
AA(j,:)=AA(j,:)-mult*AA(i,:)
    Zpetny chod Gaussovy eliminace 
    
x=zeros(size(bb));
i=4; 

bb=AA(:,5);


i=4
x(i)=(bb(i)-AA(i,(i+1):n)*x((i+1):n))/AA(i,i)
i=3
x(i)=(bb(i)-AA(i,(i+1):n)*x((i+1):n))/AA(i,i)
i=2
x(i)=(bb(i)-AA(i,(i+1):n)*x((i+1):n))/AA(i,i)
i=1
x(i)=(bb(i)-AA(i,(i+1):n)*x((i+1):n))/AA(i,i)
    Zkouska 
    A*x
b
A*x-b
  2. Gaussova eliminace (algoritmizace)
    Dopredny chod
    Zapiseme dopredny chod pomoci cyklu (nejprve v j) 
    
AA=[A ,b];

n=size(A,1);

% Eliminujeme v 1. sloupci a j-tem radku
i=1; 
for j=(i+1):n, 
  mult=AA(j,i)/AA(i,i);
  AA(j,:)=AA(j,:)-mult*AA(i,:);
end
    a nasledne cyklus pres i 
    
AA=[A ,b];

n=size(A,1);

% Eliminujeme v i-tem sloupci a j-tem radku
for i=1:n-1, 
  for j=(i+1):n, 
    mult=AA(j,i)/AA(i,i);
    AA(j,:)=AA(j,:)-mult*AA(i,:);
  end
end
    Zpetny chod
    x=zeros(size(b));
bb=AA(:,n+1);
for i=n:-1:1,
  x(i)=(bb(i)-AA(i,(i+1):n)*x((i+1):n))/AA(i,i);
end
    Vse ve funkci
    Vytvorte funkci gaussova_elim.m 
    function x=gaussova_elim(A,b)
%
% Funkce pro Gaussovu eliminaci.
%

AA=[A ,b];

n=size(A,1);

% Eliminujeme v i-tem sloupci a j-tem radku
for i=1:n-1, 
% PIVOTACE  (DOPLNIT)
  for j=(i+1):n, 
    mult=AA(j,i)/AA(i,i);
    AA(j,:)=AA(j,:)-mult*AA(i,:);
% i,j,AA, pause ?
  end
end

x=zeros(size(b));
bb=AA(:,n+1);
for i=n:-1:1,
  x(i)=(bb(i)-AA(i,(i+1):n)*x((i+1):n))/AA(i,i);
end
  3. Vyzkousime na matici 5x5
    A=[5 2 0 3 -1 ; 2 6 1 1 0; 1 0 3 2 2 ; 1 -2 2 5 1 ; 1 0 0 0 2]

b=[ 0 ; -15;   25;  -10;   30]

x=gaussova_elim(A,b);
  4. Gaussova eliminace s pivotaci (doplnte) Zvolte matici 
      A=[1 2 0 3; 2 4 1 1; 0 1 3 2; 1 4 2 5]

  b=[1;7;-1;1]
    a zkuste 
    x=gaussova_elim(A,b);
    Vysledek vede k deleni nulou! (Proc? Vyzkousejte odkomentovat pause )

    Chyba nastava pri eliminaci ve druhem sloupci -> je treba prohodit do i-teho radku, takovy, ktery nema nulu v i-tem sloupci. Nahradime radek % PIVOTACE prikazy, pro nalezeni nejvetsiho prvku (v absolutni hodnote) v i tem radku a i-tem az n-tem sloupci. 
      val=0; k=i;
  for j=i:n,
    if (abs(AA(j,i))>val),
      val=abs(AA(j,i));
      k=j;
    end
  end
  % prohodime i-ty a k-ty radek
  pom=AA(i,:); AA(i,:)=AA(k,:); AA(k,:)=pom;
  5. Reseni soustavy rovnic pomoci operatoru \(slash) Zvolte matici 
      A=[1 2 0 3; 2 4 1 1; 0 1 3 2; 1 4 2 5]

  b=[1;7;-1;1]
    a reseni lze ziskat 
     x=A\b;
    V pripade regularni matice existuje inverzni matice (inv(A)) a "teoreticky" lze reseni vypocitat pomoci ni. Tento pristup je ale nevhodny, zvysuje numerickou chybu resni (diky zaokrouhlovacim chybam).

    Navic inverzni matici MATLAB spocita i pro matici singularni: 

    A=[1 2 3 ; 4 5 6; 7 8 9]

inv(A)

b=[2;3;4]

% srovnej

inv(A)*b

A\b
    POZOR: reseni v tomto pripade neni jednoznacne urceno, a neni vhodne hledat ho prostrednictvim pocitace. Viz nasledujici pripad: 
    b=[1;-2;1]

A\b

inv(A)*b
  6. Jacobiho iteracni metoda
    Soubor jacobi.m 
    function [x,iter,rez]=jacobi(A,b,TOL,MXITER)
%
%  function [x,iter,rez]=gas(A,b,TOL,MXITER)
%    Jacobiho iteracni metod
%
%
U=triu(A,1);
L=tril(A,-1);
D=diag(diag(A));
x=b;

err=1;
iter=1;
while (err>TOL & iter<MXITER)
  xn=D\(-(L+U)*x+b);
  reziduum=norm(A*xn-b);
  err=norm(xn-x)
  x=xn;
  rez(iter)=reziduum;
  iter=iter +1;
end
  7. Gauss-Seidelova iteracni metoda.
    Soubor gauss_seidel.m. 
    
function [x,iter,rez]=gauss_seidel(A,b,TOL,MXITER)
%
%  function [x,iter,rez]=gauss_seidel(A,b,TOL,MXITER)
%    Gauss-Seidelova iteracni metod
%    TOL - tolerance
%    MXITER - maximalni pocet iteraci
U=triu(A,1);
L=tril(A,-1);
D=diag(diag(A));
x=b;

err=1;
iter=1;
while (err > TOL & iter < MXITER)
  xn=(L+D)\(-U*x+b);
  err=norm(xn-x);
  reziduum=norm(A*xn-b);
  x=xn;
  rez(iter)=reziduum;
  iter=iter +1;
end

Úlohy na procvičení

Soustava rovnic s Hilbertovou matici. 
A=hilb(10);
xx=rand(10,1);
b=A*xx;

x1=A\b;
x2=inv(A)*b;
norm(A*x1-b)

norm(A*x2-b)

norm(x1-x2)

norm(x1-xx)

norm(x2-xx)