Algorytm metody Gaussa -Seidla jest prawie taki sam jak w metodzie Jacobi. Różni się podstawianiem najnowszego przybliżenia wartości
zmiennej dla pozostałych zmiennych w każdym kroku iteracji.
Schemat iteracyjny Gaussa-Seidla rozwiązywania układu równań:
\(\displaystyle{ \begin{cases} 7x_{1}-2x_{2} + x_{3} \ \ \ \ \ \ \ \ = 17 \\ x_{1} - 8x_{2} + 3x_{3} -x_{4} =13\\ 2x_{1} \ \ \ \ \ \ +10x_{3} +x_{4}= 15\\ x_{1} - x_{2} +x_{3} + x_{4} = \ \ 10, \end{cases}}\)
jest następujący:
\(\displaystyle{ x^{(k+1)}_{1} = \frac{(17 +2x^{(k)}_{2} -x^{(k)}_{3})}{7}, }\)
\(\displaystyle{ x^{(k+1)}_{2} = \frac{(-13 +x^{(k+1)}_{1}+3x^{(k)}_{3}-x^{(k)}_{4})}{9}, }\)
\(\displaystyle{ x^{(k+1)}_{3}= \frac{(15 - 2x^{(k+1)}_{1} -x^{(k)}_{4})}{10}, }\)
\(\displaystyle{ x^{(k+1)}_{4} = \frac{(10 -x^{(k+1)}_{1}+x^{(k+1)}_{2}-x^{(k+1)}_{3})}{6}.}\)
Podstawiając \(\displaystyle{ \vec{x}^{(0)} = [0, 0, 0, 0] }\) po prawej stronie każdego z równań, otrzymujemy
\(\displaystyle{ x^{(1)}_{1} = \frac{(17 +2\cdot 0 -0)}{7}= 2,428571, }\)
\(\displaystyle{ x^{(1)}_{2} = \frac{(-13 + 2,428571 +3\cdot 0- 0)}{9}=-1,174603, }\)
\(\displaystyle{ x^{(1)}_{3}= \frac{(15 - 2\cdot 2,428571 - 0)}{10}= 1.014286, }\)
\(\displaystyle{ x^{(1)}_{4} = \frac{(10 -2,428571 -1,174603 -1,014286)}{6}= 0,897090.}\)
Metoda ta jest szybciej zbieżna od metody Jacobi. Nie będziemy przeprowadzać analizy zbieżności tych dwóch metod.
Dokładną analizę można znaleźć na przykład w podręczniku: DAVID KINCAID WARD CHENEY analiza numeryczna. WNT Warszawa 2006.
Dodam tylko, że jeżeli macierz \(\displaystyle{ A }\) jest macierzą ściśle dominującą. tzn:
\(\displaystyle{ |a_{ii}|> \sum_{j=1\\j\neq i}^{n} |a_{ij}|, \ \ i=1,2,...,n, }\)
to metody Jacobi i Gaussa-Seidla są zbieżne niezależnie od wyboru \(\displaystyle{ \vec{x}^{(0)}.}\)
Program w MATLAB
Kod: Zaznacz cały
function seidel(A,b,x0,tol,itmax)
% Rozwiązanie układu równań Ax=b metodą Gaussa-Seidla.
n=length(b);
x=zeros(n,1);
fprintf('\n');
disp(' Macierz rozszerzona =')
Augm=[A b]
Y=zeros(n,1);
Y=x0;
for k=1:itmax+1
for i=1:n
S=0;
for j=1:i-1
S=S+A(i,j)*x(j);
end
for j=i+1:n
S=S+A(i,j)*x0(j);
end
if(A(i,i)==0)
break
end
x(i)=(-S+b(i))/A(i,i);
end
err=abs(norm(x-x0));
rerr=err/(norm(x)+eps);
x0=x;
Y=[Y x];
if(rerr<tol)
break
end
end
%Drukowanie rezultatów obliczeń
if(A(i,i)==0)
disp(' Dzielenie przez zero')
elseif(k==itmax+1)
disp(' Proces iteracyjny nie jest zbieżny')
else
fprintf('\n');
disp('Wektorami rozwiązań są')
fprintf('\n');
disp('iter # 0 1 2 3 4 ...')
fprintf('\n');
for i=1:n
fprintf('x%1.0f= ',i)
fprintf('%10.6f',Y(i,[1:k+1]))
fprintf('\n');
end
fprintf('\n')
disp(['Metoda jest zbieżna po' ,num2str(k), 'iteracjach do']);
x
end
Kod: Zaznacz cały
>> A=[7 -2 1 0;1 -9 3 -1;2 0 10 1;1 -1 1 6];
>> b=[17 13 15 10]';
>> x0=[0 0 0 0]';
>> seidel(A,b,x0,10^(-3),30)
Macierz rozszerzona =
Augm =
7 -2 1 0 17
1 -9 3 -1 13
2 0 10 1 15
1 -1 1 6 10
Wektorami rozwiązań są
iter # 0 1 2 3 4 ...
x1= 0.000000 2.428571 1.948073 2.000025 2.001849 2.000025
x2= 0.000000 -1.174603 -0.989573 -0.993877 -1.000154 -1.000130
x3= 0.000000 1.014286 1.020676 0.999300 0.999517 1.000020
x4= 0.000000 0.897090 1.006946 1.001133 0.999747 0.999971
Metoda jest zbieżna po 5 iteracjach do
x =
2.0000
-1.0001
1.0000
1.0000