%% Układy równań liniowych
% Opracował: dr inż Robert JAKUBOWSKI, Politechnika Rzeszowska
% Table of Contents
% Układy równań liniowych
% Przypadki szczególne układów równań
%     Rozwiązywanie układów równań z elementami jedynie na głównej przekątnej
%     Rozwiązywanie układów równań z macierzą trójkątną
% Rozwiązywanie układów równań w ogólnym w ogólnym przypadku
%     Eliminacja Gausa
%     Funkcja do rozwiązywania układów równań metodą eliminacji Gaussa
%     Eliminacja Gausa-Jordana
% Zadania do samodzielnej realizacji
%     Zadanie 1.
%     Zadanie 2
%     Zadanie 3
% Przygotowane funkcje
%
% W wielu przypadkach zastosowań technicznych można się spotkać z problematyką rozwiązywania układów równań liniowych. Klasyczny układ równań jest postaci:
% \begin{array}{c}
% \begin{array}{c}
% a_{11} x_1 {+a}_{12} x_2 +a_{13} x_3 +\dots +a_{1n} x_n =b_1 \\
% a_{21} x_1 {+a}_{22} x_2 +a_{23} x_3 +\dots +a_{2n} x_n =b_2
% \end{array}\\
% \vdots \\
% a_{n1} x_1 {+a}_{n2} x_2 +a_{n3} x_3 +\dots +a_{nn} x_n =b_n
% \end{array}
% gdzie a_{\mathrm{ij}} jest parametrem zaś x_i jest poszukiwaną wartością – zmienną niewiadomą.
% W macierzowym ujęciu można to przedstawić następująco:
% A\times X=B
% Gdzie A – macierz z parametrami układu równań:
%     A={\left[\begin{array}{c}
% \begin{array}{c}
% a_{11} {  a}_{12}   a_{13  } \dots   a_{1n} \\
% a_{21} {   a}_{22}    a_{23} \dots a_{2n}
% \end{array}\\
% \vdots \\
% a_{n1}  a_{n2  } a_{n3 } \dots {  a}_{nn}
% \end{array}\right]}
% X – wektor zmiennych niewiadomych:
%     X={\left[\begin{array}{c}
% \begin{array}{c}
% x_1 \\
% x_2
% \end{array}\\
% \vdots \\
% x_n
% \end{array}\right]}
% Wektor B:
%     B={\left[\begin{array}{c}
% \begin{array}{c}
% b_1 \\
% b_2
% \end{array}\\
% \vdots \\
% b
% \end{array}\right]}
% Aby układ równań był jednoznacznie określony, czyli posiadał dokładnie jedno nietrywialne rozwiązanie, dla n niewiadomych konieczne jest przygotowanie n-równań niezależnych oraz przynajmniej jeden parametr b_i musi być niezerowy.
%% Przypadki szczególne układów równań
% Zanim omówimy metody rozwiązywania takich równań, zacznijmy od znacznie prostszych przypadków szczególnych. Oczywiście przykłady będziemy przedstawiać w układzie macierzowym.
%% Rozwiązywanie układów równań z elementami jedynie na głównej przekątnej
% Najprostszym przypadkiem jest układ, w którym na głównej przekątnej występują same jedynki, czyli dość trywialnie wglądający układ w postaci:
% \begin{array}{l}
% x_1 =b_1 \\
% x_2 =b_2 \\
% \vdots \\
% x_n =b_n
% \end{array}
% Dla takiego układu macierz A będzie:
%     A={\left[\begin{array}{c}
% \begin{array}{c}
% 1\;\;0  0\dots  \mathrm{ 0}\\
% 0\;\;1\mathrm{  }0\dots \;\;\;0
% \end{array}\\
% \vdots \\
% 0\mathrm{ 0}\;\;\;\;0\dots 1
% \end{array}\right]}
% Przygotowanie takiej macieży w Matlabie możemy zrealizwać używając komendy eye, poza tym musimy przygotować wektor B o takiej samej ilości wierszy co macierz A

A=eye(5)
B=[-5; 2; 6; 8; 7];

% Rozwiązaniem takiego równania są elementy wektora B, czyli:

X=B

% Możemy sprawdzić poprawność rozwiązania wyznaczając wartość iloczynu macierzowego elementów A x X. Powinien on dać dokładnie wartość elementów wektora B
% Sprawdzenie rozwiązania:

SPR=A*X
% Dla porwnania wartosc wektora B
B

