Rozwiązujemy układ równań liniowych:
\(\displaystyle{ A\cdot\vec{x} = \vec{b} }\)
Załóżmy jego postać początkową:
\(\displaystyle{ \left[ \begin{matrix} a^{(1)}_{11} & a^{(1)}_{12} & ...& a^{(1)}_{1n}\\a^{(1)}_{21} & a^{(1)}_{22} & ...& a^{(1)}_{2n}\\
\vdots & \vdots & ... & \vdots \\
a^{(1)}_{n1} & a^{(1)}_{n2} & ...& a^{(1)}_{nn} \end{matrix} \right] \cdot \left[\begin{matrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{matrix}\right]=
\left[ \begin{matrix} b^{(1)}_{1} \\ b^{(1)}_{2} \\ \vdots \\ b^{(1)}_{n} \end{matrix} \right]. }\)
Eliminacja Gaussa-Jordana polega na stosowaniu przekształceń elementarnych na wierszach macierzy rozszerzonej układu \(\displaystyle{ [A|\vec{b}^{(1)}],}\)
w celu sprowadzenia jej do równoważnej postaci trójkątnej-górnej \(\displaystyle{ U\cdot \vec{x}= \vec{d}. }\)
Opis algortmu.
Krok 1
Zakładamy, że element główny macierzy układu \(\displaystyle{ a_{11}\neq 0 }\)
Definiujemy mnożniki wierszy:
\(\displaystyle{ m_{i1}= \frac{a^{(1)}_{i1}}{a^{(1)}_{11}} }\)
Mnożymy wiersz pierwszy przez \(\displaystyle{ m_{i1} }\) i odejmujemy od \(\displaystyle{ i-}\) tego wiersza \(\displaystyle{ (i=2,3,...,n), }\)
otrzymując
\(\displaystyle{ a^{(2)}_{ij} = a^{(1)}_{ij} - m_{i1}\cdot a^{(1)}_{1j}, \ \ j=2,3,...,n }\)
\(\displaystyle{ b^{(2)}_{i} = b^{(1)}_{i} - m_{i1}\cdot b^{(1)}_{1}.}\)
Pierwszy wiersz macierzy rozszerzonej układu pozostaje bez zmian. W pierwszej kolumnie macierzy \(\displaystyle{ A }\) poniżej elementu \(\displaystyle{
a_{11} }\) otrzymujemy zera.
W rezultacie tego przekształcenia otrzymujemy równoważny układ równań:
\(\displaystyle{ \left[ \begin{matrix} a^{(1)}_{11} & a^{(1)}_{12} & ...& a^{(1)}_{1n}\\ 0 & a^{(2)}_{22} & ...& a^{(2}_{2n}\\
\vdots & \vdots & ... & \vdots \\
0 & a^{(2)}_{n2} & ...& a^{(2)}_{nn} \end{matrix} \right] \cdot \left[\begin{matrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{matrix}\right]=
\left[ \begin{matrix} b^{(1)}_{1} \\ b^{(2)}_{2} \\ \vdots \\ b^{(2)}_{n} \end{matrix} \right]. }\)
Kontynuujemy ten proces transformacji układu
Krok k
Zakładamy, że element macierzy \(\displaystyle{ a^{(k)}_{kk} \neq 0 }\)
Definijemy mnożniki wierszy:
\(\displaystyle{ m_{ik} = \frac{a^{(k)}_{ik}}{a^{(k)}_{kk}} }\)
Mnożymy, wiersz \(\displaystyle{ k-}\) ty przez \(\displaystyle{ m_{ik} }\) i odejmujemy od wiersza \(\displaystyle{ (i=k+1, k+2,...,n),}\)
otrzymując
\(\displaystyle{ a^{(k+1)}_{ij}= a^{(k)}_{ij} - m_{ik}\cdot a^{(k)}_{kj}, \ \ j=k+1,k+2,...,n,}\)
\(\displaystyle{ b^{(k+1)}_{i} = b^{(k)}_{i} - m_{ik}\cdot b^{(k)}_{k}.}\)
W tym kroku pierwsz \(\displaystyle{ k-}\) kolumn macierzy poniżej diagonali zawiera zera i wiersze \(\displaystyle{ 1 }\) do \(\displaystyle{ k }\) są
niezmienione.
Dla \(\displaystyle{ k = n-1,}\) otrzymujemy układ trójkątny - górny:
\(\displaystyle{ \begin{cases} a^{(1)}_{11} x_{1} + a^{(1)}_{12}x_{2}+...+ a_{1n}x_{n} = b^{(1)} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ a^{(2)}_{22}x_{2}+ ...+ a^{(2)}_{2n}x_{n}= b^{(2)}_{2}\\ ...............................\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ a^{(n)}_{nn}x_{n} = b^{(n)}_{n} \end{cases} }\)
Stosując podstawienie wsteczne, otrzymujemy rozwiązanie układu równań:
\(\displaystyle{ x_{n}= \frac{b^{(n)}_{n}}{a^{(n)}_{nn}}, }\)
\(\displaystyle{ x_{n-1} = \frac{b^{(n-1)}_{n-1} -a^{(n-1)}_{n-1,n}\cdot x_{n}}{a^{(n-1)}_{n-1}}, }\)
\(\displaystyle{ x_{i} = \frac{b^{(i)}_{i} - (a^{(i)}_{i,i+1}\cdot x_{i+1} + \ \ ...\ \ +a^{i}_{in}\cdot x_{n})}{a^{(i)}_{ii}} = \frac{b^{(i)} -\sum_{j=i+1}^{n}a^{(i)}_{ij}\cdot x_{j}}{a^{(i)}_{ii}}, \ \ i = n-2, n-3, ..., 1. }\)
Nie będziemy przeprowadzać analizy błędów zaokrągleń ani rachunków kosztów tej metody - odsyłam Państwa do monografii Moich Przyjaciół:
Andrzeja Kiełbasińskiego i Huberta Schwetlicka pt. Numeryczna algebra liniowa. WNT Warszawa 1992.
Dodam tylko, że koszt tej klasycznej metody eliminacji Gaussa- Jordana dla wielkich układów równań liniowych jest rzędu \(\displaystyle{ \frac{n^3}{3} }\) operacji.
Program w OCTAVE
Kod: Zaznacz cały
function ngaussel(A,b)
n=length(b);
x=zeros(n,1);
fprintf('\n');
disp('Macierz roszerzona')
augm=[A b]
for k=1:n-1
for i=k+1:n
m=A(i,k)/A(k,k);
for j=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
A(i,k)=m;
b(i)=b(i)-m*b(k);
end
end
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
S=b(i);
for j=i+1:n
S=S-A(i,j)*x(j);
end
x(i)=S/A(i,i);
end
fprintf('\n');
disp('Macierz trojkątna-górna C =')
fprintf('\n');
for i=1:n
for j=1:n
if (j<i) A(i,j)=0;end
end
end
C=[A,b]
fprintf('\n')
disp('Rozwiązanie =')
x
>> A=[1 1 1 1;2 3 1 5;-1 1 -5 3;3 1 7 -2];
>> b=[10 31 -2 18]';
>> ngaussel(A,b)
Macierz roszerzona
augm =
1 1 1 1 10
2 3 1 5 31
-1 1 -5 3 -2
3 1 7 -2 18
Macierz trojkątna-górna C =
C =
1 1 1 1 10
0 1 -1 3 11
0 0 -2 -2 -14
0 0 0 -1 -4
Rozwiązanie =
x =
1
2
3
4
Kod: Zaznacz cały
>> A=[-1 1 -2 0 4 -2 2
1 4 -2 -4 -6 -2 2
-1 1 6 2 8 2 8
1 -1 -1 -2 6 -4 8
-1 1 1 2 4 -4 2
1 -1 -1 1 1 4 4
-1 1 1 -1 -1 1 4];
>> b=[2 4 6 8 -2 -4 -6]';
>> ngaussel(A,b)
Macierz roszerzona
augm =
-1 1 -2 0 4 -2 2 2
1 4 -2 -4 -6 -2 2 4
-1 1 6 2 8 2 8 6
1 -1 -1 -2 6 -4 8 8
-1 1 1 2 4 -4 2 -2
1 -1 -1 1 1 4 4 -4
-1 1 1 -1 -1 1 4 -6
Macierz trojkątna-górna C =
C =
-1.0000 1.0000 -2.0000 0 4.0000 -2.0000 2.0000 2.0000
0 5.0000 -4.0000 -4.0000 -2.0000 -4.0000 4.0000 6.0000
0 0 8.0000 2.0000 4.0000 4.0000 6.0000 4.0000
0 0 0 -1.2500 11.5000 -4.5000 12.2500 11.5000
0 0 0 0 10.0000 -8.0000 10.0000 6.0000
0 0 0 0 0 15.2800 2.8000 2.0400
0 0 0 0 0 0 7.0838 -10.6675
Rozwiązanie =
x =
4.0109
2.0457
0.9690
-3.0440
2.4335