Praktická cvičení z Numerické matematiky

Soustavy lineárních a nelineárních algebraických rovnic - 3. týden

Jacobiho iterační metoda v MATLABu

Tento zapis je pro matici 3x3 (plnou), v praxi se ale uziva pro matice specialni (ridke) se znamou strukturou. Zapis se tedy vyrazne lisi.

  A = [4 -1 -2; -1 3 -1; -2 -1 4];
  b = [2 ; 0;  2];
  
  % pocatecni podminka napr. prava strana
  X = b; 
  n = size(A, 1);

  % jedna iterace Jacobi (opakujeme)
  X1 = X;
  for i = 1:n, X1(i) = X(i) + (b(i) - A(i,:) * X) / A(i,i); end
  % spocteme rozdil a priradime do X
  norm(X1 - X, 2), X = X1;

Gaussova-Seidelova iterační metoda v MATLABu

Plati obdobna poznamka jako u Jacobiho metody.

  A = [4 -1 -2; -1 3 -1; -2 -1 4];
  b = [2 ; 0;  2];
  
  % pocatecni podminka napr. prava strana
  X = b; 
  n = size(A, 1);

  % jedna iterace GS (Xold - predchozi iterace), opakujeme
  Xold =X; 
  for i = 1:n, X(i) = X(i) + (b(i) - A(i,:) * X) / A(i,i); end, 
  norm(X - Xold, inf)

Maticový zápis Jacobiho iterační metody v MATLABu

Zapis pomoci matic U. Toto pouziti je zavadejici (!), v praxi se to takto nikdy (!) nepouziva. Muze ale slouzit pro overeni podminek konvergence z prikladu na cviceni.

  A = [4 -1 -2; -1 3 -1; -2 -1 4];
  b = [2 ; 0;  2];
  
  D = diag(diag(A)); L = tril(A, -1); P = triu(A, 1);

  U  = -D \ ( L + P );  V = D \ b;
  % pocatecni priblizeni
  X = V;
  % opakujeme
  X = U * X + V

Maticový zápis Gauss-Seidelovy iterační metody v MATLABu

Zapis pomoci matic U. Toto pouziti je zavadejici (!), v praxi se to takto nikdy (!) nepouziva. Muze ale slouzit pro overeni podminek konvergence z prikladu na cviceni.

  A = [4 -1 -2; -1 3 -1; -2 -1 4];
  b = [2 ; 0;  2];
  
  D = diag(diag(A)); L = tril(A, -1); P = triu(A, 1);

  U = -(D + L)\ P;    V = (D + L) \b;

  % pocatecni priblizeni
  X = V;
  % opakujeme
  X = U * X + V

Příklad 3.7: Jacobiho iterační metoda (Program v MATLABu)

Program pro řešení třídiagonální soustavy rovnic z příkladu 3.7. Matici vůbec neukládáme, jen provádíme výpočet. Demonstrace v MATLABu, jen pro snadnější pochopení. Indexy jsou posunuty o jedna, MATLAB totiz nepracuje s indexem 0.

n = 20; h = 1. / (n + 1); 
b = h * h * ones(n + 2, 1);
b(1) = 0; b(end) = 0;
% hledame reseni A * x = b, kde A vubec neukladame
X = b; XN = X;

% Jacobi: vypocet nove iterace XN z X, opakujeme:
for i = 2:(n+1), XN(i)=(b(i) + X(i-1) + X(i+1)) / 2; end

norm(XN - X)

X = XN;

% opakujeme (uschovame starou iteraci)
Xold = X;
for i = 2:(n+1), X(i)=(b(i) + X(i-1) + X(i+1)) / 2;  end,
norm(X - Xold)

Příklad 3.7: Jacobiho iterační metoda (Program v C)

Program pro řešení soustavy rovnic z příkladu 3.7. Matici vůbec neukládáme, jen provádíme výpočet.

jacobi.c

Funkce pro iterace:

void jacobi(int n, double *Xnew, double *X, double *b, int niter)
{
   int i, iter;

   for(iter = 0; iter < niter; iter++)
   {
     for(i = 1; i < n; i++)
         Xnew[i] =  (B[i] + X[i - 1] + X[i + 1]) /  2.;
   }
}

Příklad 3.7: Gauss-Seidelova iterační metoda (Program v C)

Program pro řešení soustavy rovnic z příkladu 3.7. Matici vůbec neukládáme, jen provádíme výpočet.

gaussseidel.c

double gs(int n, double *x, double *b, int niter)
{
   int i, iter;
   for(iter = 0; iter < niter; iter++)
   {
       for(i = 1; i < n; i++)
           x[i] =  (b[i] + x[i - 1] + x[i + 1]) / 2;
   }
   return rez;
}

sor.c

Příklad 3.8: Jacobi/Gauss-Seidelova iterační metoda (Program v C)

Program pro řešení soustavy rovnic z příkladu 3.8. Matici vůbec neukládáme, jen provádíme výpočet.

jac2d.c

gs2d.c

Část kódu, která provádí řešení (Vektory jsou ulozeny jako dvourozmerne pole U[][] a B[][] , globální proměnné)


void gs2d(int n, int niter)
{
 int i, j, iter;
 for(iter = 0; iter < niter; iter++)
 {
  for(i = 1; i < n; i++)
   for(j = 1; j < n; j++)
    U[i][j] = (B[i][j] - U[i - 1][j]  - U[i + 1][j]
                       - U[i][j - 1]  - U[i][j + 1]) / 4;
 }
}

Zobrazeni vysledku v MATLABu

load sol.dat;
mesh(sol);

Zobrazeni vysledku v GNUPLOTu

sp "sol.dat" matrix with linepoints