%
% Sprawa się nieco tylko komplikuje, gdy na głównej przekątnej macierzy A występują elementy inne niż 1. Można je otrzymać mnożąc macierz z jedynkami na głównej przekątnej przez wektor elementów, występujących na głównej przekątnej. Załóżmy że nadal będziemy rozwiązywać układ równań z pięcioma równaniami. Wektor elementów głównej przekątnej  przygotujemy jako wektor kolumnowy, a do mnożenia wykorzystamy mnożenie element przez element czyli z operatorem ".*"

% Wektor elementow glownej przekatnej EGP
EGP=[5; 3;  1; -2; 3];
A=eye(5).*EGP

% Rozwiązaniem w takim przypadku jest wektor X otrzymany przez podzielenie elementów wektora B przez elementy stojące na gwnej przekątnej A(i,i) (wektor EPG) Możemy to uzyskać komendą diag. W celu wykonania dzielenia element wektora przez element wektora drugiego zastosujemy operator "./". Przyjmując wartość wektora B jak w poprzednim przypadku. rozwiązanie będzie wyglądać następująco:

pom=diag(A)
X=B./diag(A)

% Sprawdzenie poprawności rozwiazania
SPR=A*X
% Wektor B dla porwownania
B

%
%% Rozwiązywanie układów równań z macierzą trójkątną
% Poniżej przedstawiono przykład układu równań, gdzie wektor A stanowi macierz trójkątną dolną:
% \begin{array}{l}
% {a_{11} x}_1 =b_1 \\
% {a_{21} x}_1 +{a_{22} x}_2 =b_2 \\
% {a_{31} x}_1 +{a_{32} x}_2 +{a_{33} x}_3 =b_3 \\
% \vdots \\
% {a_{n1} x}_1 +{a_{n2} x}_2 +{a_{n3} x}_3 +\cdots +{a_{\mathrm{nn}} x}_n =b_n
% \end{array}
% Macierz A w tym przypadku posiada powyżej gwłónej przkątnej jedynie zerowe elementy. Przykładowa macierz trójkątna o wymiarze 4 na 4

A=[2  0 0 0;...
   3  2 0 0;...
   -2 2 3 0;...
   0 -5 2 1]

% Rozwiązanie układu równań z macierzą trójkątną dolną prowadzi się metodą przez podstawienie w przód czyli wyznacza się od pierwszego x1do ostatniego xn przyjmując, że:
% x_1 =\frac{b_1 }{a_{11} }
% Zaś kolejne elementy wyznacza się z zależności:
% x_k =\frac{b_k -\sum_{i=1}^{k-1} a_{\textrm{ki}} x_i }{a_{\textrm{kk}} }
% Tworząc specjalny wektor kolumnowy X zawierający same zera można ten zapis przedstawić z uwzględnieniem rachunku wektorowego:
% x_{i=1\ldotp \ldotp n} =\frac{b_i -A_{k,\ldotp \ldotp } \mathrm{xX}}{a_{\mathrm{kk}} }
% Przykład zastosowania tego podejścia w zapisie Matlaba przedstawiono w funkcji do obliczeń układów równań z macierzą trójkątną dolną przedstawioną dalej.
% Dla wektora B zdefiniowanego poniżej oraz wcześniej przyjętej macierzy A, wyznaczymy wartości wektora X stanowiącego rozwiązanie z wykorzystaniem podanej funkcji:

B=[1; 3; -1; 5]
[X]=uklad_rown_MTD(A,B)

% Podobnie jak poprzednio spawdzimy otrzymane rozwiązanie i porównamy wynik z wektorem B

SPR=A*X

% Jak ma się otrzymany wynik SPR do wektora B?
% Układ równań z macierzą trójkątną górną
% Bardzo podobnie wygląda proces rozwiązywania układu równań, dla którego wektor parametrów równania stanowi macierz trójkątna górna. Sam układ równań w tym wypadku wygląda następująco
% \begin{array}{l}
% {a_{11} x}_1 +{a_{12} x}_2 +{a_{13} x}_3 +\cdots +{a_{1\textrm{n}} x}_n =b_1 \\
% {a_{22} x}_2 +{a_{23} x}_3 +\cdots +{a_{2\textrm{n}} x}_n =b_2 \\
% \vdots \\
% {a_{\mathrm{nn}} x}_n =b_n
% \end{array}
% Dla takiego układu równań wektor parametrów równania A cechuje się tym, że poniżej głównej przekątnej występują tylko 0, vektor b jest definiowany jak poprzednio
% A=\left\lbrack \begin{array}{ccccc}
% a_{11}  & a_{12}  & a_{13}  & \cdots  & a_{1n} \\
% 0 & a_{22}  & a_{23}  & \cdots  & a_{2n} \\
%  &  &  & \vdots  & \\
% 0 & 0 & 0 & \cdots  & a_{\mathrm{nn}}
% \end{array}\right\rbrack
% Załóżmy, że tym razem mamy układ 5-u równań podany następująco
% \begin{array}{l}
% 5x_1 +2x_2 -3x_3 -x_4 +x_5 =-2\\
% 3x_1 +2x_4 -2x_5 =-1\\
% x_1 +x_2 -x_4 +2x_5 =2\\
% x_4 -3x_5 =-6\\
% \mathrm{x5}=2
% \end{array}
% Macierz parametrów jest tutaj następująca

