Iteracyjne metody rozwiązywania układu równań liniowych:
\(\displaystyle{ A\cdot \vec{x} = \vec{b} }\)
polegają na przekształceniu układu równań do postaci:
\(\displaystyle{ \vec{x} = \vec{b}' - B\cdot \vec{x} }\)
i po przyjęciu rozwiązania początkowego \(\displaystyle{ \vec{x}_{0} }\)
tworzeniu ciągu rozwiązań aproksymujących:
\(\displaystyle{ \vec{x}^{(k)} = \vec{b}' - B\cdot \vec{x}^{(k-1)} }\)
Jako przykład wykonamy jedną iterację metodą Jacobi następującego 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}}\)
przyjmując rozwiązanie początkowe \(\displaystyle{ \vec{x}^{(0)} = \vec{0}.}\)
Układ równań przekształcamy do postaci:
\(\displaystyle{ x_{1} = \frac{(17 +2x_{2} -x_{3})}{7}, }\)
\(\displaystyle{ x_{2} = \frac{(-13 + x_{1}+3x_{3} - x_{4})}{9}, }\)
\(\displaystyle{ x_{3}= \frac{(15 - 2x_{1} - x_{4})}{10}, }\)
\(\displaystyle{ x_{4} = \frac{(10 - x_{1} + x_{2} - x_{3})}{6}.}\)
i zapisujemy w postaci schematu iteracyjnego Jacobi:
\(\displaystyle{ x^{(k+1)}_{1} = \frac{(17 +2x^{(k)}_{2} -x^{(k)}_{3})}{7}, }\)
\(\displaystyle{ x^{(k+1)}_{2} = \frac{(-13 +x^{(k)}_{1}+3x^{(k)}_{3}-x^{(k)}_{4})}{9}, }\)
\(\displaystyle{ x^{(k+1)}_{3}= \frac{(15 - 2x^{(k)}_{1} -x^{(k)}_{4})}{10}, }\)
\(\displaystyle{ x^{(k+1)}_{4} = \frac{(10 -x^{(k)}_{1}+x^{(k)}_{2}-x^{(k)}_{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 +0 +3\cdot 0- 0)}{9}= -1,444444, }\)
\(\displaystyle{ x^{(1)}_{3}= \frac{(15 - 2\cdot 0 - 0)}{10}= 1.500000, }\)
\(\displaystyle{ x^{(1)}_{4} = \frac{(10 -0 + 0 - 0)}{6}= 1.666667.}\)
Program w MATLAB
Kod: Zaznacz cały
function jacobi(A,b,x0,tol,itmax)
% Rozwiązanie układu równań Ax=b metodą Jacobi.
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:n
if(j~=i)
S=S+A(i,j)*x0(j);
end
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
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]';
>> jacobi(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.801587 2.061829 1.977515 2.008259 1.997187 2.001020 1.999639 2.000127
x2= 0.000000 -1.444444 -0.859788 -1.047413 -0.981367 -1.006166 -0.997771 -1.000802 -0.999721 -1.000100
x3= 0.000000 1.500000 0.847619 1.062566 0.979451 1.007360 0.997320 1.000926 0.999667 1.000118
x4= 0.000000 1.666667 0.771164 1.081834 0.971365 1.010278 0.996369 1.001287 0.999542 1.000162
Metoda jest zbieżna po 9 iteracjach do
x =
2.0001
-1.0001
1.0001
1.0002