%% Rozwiązywanie układów równań liniowych
%
% Zapisujemy układ równań:   x+3y-4z=2,\;\;\;2x+4y+z=-2,\;\;\;\;x-y+z=0\ldotp \;

syms x y z 
uklad=[x+3*y-4*z==2,2*x+4*y+z==-2, x-y+z==0]

%
%% Wizualizacja:
% Gdy mamy równania z trzema niewiadomymi, możemy narysować płaszczyzny odpowiadające poszczególnym równaniom.
%

p1=fimplicit3(uklad(1), [-10 10],'FaceColor','g'); % zielona płaszczyzna - wykres zbioru rozwiązań I równania 
% dla x,y,z z przedzialu [-10,10]
p1.FaceAlpha=0.3;   % ustalamy poziom przezroczystości 

%

hold on
p2=fimplicit3(uklad(2), [-10 10],'FaceColor','b'); % niebieska płaszczyzna - dla rownania II
p2.FaceAlpha=0.4; 

% Przecięcie dwóch płaszczyzn nierównoległych jest prostą.
%


p3=fimplicit3(uklad(3), [-10 10],'FaceColor','y'); % żółta płaszczyzna dla równiania III
p3.FaceAlpha=0.4;
hold off

% W typowej sytuacji trzy płaszczyzny przecinają się w jednym punkcie.
%
%% Przykład 1.
%% Rozważmy układ oznaczony (3 równania i 3 niewiadome)
%
% Zapisujemy macierz współczynników i wektor wyrazów wolnych związanych z danym układem równań

[A,b]=equationsToMatrix(uklad)    % A - macierz współczynników, b - wektor wyrazów wolnych

%
% Sprawdźmy, czy istnieje rozwiązanie. W tym celu wyznaczymy rząd macierzy A i macierzy [A|b].

n=3  % zapisujemy liczbę niewiadomych
rA=rank(A)
rAb=rank([A,b])

%
% Na podstawie twierdzenia Kroneckera-Capellego wnioskujemy, że układ posiada dokładnie jedno rozwiązanie.
% Rozwiązanie układu Cramera otrzymamy wpisując polecenie linsolve(A,b)  lub A\b.

linsolve(A,b) 
A\b

%
%% Przykład 2.
%% Rozważmy układ z dwoma równaniami o trzech niewiadomych

rownania1i2=[uklad(1),uklad(2)] % wybieramy dwa równania z początkowego układu
[A1,b1]=equationsToMatrix(rownania1i2)
macierzAB=[A1,b1] % macierz rozszerzona układu

%
% Sprawdźmy, czy istnieje rozwiązanie. W tym celu wyznaczymy rząd macierzy A1 i macierzy [A1|b1].

n=3
rA1=rank(A1)
rA1b1=rank(macierzAB)

%
% Na podstawie twierdzenia Kroneckera-Capellego wnioskujemy, że układ posiada nieskończenie wiele rozwiązań zależnych od 1 parametru.
% Użyjemy polecenia linsolve.

X0=linsolve(A1,b1)

% Układ ten posiada nieskończenie wiele rozwiązań zależnych od 1 parametru, ale  Matlab podaje tylko jedno z rozwiązań (rozwiązanie szczególne).
%
% Za pomocą polecenia rref(macierz rozszerzona ukladu)
% wyznaczymy postać zredukowaną macierzy rozszerzonej układu.
% Jest to postać macierzy, którą można otrzymać wykonując operacje na wierszach.
% Zawiera maksymalną liczbę niezależnych kolumn jednostkowych.

nowaA1b1=rref([A1 b1])

%
% Z tej postaci możemy już odczytać rozwiązanie.
% Zapisujemy równania

nowerownania = nowaA1b1*[x;y;z;0] == nowaA1b1*[0;0;0;1]

%
% Teraz zapiszemy rozwiązanie ogólne układu niejednorodnego.
% Dostaniemy je przenosząc w równaniach  wyrażenia ze zmienną z na prawą stronę równań.

rozw = nowaA1b1*[x;y;0;0] == nowaA1b1*[0;0;-z;1]

% Zmienne x i y zależą od z, gdzie \mathit{\mathbf{z}}\in R\ldotp
%
%% Przykład 3.
%% Rozważmy układ 3 równań z czterema niewiadomymi

syms x y z u
eq=[x+3*y-4*z-u==2,2*x+4*y+z-u==-2, x+y+5*z==-4] % zapisujemy układ
[Aa,Bb]=equationsToMatrix(eq)
[AA,bb]=equationsToMatrix(eq,[x y z u]) % zapisujemy macierze związane z układem (określamy kolejność zmiennych)
% i wyznaczamy ich rzędy
rAA=rank(AA)
rAAbb=rank([AA,bb])

%
% Na podstawie twierdzenia Kroneckera-Capellego wnioskujemy, że układ posiada nieskończenie wiele rozwiązań zależnych od 2 parametrów.
% Znajdźmy rozwiązanie szczególne.
%

X0=AA\bb

%
% Dalej wykorzystamy fakt, że każde rozwiązanie może być zapisane jako suma rozwiązania szczególnego i kombinacji liniowej wektorów z układu  fundamentalnego rozwiązań układu jednorodnego.
% Polecenie null(AA) zwraca macierz Ker(AA) - jądro macierzy AA, która w kolumnach ma zapisane  wektory ukladu fundamentalnego dla równania AA*X=0.

fund=null(AA)
X1=fund(:,1) % zapisujemy pierwszy wektor
X2=fund(:,2) % zapisujemy drugi wektor