A0=[5 2 -3 -1  1; ... 
    3 0  0  2 -2; ...
    1 1  0 -4  2; ...
    0 0  0  1 -3; ...
    0 0  0  0  1]

% Wektor B

B0=[-2; -1; 2; -6; 2]

% Czy ten układ równań stanowi macież trójkątną górną? Pewnie na pierwszy rzut oka nie, ale jak się lepiej przyjrzy i trochę przeorganizuje, to można otrzymać układ, w którym rzeczywiście macierz parametrów równania jest macierzą trójkątną górną. Po zamianie drugiego wiersza z trzecim i odpowiednim przestawieniu kolejności zmiennych x otrzymuje się układ:
% \begin{array}{l}
% -3x_3 +2x_2 +5x_1 -x_4 +x_5 =-2\\
% x_2 +x_1 -x_4 +2x_5 =2\\
% 3x_1 +2x_4 -2x_5 =-1\\
% x_4 -3x_5 =-6\\
% {\textrm{x}}_5 =2
% \end{array}
% W posaci macierzowej zapiszemy go odpowiednio:
% \left\lbrack \begin{array}{ccccc}
% -3 & 2 & 5 & -1 & 1\\
% 0 & 1 & 1 & -1 & 2\\
% 0 & 0 & 3 & 2 & -2\\
% 0 & 0 & 0 & 1 & -3\\
% 0 & 0 & 0 & 0 & 1
% \end{array}\right\rbrack *\left\lbrack \begin{array}{c}
% x_3 \\
% x_2 \\
% x_1 \\
% x_4 \\
% x_5
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% -2\\
% 2\\
% -1\\
% -6\\
% 2
% \end{array}\right\rbrack
% Zapiszmy to w Matlabie:

A=[-3 2 5 -1  1; ...
    0 1 1 -1  2; ...
    0 0 3  2 -2; ...
    0 0 0  1 -2; ...
    0 0 0  0  1]
B=[-2; 2; -1; -6; 2]

% W wektorze x musimy pamiętać, że zmienne x są zamienione, choć na razie na etapie rozwiązywania nie ma to większego znaczeni. Znaczenie to będzie miało dopiero, gdy będziemy chcieli powrócić do wyjściowej postaci układu równań.
% Rozwiązanie takiego układu prowadzi się w bardzo podobny sposób ja przedstawiono powyżej, z tym że zaczyna się od rozwiązania ostatniego równania i następnie po kolei wyznacza się kolejne niewidome w kolejności od ostatniej do pierwszej, stąd
% x_n =\frac{b_n }{a_{\mathrm{nn}} }
% Zaś kolejne elementy wyznacza się z zależności:
% x_k =\frac{b_k -\sum_{i=k+1}^n a_{\textrm{ki}} x_i }{a_{\textrm{kk}} }
% Po przyjęciu metodyki przedstawionej poprzednio funkcję do rozwiązywania takich układów równań zawarto jako uklad_rown_MTG i przedstawiono w rozdziale Przygotwane funkcje.
% Wyznaczmy wartości wektora X

[X]=uklad_rown_MTG(A,B)

% Sprawdźmy poprawność rozwiązania i porównajmy je z wektorem B:

SPR=A*X
B

% Wynik sprawdzenia jest taki sam jak wektor B, czyli uzyskaliśmy poprawne rozwiązanie.
% Możemy to także sprawdzić względem zadanego na początku układu równania. Wektor parametrów macierzy został zapamiętany jako A0 oraz odpowiadał mu wektor B0. Przed rozpoczęciem sprawdzenia musimy przywrócić wektor x do właściwej postaci. Przypomnijmy, że po reorganizacji wektor X wyglądał następująco:
%     X=\left\lbrack \begin{array}{c}
% x_3 \\
% x_2 \\
% x_1 \\
% x_4 \\
% x_5
% \end{array}\right\rbrack
% Zatem, aby powrócić do jego pierwotnego układu należy przestawić x1 z x3. Zrobimy to z dodtkową pomocniczą zmienną POM

POM=X(1);
X(1)=X(3);
X(3)=POM

% Teraz możemy sprawdzić rozwiązanie porównując iloczyn A0 razy zmodyfikowany X z B0

SPR=A0*X
B0

% Otrzymany wynik iloczynu jest zgodny z wektorem B0, co także potwierdza poprawność uzyskanego rozwiązania oraz wykonanych na wstępie przekształceń
%% Rozwiązywanie układów równań w ogólnym w ogólnym przypadku
% W praktyce inżynierskiej rzadko kiedy, albo nawet wcale nie spotykamy się ze szczególnymi rodzajami układów równań, jakie przedstawiono powyżej. Najczęściej mamy do czynienia z ogólnymi przypadkami, gdzie występuje większość parametrów stojących przy niewiadomej x. Rozwiązanie zadania w takim przypadku nie jest proste, gdy liczba niewiadomych i równań jest większa od 4.
% Dlatego wykorzystuje się metody, które pozwalają przekształcić układ równań o ogólnej postaci do jednej z prostszych omówionych powyżej. Występuje wiele metod rozwiązywania układów równań. Kilka metod zostanie przedstawione i omówione w tym opracowaniu.
%% Eliminacja Gausa
% Jedną z podstawowych metod rozwiązywania układów równań jest metoda eliminacji Gausa. Polega ona na sprowadzeniu układu równań z postaci ogólnej do układu równań z macierzą trójkątną dolną. Wykonując operacje na macierzy parametrów i równocześnie na wektorze B, krok po kroku przekształca się macierz parametrów, tworząc z niej macierz trójkątną górną. Samo rozwiązanie końcowe otrzymuje się stosując metodę rozwiązania dla macierzy trójkątnej górnej.
% Metoda eliminacji Gaussa zostanie omówiona na przykładzie przekształceń wykonanych na układzie równań:
% \left\lbrack \begin{array}{ccccc}
% 1 & 2 & 5 & -1 & 1\\
% 2 & 4 & 1 & -1 & 2\\
% 1 & 2 & 3 & 2 & -2\\
% -2 & 5 & 0 & 1 & -3\\
% 0 & 1 & 1 & 0 & 1
% \end{array}\right\rbrack *\left\lbrack \begin{array}{c}
% x_1 \\
% x_2 \\
% x_3 \\
% x_4 \\
% x_5
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% -2\\
% 2\\
% -1\\
% -6\\
% 2
% \end{array}\right\rbrack
% Stąd:

A0=[1  2  5 -1  1; ...
   2  4  1 -1  2; ...
   1  2  3  2 -2; ...
  -2  5  0  1 -3; ...
   0  1  1  0  1]
B0=[-2; 2; -1; -6; 2]
A=A0;
B=B0;

% Obliczenia będą prowadzone na zmiennych A i B, aby elemny wyjściowe pozostawić niezmienione, w celu sprawdzenia na końcu poprawności wyniku.
% Przygoptujmy wektor kolumnowy mnożnika dla pierwszej kolumny C, przypisując mu 0 jako pierwszy element oraz elementy od 2 do 5 z pierwszej kolumny macierzy A

C=[0;A(2:5,1)]

% Modyfikujemy ten wektor dzieląc go przez A(1,1)

C1=C./A(1,1);

% Wyznaczamy odjemnik M dla macierzy parametrów A mnożąc wektor C1 przez pierwszy wiersz macierzy A

OA=C1*A(1,:)

% Odejmujemy od macierzy A macierz M, wyznaczając w ten sposób nową wartość macierzy A w której w pierwszej kolumnie poza pierwszym niezerowym parametrem pozostałe są zerami:

A=A-OA

%  Podobne operacje wykonuje się na wektorze B. Najpierw wznacza się odjemnik OB mnożąc wektor C1 przez pierwszy element wektora B, anastępnie od wektora B odejmuje się wektor OB

OB=C1*B(1)
B=B-OB