%
% Ostatecznie wektor rozwiązań zapiszemy jako XX, parametry to a i b.
% Będzie to rozwiązanie ogólne układu niejednorodnego.

syms a b % rozwiazanie zapiszemy z wykorzystaniem parametrow a i b
XX=X0+a*X1+b*X2 % kolejne wyrazy tego wektora odpowiadają zmiennym x, y, z, u.

%
% Sprawdźmy poprawność rozwiązania

isequal(AA*XX,bb)

%
%% Przykład 4.

 clear;
close all;

%
%% Zbadamy rozwiązalność układu w zależności od parametru a.
% Zajmiemy się układem trzech równań z czterema niewiadomymi.

syms x y z t a
eq=[6*x+2*y+(a+4)*z-2*t==a,3*x+y+3*z-t==a-1, 3*x+a*y+2*z-t==a-1] % zapisujemy układ
[A,B]=equationsToMatrix(eq,[x y z t]) % w nawiasie kwadratowym zapisano, które zmienne są niewiadomymi (pominięto a)
AB=[A,B] % macierz rozszerzona układu

%
% Standardowe polecenia np. linsolve nie dają pełnej informacji o rozwiązaniach układu.
% Liczba niewiadomych n=4, a liczba równań to 3, więc nie dostaniemy jednoznacznego rozwiązania.
% Wykorzystamy twierdzenie Kroneckera Capellego.
% Uwaga: Polecenie rank(A) nie zwróci informacji, jaki będzie rząd macierzy dla różnych a, poda  maksymalny możliwy rząd.
% Wyznaczymy rząd macierzy, badając minor macierzy A maksymalnego stopnia (3).
% Wybieram podmacierz macierzy A zbudowaną z jej trzech kolumn i obliczę jej wyznacznik.

mA=A(:,[1,2,3]) % wybieram podmacierz macierzy A skladajaca sie z kolumn 1,2,3.
minor=det(mA)
solve(minor==0,a) % znajdujemy pierwiastki otrzymanego wielomianu

%
%% Wybrany minor jest niezerowy, gdy a\not= 1 i gdy a\not= 2. Wtedy r(A)=r(A|B)=3, więc istnieje nieskończenie wiele rozwiązań zależnych od 1 parametru.
%% Osobno wyznaczymy rzędy macierzy, gdy a=1 lub a=2.
%
%% Badamy pierwszy wyjątek: a=1

A1=subs(A,a,1) % Podstawiam a=1 do macierzy A (polecenie subs w obiekcie A wstawia a=1)
AB1=subs(AB,a,1) % Podstawiam a=1 do macierzy AB
B1=AB1(:,5) % zapisuję kolumnę wyrazów wolnych
rank(A1) % sprawdzam rząd macierzy A dla a=1
rank(AB1) % sprawdzam rząd macierzy AB dla a=1

% W przypadku a=1 rzędy macierzy są różne, więc układ nie ma rozwiązań.
%
% Sprawdźmy jaką informację poda Matlab w przypadku układu sprzecznego.

linsolve(A1,B1) % w przypadku macierzy symbolicznych, poszukiwane jest rozwiązanie dokładne

% Inaczej będzie w przypadku macierzy numerycznych. Wtedy szukane będzie rozwiązanie przybliżone.

num_A1=eval(A1) % tworzona jest numeryczna wersja macierzy A1
num_B1=eval(B1) % tworzona jest numeryczna wersja macierzy B1
num_x=linsolve(num_A1,num_B1) % otrzymamy rozwiązanie przybliżone

%
%
%% Badamy drugi wyjątek.
% Powtarzamy procedurę badania rzędów, przyjmując a=2.

A2=subs(A,a,2)
AB2=subs(AB,a,2)
B2=AB2(:,5)
rank(A2)
rank(AB2)

% W przypadku a=2 rzędy macierzy są jednakowe, równe 2, więc układ ma nieskończenie wiele rozwiązań zależnych od 2 parametrów.
%
% Rozwiążemy ten układ metodą eliminacji wykonując operacje na wierszach macierzy rozszerzonej.

AB2 % macierz rozszerzona układu

%

AB2(1,:)=AB2(1,:)-2*AB2(2,:)

%
% Pomijamy dalej wiersz 1. o wyrazach zerowych

AB2=AB2([2,3],:)
% pozostawiamy jedynie wiersz 2 i 3

%
% Wykonuję operacje na wierszach, takie by uzyskać w macierzy A dwie kolumny jednostkowe

AB2(1,:)=AB2(1,:)-AB2(2,:) 

%

AB2(2,:)=AB2(2,:)-2*AB2(1,:) 

%

AB2(2,:)=-AB2(2,:) 

%
% Kolumny jednostkowe to kolumny 3 i 4. Odpowiadają im zmienne z i t.
% Wyznaczymy te zmienne, traktując x i y jako parametry.
% Zapisuję otrzymany układ równań..

rownania=AB2*[x; y; z; t; 0]==AB2*[0; 0; 0; 0; 1]

%
% Teraz zapiszemy rozwiązanie ogólne tego układu.
% Dostaniemy je przenosząc w równaniach  wyrażenia ze zmienną x i y na prawą stronę równań.

rozw = AB2*[0;0;z;t;0] == AB2*[-x;-y;0;0;1]
XX=[x;y;AB2*[-x;-y;0;0;1]]  % wektor rozwiązań
isequal(A2*XX,B2) % sprawdzenie poprawności rozwiązania