% W ten sposób został zakończony pierwszy krok obliczeń. W następnym kroku eliminacji w kolumnie drugiej wyzerowane zostaną elementy poniżej głównej przekątnej. Jednakże ze względu na to, że na głównej przekątnej występuje 0 (macierz A pozycja 2,2), zatem nie jesteśmy w stanie wyznaczyć wektora C. W celu wyeliminowania tego problemu dokonuje się przestawienia wierszy macierzy poniżej. Wyszukuje się w drugiej kolumnie poniżej drugiego wiersza elementu największego, bądź równego 1. W przedstawionym przypadku w 5-tym wierszu występuje 1, zatem wiersz piąty zamienimy z wierszem 2. Równocześnie z tą operacją trzeba w wektorze B zamienić element 5 z drugim. Uwaga – Gdyby była dokonywana zamiana kolumn, to wtedy zamienia się pozycje odpowiednich elementów w wektorze X (operacja mniej zalecana). Po zaproponowanej zmianie macierz A i wektor B będą miały postać jak poniżej, zaś elementy wektora X się nie zmienią:

% Wprowadzamy zmienną pomocniczą POM ktrej przypisuje się wartosc 2-go wiersza 
% macierzy A
POM=A(2,:);
% Podstawiamy najpierw piąty wiersz w miejsce pierwszego, a następnie zmienną POM
% za piąty wiersz
A(2,:)=A(5,:);
A(5,:)=POM
%Podobnie postąpimy z wektorem B:
POM=B(2);
B(2)=B(5);
B(5)=POM

% Teraz można przystąpić do wyznaczenia mnożnika C, który będzie zawierał elementy 2-giej kolumny wektora A poniżej głównej przekątnej, a od góry do drugiego elementu będzie zawierał 0.

C=[0;0;A(3:5,2)]

% Powtarzamy następne operacje. Tym razem modyfikujemy wektor dzieląc go przez A(2,2)

C1=C./A(2,2);

% Wyznaczamy odjemnik M od macierzy parametrów A mnożąc wektor C1 przez tym razem drugi wiersz macierzy A

OA=C1*A(2,:)

% Odejmujemy od macierzy A macierz M, wyznaczając w ten sposób nową wartość macierzy A w której w pierwszej kolumnie poza pierwszym niezerowym parametrem pozostałe są zerami:

A=A-OA

%  Podobne operacje wykonuje się na wektorze B. Najpierw wznacza się odjemnik OB mnożąc wektor C1 przez drugi element wektora B, a następnie od B odejmuje się OB

OB=C1*B(2)
B=B-OB

% W ten sposób został zakończony drugi krok obliczeń. Poniewarz tym razem na trzeciej pozycji głównej przekątnej macierzy A nie występuje 0, dlatego możemy kontynuwać eliminację bez zamiany wierszy. Powtarzamy operację jak wyżej, tym razem działania mają doprowadzić do wyzerowania elementów poniżej głównej przekątnej w trzeciej kolumnie macierzy A.

C=[0;0;0;A(4:5,3)]

C1=C./A(3,3);
OA=C1*A(3,:)
A=A-OA
OB=C1*B(3)
B=B-OB

% Ostatni krok operacji jest realizowany dla przedostatniej kolumny. Ponieważ w tym wypadku też nie trzeba przestawiać elementów, bo A(4,4) jest różne od 0. Więc jak wcześniej liczymy:

C=[0;0;0;0;A(5,4)]
C1=C./A(4,4)
OA=C1*A(4,:)
A=A-OA

OB=C1*B(4)
B=B-OB

% Otrzymaliśmy zmodyfikowaną w wyniku eliminacji Gausa macierz A i odpowiednio do tego przerobiony wektor B. Otrzymano to po wykonaniu n-1 kroków (n- jest ilością rwównań). Na razie nie mamy jeszcze rozwiązania, ale otrzymana macierz trójkątna górna, nie stanowi już dla nas problemu. Rozwiążemy ją metodą przez podstawienie w tył. Do tego wykorzystamy odpowiednią funkcję uklad_rown_MTG przedstawioną poniżej. Rozwiązanie stanowią wartości wektora X

[X]=uklad_rown_MTG(A,B)

% Sprawdźmy podobnie jak poprzednio rozwiązanie. Zrobimy to mnożąc macierz wyjściową A0 przez wektor X, czego wynik porównamy z wartością wektora wyjściowego B0

SPR=A0*X
B0

% Wyznaczone wartości wektora SPR są takie same jak wartość B0. Stąd wynika, że uzyskany wynik rowiązania X jest poprawny.
%% Funkcja do rozwiązywania układów równań metodą eliminacji Gaussa
% Opracowując funkcję do rozwiązywania układów równań liniowych z wykorzystaniem eliminacji Gaussa, zaczniemy od opracowania funkcji realizującej eliminację Gaussa.
% Funkcję realizującą eliminację Gaussa przedstawiono w funkcjach poniżej jako eliminacja_Gaussa. Funkcję podzielono na dwie części. Piersza część odpowiada za wyeliminowanie 0 z głównej przekątnej macierzy A. Przy czym jak można zauważyć nie ma w funkcji żadnej komendy wskazującej na sprawdzanie wartości elementu na głównej przekątnej. W funkcji jest zawarta z automatu opcja, że za każdym razem następuje przeszukiwanie i-tej kolumny poniżej głównej przekątnej i podmiana wierszy tak aby na przekątnej był element największy co do wartości bezwzględnej z tej kolumny z elementów poniżej. Takie rozwiązanie eliminuje błędy jakie mogą się pojawić podczas dzielenia przez liczby bliskie zeru. Druga część przedstawionego kodu zawiera proce obliczeń odpowiadający właściwej eliminacji Gagussa.
% W wyniku wykonania eliminacji Gaussa uzyskuje się macierz trójkątną górną i odpowiednio zmodyfikowany wektor B. Dopiero w kolejnym roku następuje rozwiązanie przez rozwiązanie macierzy trójkątnej górnej.
% Przyjrzyjmy się jak to działa. Weżmiemy układ równań jak na wstępie czyli A0*X=B0, gdzie:

A0
B0

% Po podstawienia tej macierzy i wektora do funkcji realizującej eliminację Gausa otrzymamy macierz trójkątną górną i zmodyfikowany wektor:

[A,B]=eliminacja_Gaussa(A0,B0)

% Teraz wykorzystując funkcję do rozwiązania macierzy trójkątnej górnej uzyskujemy ostateczne rozwiązanie

X=uklad_rown_MTG(A,B)

% A zatem funkcja do rozwiązywania układów równań będzie złożeniem dwóch opisanych wyżej funkcji. Zawarto to w funkcji  rozw_ukl_rown_EG. Działanie funkcji daje od razu wynik rozwiązania:

X=rozw_ukl_rown_EG(A0,B0)

%
%% Eliminacja Gausa-Jordana
% Modyfikacją metody Gaussa jest metoda Gaussa-Jordana. Polega ona na tym, że w wyniku realizacji przekształceń z macierzy ogólnej parametrów otrzymuje się macierz jednostkową (jedynki na przekątnej głównej)
% Kolejne kroki w realizacji tej metody polegają na tym, że na początku dzieli się i-ty (1-szy, w kolejnym kroku 2-gi itd.) wiersz macierzy A i wektora B przez element na głównej przekątnej macierzy A w i-tym wierszu. W ten sposób na głwónej przekątnej otrzymuje się jedynkę. UWAGA najpierw należy dzielić wektor B, a dopiero później macierz A. Gdy nie zachowa się tej kolejności, to wektor B będzie podzielony przez nową zmodyfikowaną wartość A(i,i) inną od tej przez którą dzieli się macierz A.

A=A0
B=B0
B(1)=B(1)/A(1,1)
A(1,:)=A(1,:)/A(1,1)

% Oczywiście w pokazanym przykładzie, wykonana operacja niczego nie zmieniła, bo element A(1,1)=1
% Następnie wyznacza się wektor mnożnika C tak, że w i-tym (w tym wypadku pierwszym) wierszu jest 0, a w pozostałych wierszach są elementy z i-tej (pierwszej) kolumny macierzy A

C=[0; A(2:5,1)]

% Wyznacza się odjemnik dla macierzy A mnożac wektor C przez i-ty (pierwszy) wiersz macierzy A. Po odjęciu wyznaczonego odjemnika od A, macirz A powinna zawierać 1 na pozycji A(1,1) i 0 w pozstałych pozycjach pierszyego wiersza.

OA=C*A(1,:)
A=A-OA

% Podobnie pstępuje się dla wektora B:

OB=C*B(1)
B=B-OB

% Po tej operacji podobnie jak poprzednio otrzymano, że w macierzy A na pozycji 2,2 występuje 0. Aby kontynuować obliczenia, trzeba zamienić wiersze. Jeżeli w tym wypadku zamienimy wiersz 2-gi z 5-tym wierszem, to na głównej przekątnej w pozycji A(2,2) pojawi się 1, co spowoduje, że nie będzie trzeba dzielić przez A(2,2). Od razu można będzie przejść do wyznaczania mnożnika, który w tym przypadku będzie zawierał elementy w 2-ej kolumny macierzy A, i tylko na drugiej pozycji będzie miał wartość 0.

% Zamiana elementow 2 i 5 macierzy A
POM=A(2,:);
% Podstawiamy najpierw piąty wiersz w miejsce pierwszego, a następnie zmienną POM
% za piąty wiersz
A(2,:)=A(5,:);
A(5,:)=POM
%Zamiana elementw 2 i 5 wektora B:
POM=B(2);
B(2)=B(5);
B(5)=POM
%Wyznaczenie mnożnika dla 2-giego wiersza
C=A(:,2);
C(2)=0

% Teraz wyznaczymy odjemnik dla macierzy A i wektora B

%Odjemnik macierzy A
OA=C*A(2,:)
%Odjemnik wektora B
OB=C*B(2)

% Po odjęciu OA od A w drugiej kolumnie zostanie 1 na głównej przekątnej i reszta będą 0.

A=A-OA
B=B-OB

% Przyglądając się macierzy A, można zauważyć, że zamieniając wiersz trzeci z czwartym dostalibyśmy na głównej przekątnej 1, co wyeliminowało by konieczność dzielenia przez A(3,3). Ponieważ nie będziemy tego robić, dlatego tym razem znowu podzielmy 3 wiersz macierz macierzy A i wektora B przez A(3,3)

B(3)=B(3)/A(3,3)
A(3,:)=A(3,:)/A(3,3)

% Znowu wyznaczamy mnożnik i odjemniki, które wylicza się odpowiednio

C=A(:,3);
C(3)=0;
OA=C*A(3,:)
OB=C*B(3)

% Po odjęciu A i B:

A=A-OA
B=B-OB

% Do wykonania zostały jeszcze dwa kroki. Pierwszy modyfikujący wiersz 4 i kolejny modyfikujący wiersz 5. W metodzie Gaussa-Jordana należy zatem wykonać o jeden krok więcej niż w klasycznej eliminacji Gausa. Krok czwarty:

B(4)=B(4)/A(4,4)
A(4,:)=A(4,:)/A(4,4)
C=A(:,4);
C(4)=0;
OA=C*A(4,:)
OB=C*B(4)
A=A-OA
B=B-OB

% Krok 5

B(5)=B(5)/A(5,5)
A(5,:)=A(5,:)/A(5,5)
C=A(:,5);
C(5)=0;
OA=C*A(5,:)
OB=C*B(5)
A=A-OA
B=B-OB

% Wyznaczony wektor B jest w tym wypadku rozwiązaniem równania, czyli poszukiwaną wartością X.
% Funkcję eliminację Gaussa-Jordana przedstawiono w wykazie funkcji eliminacja_G_J. Pierwsza część w przedstawionej funkcji dotyczy zamiany elementów, jak w eliminacji Gaussa, natomiast druga część zawiera właściwy kod eliminacji Gausa-Jordana. Następnie funkcja ta została wykorzystana w funkcji do rozwiązywania układów równań w funkcji ozw_ukl_rown_EG

[X]=rozw_ukl_rown_EG_J(A0,B0)

% Sprawdźmy otrzymane rozwiązanie:

SPR=A0*X
B0

% Widzimy, że wektor SPR jest rwny B0
%
%% Zadania do samodzielnej realizacji
%% Zadanie 1.
% Dla przedstawionego układu równań zapisz odpowiednią macierze i wektor a następnie wykorzystaj jedną z przedstawionych funkcji do rozwiązania układu. Sprawdź czy uzyskany wynik jest poprawny. Jaką najprostszą metodę można zastosować do rozwiązania tego układu?
% \begin{array}{l}
% -2x_1 +3x_2 +5x_3 -x_4 +x_5 =2\\
% x_1 +{2x}_2 -x_3 +2x_4 =3\\
% 3x_1 +x_2 -2x_3 =-1\\
% x_1 -3x_2 =-6\\
% {\textrm{x}}_1 =2
% \end{array}






%
%% Zadanie 2
% Dla podanego poniżej układu równań zapisz macierz A0 i wektor B0, a następnie przeorganizuj układ tak aby otrzymać macierz trójkątną górną. Dla tak zmodyfikowanego układu przygotuj macierz A i wektor B, a następnie rozwiąż układ równań i sprawdź wynik względem początkowego układu równań. Zastanów się, czy wykonane przestawienia wpłynęły na zmianę układu wektora X?
% \begin{array}{l}
% -2x_1 +3x_2 +5x_3 -x_4 +x_5 =2\\
% x_1 +2x_2 -{3x_5 +x}_6 =-1\\
% x_1 +{2x}_2 -x_3 +2x_4 =3\\
% 3x_1 +x_2 -2x_3 =-1\\
% x_1 -3x_2 =-6\\
% {\mathrm{x}}_1 =2
% \end{array}
%





%
%% Zadanie 3
% Rozwiąż poniższy układ równań, Sprawdź poprawność rozwiązania.
% \begin{array}{l}
% -3x_1 +2x_2 +5x_3 -x_4 +x_5 =-2\\
% x_1 +x_2 -x_4 +2x_5 =2\\
% 3x_1 -2x_2 +x_3 +2x_4 -2x_5 =-1\\
% {x_1 -2x_2 +x_3 +x}_4 -3x_5 =-6\\
% {{-x}_1 +2x_2 -x_3 -4x}_5 =2
% \end{array}
%




%
%% Przygotowane funkcje
% Funkcja do rozwiązywania układów równań z macierzą trójkątną dolną

function [X]=uklad_rown_MTD(U,B)
    [i,~]=size(U);
    
    X=zeros(i,1);
    
    for k=1:i
        c=U(k,:)*X;
        X(k)=(B(k)-c)/U(k,k);
    end
end


% Funkcja do rozwiązywania układów równań z macierzą trójkątną górną

function [X]=uklad_rown_MTG(U,B)

    [i,~]=size(U);
    
    X=zeros(i,1);
    
    for k=i:-1:1
        c=U(k,:)*X;
        X(k)=(B(k)-c)/U(k,k);
    end
end

% Funkcja realizująca eliminację Gaussa

function [a1,b1]=eliminacja_Gaussa(a,b)
n=length(b);
    for i=2:n
        % Częśc funkcji realizująca zamianę elementu z gwnej przekątnej w i-tej kolumnie
        % na element o najwyżsej wartości bezwzględnej
        [~,im]=max(abs(a(i-1:n,i-1)));
        im=im+i-2;
        if im~=i-1
            am=a(im,:);
            a(im,:)=a(i-1,:);
            a(i-1,:)=am;
            am=b(im);
            b(im)=b(i-1);
            b(i-1)=am;
        end
            
        % Waściwa częśc eliminacji Gaussa  
        c(1:i-1,1)=0;
        c(i:n,1)=a(i:n,i-1);
        c=c/a(i-1,i-1);
        a=a-c*a(i-1,:);
        b=b-c*b(i-1,1);
    end
    a1=a;
    b1=b;
 end

% Funkcję do rozwiązywania układów równań liniowych z wykorzystaniem eliminacji Gaussa

function x=rozw_ukl_rown_EG(a,b)
% a - pinnna byc macierzą kwadratową parametrow rownania
% b - jest wektorem zawierającym elementy drugiej strony rownania
% a*x=b;

    [x1,~]=size(b);
    if x1==1
        b=b';
    end
    [a,b]=eliminacja_Gaussa(a,b);
    [x]=uklad_rown_MTG(a,b);
end
    

% Funkcja realizująca eliminację Gaussa-Jordana

function [a1,b1]=eliminacja_G_J(a,b)
n=length(b);
    for i=1:n
        % Częśc funkcji realizująca zamianę elementu z gwnej przekątnej w i-tej kolumnie
        % na element o najwyżsej wartości bezwzględnej
        [~,im]=max(abs(a(i:n,i)));
        im=im+i-1;
        if im~=i
            am=a(im,:);
            a(im,:)=a(i,:);
            a(i,:)=am;
            am=b(im,:);
            b(im,:)=b(i,:);
            b(i,:)=am;
        end
            
        % Waściwa częśc eliminacji Gaussa-Jordana 
        b(i,:)=b(i,:)/a(i,i);
        a(i,:)=a(i,:)/a(i,i);
       
        c=a(:,i);
        c(i,1)=0;
        a=a-c*a(i,:);
        b=b-c*b(i,:);
    end
    a1=a;
    b1=b;
end


 

% Funkcję do rozwiązywania układów równań liniowych z wykorzystaniem eliminacji Gaussa Jordana

function x=rozw_ukl_rown_EG_J(a,b)
% a - pinnna byc macierzą kwadratową parametrow rownania
% b - jest wektorem zawierającym elementy drugiej strony rownania
% a*x=b;

    [x1,~]=size(b);
    if x1==1
        b=b';
    end
    [a,x]=eliminacja_G_J(a,b);
end
    