Rozdział 1.Wiadomości wstępne
Prawidłowa eksploatacja istniejących systemów gazowniczych, jak również ich rozbudowa nie może być realizowana bez pomocy programów obliczeniowych. Prawidłowa eksploatacja sieci to przede wszystkim realizacja dostaw gazu odbiorcom zgodnie z podpisanymi umowami, przy jednoczesnej minimalizacji szeroko rozumianych kosztów eksploatacji. Z kolei prawidłowa rozbudowa sieci to wybór wariantu, który spełni warunki przyjętego kryterium kosztów eksploatacji oraz budowy.
Oprogramowanie wykorzystujące modele matematyczne elementów sieci powinno zatem stanowić podstawowe narzędzie pracy służb odpowiedzialnych za transport gazu.
Model matematyczny jest układem równań, który w określony sposób i z określoną dokładnością reprezentuje rzeczywisty system. Model matematyczny jest oczywiście uproszczeniem systemu rzeczywistego, ponieważ zawiera tylko te aspekty systemu, które odpowiadają wymaganiom stawianym modelowi.
Modele można klasyfikować zgodnie z założeniami przyjętymi przy ich tworzeniu. Jeżeli system ma zmienne, które zmieniają się w sposób ciągły w czasie, to należy używać modeli dynamicznych opisywanych równaniami ogólnie różniczkowymi. Jeżeli zmienne nie są funkcją czasu, wtedy należy stosować model statyczny opisywany równaniami algebraicznymi. Wybór modelu zależy od konkretnej sytuacji. Na przykład, w sieci gazowej niskociśnieniowej zmiany ciśnienia i przepływu są bardzo szybkie. Są one pomijane w większości rozpatrywanych przypadków, stąd stosowanie modeli stanu ustalonego - równań algebraicznych. W przypadku sieci wysokiego ciśnienia dynamika przepływającego gazu jest znacznie wolniejsza ze względu na ogromne ilości gazu zmagazynowane w rurociągach. Pominięcie jej prowadziłoby do dużych błędów w opisie zjawiska. Stąd konieczność stosowania modeli matematycznych w postaci równań różniczkowych.
Głównym celem stosowania programów symulacyjnych jest badanie zachowania się sieci gazowych w określonych warunkach. Rozwiązanie określonego układu równań algebraicznych lub różniczkowych opisujących sieć dla założonych parametrów obciążenia i zasilania sieci da nam odpowiedź na pytania dotyczące wartości ciśnienia lub przepływu w wybranych punktach lub fragmentach sieci. Dzięki algorytmom symulacji możemy znacznie zmniejszyć nakłady na wyposażenie pomiarowe systemów gazowniczych, a także, zakładając poprawność modeli, uzyskiwać pełną i dokładną informację o systemie. Efektywne zarządzanie takim systemem wymaga zastosowania zaawansowanych metod numerycznych, algorytmów obliczeniowych o dużym stopniu złożoności i odpowiedniej jakości środków technicznych.
W zależności od wykorzystywanych modeli matematycznych opisujących przepływ gazu w gazociągu symulację dzielimy na statyczną i dynamiczną.
Symulacja statyczna sieci gazowej [4]
Algorytmy symulacji statycznej opracowywane są w oparciu o:
- I prawo Kirchhoffa
[1.1]
gdzie:
- wektor obciążenia w węzłach,
- wektor przepływu w łukach grafu sieci,
- zredukowana macierz incydencji węzłów i łuków,
n - liczba węzłów,
m - liczba łuków,
n1 - liczba węzłów wyróżnionych (źródła);
- II prawo Kirchhoffa
[1.2]
gdzie:
- wektor spadków ciśnienia w łukach,
- macierz incydencji oczek podstawowych i łuków;
k - liczba niezależnych oczek;
- równanie przepływu, które jest nieliniowym równaniem algebraicznym Q = f (?p), a dla całej sieci jest postaci
Q = ?(?P)
[1.3]
gdzie:
?(?P) - wektor funkcji spadku ciśnienia w łukach,
- wektor spadku ciśnienia w łukach.
Danymi wejściowymi do programów symulacji statycznej są:
- struktura topologiczna sieci,
- wymiary geometryczne poszczególnych rurociągów,
- jakość wewnętrznej powierzchni rurociągów,
- wartości obciążenia oraz parametry zasilania w określonych węzłach sieci.
Rezultatem symulacji statycznej są wartości przepływu w łukach spełniające I prawo Kirchhoffa, wartości ciśnienia w węzłach zapewniające spadki ciśnienia w każdym łuku spełniające II prawo Kirchhoffa w każdym oczku, a relacje między spadkiem ciśnienia a przepływem w każdym łuku spełniają przyjęte równanie przepływu.
Symulacja dynamiczna sieci gazowej [3],[5]
Stany nieustalone w sieci występują w sytuacji jej niezbilansowania. Ilości gazu dostarczane do sieci w określonym przedziale czasu są różne od ilości gazu odbieranych z sieci. Mamy wtedy do czynienia ze stanem sczerpywania sieci - gdy ilości dostarczane są mniejsze od odbieranych, lub ze stanem ładowania w przypadku odwrotnym. Wartość ciśnienia w dowolnym punkcie sieci jest funkcją drogi x (przebytej przez strumień gazu) oraz czasu t. Nieustalony przepływ gazu może być opisany układem równań różniczkowych zwyczajnych lub cząstkowych typu hiperbolicznego lub parabolicznego. Model matematyczny nieustalonego przepływu gazu wyprowadzany jest przy założeniu pewnych uproszczeń, z których najważniejszym jest założenie jednowymiarowości przepływu. Przyjęcie przepływu gazu jako jednowymiarowy jest zasadne, jeśli długość przewodu jest wielokrotnie dłuższa od jego średnicy, szybkość zmian takich parametrów jak ciśnienie, gęstość, temperatura czy prędkość w kierunku normalnym do kierunku linii prądu jest pomijalna w stosunku do szybkości zmian tych parametrów wzdłuż linii prądu. Oznacza to, że parametry przepływu gazu w dowolnym przekroju poprzecznym rurociągu traktowane są jako stałe, a przepływ jako jednorodny. Poza tym przyjmuje się, że promień krzywizny gazociągu jest duży w porównaniu z jego średnicą - przekrój poprzeczny gazociągu jest stały (przedziałami stały), rurociąg jest sztywny (nie odkształca się pod wpływem wewnętrznego ciśnienia panującego w rurociągu) oraz kształt profili prędkości i temperatury (w przypadku przepływu nieizotermicznego) są w przybliżeniu stałe wzdłuż gazociągu.
Model matematyczny nieustalonego przepływu gazu
Do tworzenia modelu matematycznego nieustalonego przepływu gazu w rurociągu stosuje się następujące równania [3]:
- równanie ciągłości,
- równanie ruchu,
- równanie energii,
- równanie stanu.
Równanie ciągłości
Równanie ciągłości jest zasadą zachowania materii w odniesieniu do zjawiska przepływu płynu. Masa, która jest miarą ilości materii, nie może powstać ani zniknąć, stąd zmiana masy w objętości kontrolnej może być spowodowana tylko różnicą wpływu i wypływu do i z objętości. Równanie (1.4) jest prawem zachowania masy w postaci różniczkowej
[1.4]
gdzie: ? - gęstość gazu (kg/m3), w - prędkość strumienia gazu (m/s2).
Mnożąc obie strony równania przez pole przekroju poprzecznego A oraz uwzględniając, że M = ? - w - A, gdzie M - przepływ masowy, otrzymamy
[1.5]
Zakładając izotermiczny proces przepływu gazu, równanie stanu możemy zapisać
? c2 = p
[1.6]
gdzie:
c - izotermiczna prędkość dźwięku w gazie,
p - ciśnienie gazu.
Podstawiając równanie (1.4) do (1.6), otrzymujemy
[1.7]
Równanie ruchu
Równanie ruchu (II prawo dynamiki Newtona) mówi, że algebraiczna suma sił działających na element płynu o określonej masie w określonej chwili jest równa zmianie pędu tego elementu w tej chwili
[1.8]
Lewa strona równania (1.8) może być wyrażona następująco:
[1.9]
W przypadku transportu gazu rurociągiem wyróżnia się trzy składowe siły działającej na element gazu
[1.10]
gdzie:
F1 - siła od ciężaru gazu znajdującego się wewnątrz powierzchni kontrolnej,
F2 - siła od oporu hydraulicznego,
F3 - siła od ciśnienia gazu,
? - kąt nachylenia osi rurociągu do poziomu,
? - współczynnik oporu hydraulicznego Darcy'ego-Weisbacha,
D - średnica rurociągu.
Po przekształceniach otrzymujemy
[1.11]
Równanie (1.11) jest ogólną postacią równania ruchu Newtona dla jednowymiarowego izotermicznego przepływu gazu.
W przedstawionych powyżej równaniach ciągłości oraz ruchu występują cztery zmienne: ciśnienie, gęstość, przepływ oraz temperatura, które są funkcjami czasu i zmiennej przestrzennej. Określenie wartości tych zmiennych będzie możliwe po dodaniu dodatkowych dwóch równań. Tymi równaniami są równanie stanu oraz równanie energii. Należy przy tym zauważyć, że jeśli przyjmiemy przemianę izotermiczną i odpowiadające jej równanie stanu uwzględnimy w równaniu ciągłości, to równanie (1.7) wraz z równaniem (1.11) tworzą model matematyczny nieustalonego izotermicznego przepływu gazu w gazociągu.
Równanie stanu
Równanie stanu opisuje wzajemną zależność pomiędzy parametrami stanu gazu (ciśnieniem, gęstością, temperaturą)
[1.12]
gdzie:
R - indywidualna stała gazowa, wyrażona w m2s-2K-1,
Z - współczynnik ściśliwości, zależny od wartości krytycznych temperatury i ciśnienia oraz ich wartości rzeczywistych,
T - temperatura gazu.
W przemianie izotermicznej (T = const) zależność pomiędzy ciśnieniem a gęstością przyjmuje postać
[1.13]
Przekształcając powyższe równanie oraz korzystając z definicji izotermicznej prędkości dźwięku w gazie, otrzymujemy
[1.14]
Równanie energii
Pierwsze prawo termodynamiki w odniesieniu do układu znajdującego się wewnątrz objętości kontrolnej jest postaci
[1.15]
gdzie:
Q - ilość ciepła dostarczonego do układu,
W - praca wykonana przez układ,
?E - zmiana energii układu.
Energia związana z masą układu wyrażana jest następująco:
[1.16]
gdzie:
U - energia wewnętrzna związana z oddziaływaniami molekularnymi w układzie,
(1/2)mw2 - energia kinetyczna układu,
mgz - energia potencjalna związana z położeniem układu w polu grawitacyjnym.
Po przekształceniach i pewnych uproszczeniach wynikających z definicji warunków, w jakich odbywa się przepływ, równanie energii w postaci różniczkowej przyjmie postać
[1.17]
gdzie:
u - energia wewnętrzna właściwa, tzn. wartość na jednostkę masy gazu;
- człon określający ilość ciepła wywiązującą się w płynie wewnątrz powierzchni kontrolnej w jednostce czasu;
- człon dotyczący lokalnej zmiany energii odpowiednio - wewnętrznej, kinetycznej i potencjalnej w czasie;
- człon opisujący strumień energii przez powierzchnię kontrolną;
- człon dotyczący zmiany strumienia energii wskutek konwekcji;
- człon reprezentujący strumień energii związany z pracą naprężeń normalnych.
W przypadku przepływu adiabatycznego oraz izotermicznego człon wyrażający ilość ciepła wywiązującego się w płynie na jednostkę czasu ?Q /dt przyjmuje odpowiednio wartości:
- ?Q /dt = 0 - przepływ adiabatyczny,
- ?Q /dt ? 0, T = const - przepływ izotermiczny.
Przepływ adiabatyczny dotyczy procesów szybkich, w których zjawiska związane z wymianą ciepła są pomijalne. Przepływy izotermiczne (temperatura gazu wzdłuż gazociągu nie ulega zmianie, T = const) dotyczą wolnych zmian przepływu, w których zjawiska związane z wymianą ciepła pomiędzy gazem a otaczającym gazociąg gruntem nie są uwzględnione. Równaniami opisującymi nieustalony, izotermiczny przepływ gazu w gazociągu jest równanie ciągłości (1.7) oraz równanie ruchu (1.11). W przypadku przepływu nieizotermicznego (T ? const) następuje wymiana ciepła między gazem a otoczeniem. W takim przypadku należy uzupełnić układ równań (1.5) i (1.11) o równanie energii (1.17) oraz równanie stanu (1.12).
Modele nieizotermiczne w odróżnieniu od izotermicznych pozwalają uwzględniać zmiany temperatury gazu w sieci. Stabilizacji temperatury w strumieniu służy dostatecznie duża długość gazociągu, przy niezmiennych warunkach termicznych panujących w gruncie. Jednak dla małych odległości od stacji przetłoczonych, gdzie następuje wzrost temperatury sprężanego gazu, zmiany temperatury gazu T(x) są wyraźne.
Przedstawione powyżej modele matematyczne nieustalonego nieizotermicznego i izotermicznego przepływu gazu wymagają dużych nakładów obliczeniowych. Należy pamiętać, że nie w każdym przypadku istnieje potrzeba korzystania z pełnego modelu, aby zapewnić oczekiwaną dokładność obliczeń. W praktyce szukamy modeli, które są kompromisem między dokładnością opisu zjawiska symulowanego a niezbędnym nakładem obliczeniowym wymaganym do jego rozwiązania.
Takie modele otrzymuje się poprzez odrzucenie niektórych członów z równań dokładnych, w oparciu o ocenę ich ilościowego wpływu na wynik obliczeń. Oznacza to, że wykorzystywany w procesie symulacji model nieustalonego przepływu gazu musi być odpowiedni do jasno zdefiniowanych granicznych wartości wybranych parametrów eksploatacji gazociągu. Analiza tych parametrów jest warunkiem koniecznym odpowiedniego wyboru modelu.
Przykładowe uproszczone modele nieustalonego przepływu gazu [3], [5]
Dokładny model matematyczny nieustalonego, izotermicznego przepływu gazu w gazociągu poziomym ma postać
[1.18]
Zakładając, że mamy do czynienia z gazociągiem płaskim, a zmiany prędkości przepływającego gazu są powolne, równanie ruchu uprościć można do postaci
[1.19]
W przypadku szybkich zmian prędkości przepływu gazu, spowodowanych na przykład dużymi zmianami obciążenia gazociągu lub nagłymi zmianami parametrów dostawy gazu, w równaniu ruchu koniecznym staje się uwzględnienie członu inercyjnego. Wówczas
[1.20]
Uzupełniając równanie (1.19) o równanie ciągłości w procesie izotermicznym (1.7), otrzymamy model nieustalonego izotermicznego przepływu gazu opisany układem dwóch równań różniczkowych cząstkowych typu parabolicznego
[1.21]
Jeśli z kolei przyjmiemy równanie ruchu (1.20), to otrzymamy
[1.22]
Układ (1.22) jest układem dwóch równań różniczkowych cząstkowych typu hiperbolicznego.
Dokładny model matematyczny nieustalonego, nieizotermicznego przepływu gazu w gazociągu poziomym ma postać
[1.23]
Przy założeniu, że wymiana ciepła może nastąpić jedynie jako przenikanie przez ściankę gazociągu lub jako efekt podłużnego przewodnictwa cieplnego (wszelkie możliwe reakcje chemiczne, dysocjację, jonizację, pochłanianie lub emisję promieniowania zaniedbujemy, przy czym gaz traktujemy jako ośrodek jednorodny o stałej przewodności cieplnej), z przeprowadzonej analizy wartości poszczególnych członów dokładnego równania energii dla określonych warunków eksploatacyjnych i geometrii gazociągu wynika, że równanie energii można uprościć do postaci [6]
[1.24]
Oznacza to, że uproszczona postać modelu nieustalonego nieizotermicznego przepływu gazu w rurociągu poziomym może być opisana poniższym układem równań
[1.25]
gdzie:
- lokalna zmiana gęstości w czasie,
- strumień masy przypadający na jednostkę powierzchni,
- zmiana ciśnienia gazu na długości,
- siły bezwładności w układzie,
- zmiana energii kinetycznej (zmiana pola prędkości),
- składowa sił grawitacji,
- siła od oporu hydraulicznego,
- zmiana ilości ciepła w strumieniu gazu na długości,
- rozprężanie gazu na długości,
- zmiana energii potencjalnej na długości,
- zmiana ilości ciepła w jednostce objętości gazu w czasie,
- przenikanie ciepła przez ściankę gazociągu.
Bardzo istotnym zagadnieniem jest odpowiedź na pytanie, w jakich warunkach eksploatacji gazociągu/systemu przesyłowego należy stosować model hiperboliczny, a w jakich model paraboliczny. Wiadomo, że model hiperboliczny ze względu na swoją strukturę dokładniej opisuje szybko zachodzące zmiany parametrów gazu spowodowane nagle zmieniającymi się warunkami brzegowymi. Oznacza to, że istnieje pewna graniczna częstotliwość zmian warunków brzegowych, powyżej której zaleca się ze względu na dokładność opisu zjawiska stosowanie modelu hiperbolicznego. Należy jednocześnie pamiętać, że to, jaka będzie reakcja gazociągu na zakłócenie, zależy nie tylko od jakości zakłócenia, ale także od pojemności akumulacyjnej gazociągu, czy też od średniego ciśnienia w gazociągu. Sformułowanie jasnych kryteriów doboru modeli do warunków pracy gazociągu jest bardzo istotne, jako że model paraboliczny jest łatwiejszy do numerycznego rozwiązania, istnieje wiele efektywnych obliczeniowo numerycznych algorytmów, a także, co jest nie bez znaczenia, znacznie łatwiej sformułować algorytm symulacji sieci o dowolnej strukturze z elementami nierurowymi (m.in. tłocznie, sprężarki, zawory, reduktory, stacje gazowe) o różnych wariantach pracy tych elementów. Modele paraboliczne aczkolwiek są znacznie łatwiejsze do numerycznego rozwiązywania i wykorzystywania do tworzenia programów symulacyjnych pozwalających na obliczanie sieci o dowolnych strukturach i dowolnej liczbie elementów nierurowych, jednak w pewnych warunkach mniej dokładnie opisują zmiany zachodzące w sieci w stosunku do modeli hiperbolicznych. Te ostatnie z kolei są bardziej pracochłonne obliczeniowo, a także trudniej je wykorzystać do stworzenia efektywnego symulatora sieci o dowolnej strukturze.
Algorytmy symulacji nieustalonego przepływu gazu w sieci o dowolnej strukturze
Algorytm symulacji dynamicznej oblicza wartości ciśnienia lub przepływu w wybranych punktach (węzłach) sieci w określonych chwilach. W porównaniu z algorytmem stanów ustalonych liczba danych wejściowych niezbędnych do wykonania obliczeń jest znacznie większa, ponieważ zarówno parametry obciążenia (odbiór gazu z sieci), jak również parametry zasilania systemu (parametry źródeł) są funkcją czasu. Rośnie w sposób znaczący w porównaniu z symulacją statyczną czas obliczeń, ponieważ:
- problem polega na rozwiązaniu układu równań różniczkowych cząstkowych lub zwyczajnych zamiast algebraicznych,
- symulacja dynamiczna dotyczy określonego przedziału czasu, a więc wielokrotnego rozwiązywania układu równań w tym przedziale czasu.
Wnioski
Dzięki algorytmom symulacji możemy znacznie zmniejszyć nakłady na wyposażenie pomiarowe systemów gazowniczych, a także, zakładając poprawność modeli, uzyskiwać pełną i dokładną informację o systemie. Algorytmy symulacji powinny umożliwiać realizację następujących funkcji:
- Przeprowadzanie symulacji sieci w całości lub wybranych jej fragmentów w trybie off-line w celu analizy reakcji sieci na określone wymuszenia lub zmiany jej konfiguracji.
- Symulacja on-line na bazie danych dostarczonych przez system telemetryczny.
- Ocena jakości pracy systemu telemetrycznego i urządzeń pomiarowych.
- Symulacja w celu monitorowania pracy sieci.
- Śledzenie zmian akumulacji gazu w systemie.
- Wydawanie zapewnień dostawy gazu przyszłym odbiorcom.
- Sprawdzenie jakości funkcjonowania sieci podczas prowadzenia remontu wybranej jej części.
- Określenie masy gazu traconej w systemie na skutek nieszczelności.
- Ocena możliwości przepustowych sieci.
- Śledzenie procesu mieszania się gazów o różnych składach.
- Lokalizacja nieszczelności.
Literatura do rozdziału pierwszego
1. Kołodziej W., Analiza matematyczna, Państwowe Wydawnictwo Naukowe, Warszawa 2020.
2. Krysicki W., Włodarski L., Analiza matematyczna w zadaniach, Część I, Część II, Państwowe Wydawnictwo Naukowe, Warszawa 1999.
3. Osiadacz A.J., Simulation and Analysis of Gas Networks, E.@F.N. Spon, London 1987.
4. Osiadacz A.J., Statyczna symulacja sieci gazowych, Fluid Systems, Warszawa 2001.
5. Osiadacz A.J., Symulacja i optymalizacja sieci gazowych, [w:] Vademecum Gazownika, Tom IV, Rozdział 6, SITPNIG, Kraków 2012.
6. Osiadacz A.J., Chaczykowski M., Comparison of isothermal and non-isothermal pipeline flow models, Chemical Engineering Journal, vol. 81, no 1-3, 2001, s. 41-51.
Rozdział 2.Podstawy rachunku macierzowego
2.1. Macierze oraz działania na macierzach
Macierzą A (dokładniej macierzą dwuskładnikową) nazywamy funkcję dwóch zmiennych, która każdej parze liczb naturalnych (i,j) (i = 1, 2, ..., m; j = 1, 2, ..., n) przyporządkowuje dokładnie jeden element ai,j będący liczbą rzeczywistą lub zespoloną, lub operatorem (np. różniczkowania), bądź wielomianem. Macierz zapisujemy w postaci tablicy
[2.1]
lub krócej jako A = [ai,j], i = 1, 2, ..., m; j = 1, 2, ..., n; Am×n. Element ai,j znajduje się w i-tym wierszu oraz w j-tej kolumnie macierzy.
Macierz prostokątną, której liczba kolumn jest równa jeden (n = 1), nazywamy macierzą kolumnową, którą można traktować jako m-wymiarowy wektor kolumnowy. Macierz prostokątną, której liczba wierszy wynosi jeden (m = 1), nazywamy macierzą wierszową, którą można traktować jako n-wymiarowy wektor wierszowy. Traktując kolumny macierzy jako m-wymiarowe wektory
( j = 1, 2, ..., n)
[2.2]
można macierz (2.1) zapisać jako
[2.3]
Jeśli wiersze macierzy potraktujemy jako n-wymiarowe wektory wierszowe postaci Ai = [ai1 ai2 ... ... ain] (i = 1, 2, ..., m), to macierz (2.1) zapiszemy jako
Wymiarem macierzy (2.1) nazywamy parę uporządkowaną zawierającą liczbę wierszy m i kolumn n, którą oznaczamy przez m × n.
Macierz, której liczba wierszy m jest równa liczbie kolumn n, nazywamy macierzą kwadratową. W przeciwnym wypadku mamy do czynienia z macierzą prostokątną.
Macierz kwadratową, w której różne od zera elementy, dii ? 0 (i = 1, 2, ..., n), występują tylko na głównej przekątnej, nazywamy macierzą diagonalną i oznaczamy symbolem diag(d11, d22, ..., dnn)
[2.4]
Jeżeli w macierzy diagonalnej wszystkie dii ? 0 są sobie równe, to taką macierz nazywamy macierzą skalarną. Szczególnym przypadkiem macierzy D jest macierz jednostkowa In stopnia n, określona wzorem
[2.5]
Wprowadzając symbol Kroneckera
[2.6]
macierz jednostkową zapiszemy w postaci I = [?ij].
Macierz, której wszystkie elementy są równe zero, nazywamy macierzą zerową.
Macierz kwadratowa, której elementy spełniają warunek aij = aji dla i = 1, 2, ..., m; j = 1, 2, ..., n, nazywamy macierzą symetryczną.
Macierz kwadratową nazywamy trójkątną dolną - L oraz górną - U, jeśli wszystkie elementy położone poniżej/powyżej głównej przekątnej są równe zero.
Macierz trójkątna dolna stopnia n jest postaci
a macierz trójkątna górna
Podstawowe działania na macierzach
Dwie macierze A = [aij]mn i B = [bij]mn są równe wtedy i tylko wtedy, gdy ich odpowiednie elementy są sobie równe, tzn. A = B, gdy [aij] = [bij] dla i = 1, 2, ..., m; j = 1, 2, ..., n.
Sumą (różnicą) dwóch macierzy A = [aij]mn i B = [bij]mn o jednakowych wymiarach jest macierz C o elementach równych sumom (różnicom) odpowiednich elementów macierzy A i B, tzn. C = A ? B, gdy cij = [aij] ? [bij] dla i = 1, 2, ..., m; j = 1, 2, ..., n.
Z powyższej definicji wynikają następujące własności:
(A + B) + C = A + (B + C )
A + B = B + A
A + 0mn = A
Iloczyn A × k lub k × A macierzy A przez liczbę k jest macierzą o elementach równych odpowiednim elementom macierzy A pomnożonym przez k.
Niech będą dane dwie macierze
w których liczba kolumn macierzy A jest równa liczbie wierszy macierzy B.
Iloczynem macierzy A i B nazywamy macierz , której elementy określone są zależnością dla i = 1, 2, ..., m; k = 1, 2, ..., q.
Na przykład iloczynem macierzy przez jest macierz .
Mnożenie macierzy na ogół nie jest przemienne, czyli AB ? BA.
W związku z nieprzemiennością mnożenia macierzy wyróżnia się mnożenie prawostronne macierzy A przez macierz B, czyli C = AB, oraz mnożenie lewostronne, czyli D = BA.
Macierze A i B spełniające warunek AB = BA nazywamy macierzami przemiennymi.
Ogólne własności mnożenia macierzy są następujące:
Iloczyny macierzy trójkątnych dolnych są także macierzami trójkątnymi dolnymi, a iloczyny macierzy trójkątnych górnych są macierzami trójkątnymi górnymi.
Transpozycja macierzy
Macierzą transponowaną AT macierzy A = [aij] o wymiarach m × n nazywamy macierz o wymiarach m × n, którą otrzymano w rezultacie zamiany wierszy z kolumnami. Elementy [aij] macierzy A zastąpiono elementami [aji] macierzy AT.
Dla dowolnej macierzy prostokątnej
macierz transponowana ma postać:
Wektor wierszowy xT = (x1, x2, ..., xn) jest transpozycją pewnego wektora kolumnowego. Można wykazać, że:
Jeżeli A jest macierzą kwadratową, dla której A = AT, to macierz A jest macierzą symetryczną, tzn. aij = aji dla wszystkich i oraz wszystkich j. W macierzy symetrycznej elementy symetryczne względem głównej przekątnej są równe.
Iloczyn dwóch macierzy symetrycznych A i B jest symetryczny tylko pod warunkiem, że AB = BA.
Macierz symetryczna nazywa się dodatnio określoną, jeśli związana z nią forma kwadratowa xTAx spełnia warunek xTAx > 0 dla każdego rzeczywistego x ? 0.
Macierzą ortogonalną nazywamy macierz kwadratową An×n spełniającą równość AAT = AAT= In.
Macierze blokowe
Dowolną macierz A można przedstawić jako macierz blokową, zbudowaną z pewnej liczby macierzy o mniejszych wymiarach
[2.7]
gdzie Aij jest macierzą o wymiarach mi × ni (i = 1, 2, ..., q ? m; j = 1, 2, ..., p ? n).
W przypadku kiedy macierze na przekątnej Aii są kwadratowe, macierz A również musi być kwadratowa, a pi = qi (i = 1, 2, ..., n). Dodawanie i mnożenie takich macierzy blokowych wykonuje się, traktując, że bloki są liczbami.
Na przykład, dla C = AB istnieje zależność
[2.8]
Macierzą blokowo-przekqtniową nazywamy macierz, którą można zapisać w postaci
A = diag(A11, A22, ..., Ann)
przy czym macierze Aii muszą być kwadratowe.
Analogicznie definiuje się macierz blokowo-trojkątną.
Wyznaczniki
Wyznacznik macierzy kwadratowej A stopnia n oznacza się symbolem det(A) lub | A|:
[2.9]
W przypadku macierzy prostokątnej A = [aij]m×n skreślamy pewna liczbę wierszy i kolumn, tak aby z elementów nieskreślonych macierzy prostokątnej utworzyć macierz kwadratową. Wyznacznik tak otrzymanej macierzy nazywamy minorem macierzy A.
Najwyższy spośród stopni minorów macierzy A różnych od zera nazywamy rzędem macierzy A i oznaczamy r(A). Macierz stopnia m × n ma rząd spełniający nierówność 0 ? r(A) ? minor(m,n).
Rząd macierzy jest równy maksymalnej liczbie liniowo niezależnych wierszy lub kolumn. Rząd kolumnowy macierzy A to liczba liniowo niezależnych kolumn macierzy A. Rząd wierszowy macierzy A to liczba liniowo niezależnych wierszy macierzy A.
Zawsze rząd kolumnowy i wierszowy macierzy A są sobie równe i będziemy je nazywać rzędem macierzy A (oznaczenie r(A)).
W szczególnym przypadku, jeśli r = m = n, macierz A jest macierzą nieosobliwą.
Przykład
Określmy rząd macierzy
Tworzymy minory stopnia 3
Minor stopnia 2 wynosi
Oznacza to, że rząd macierzy wynosi 2.
Minorem Mij elementu aij wyznacznika det A stopnia n nazywamy wyznacznik stopnia n - 1 powstały przez usunięcie z wyznacznika det A i-tego wiersza oraz j-tej kolumny.
Wartość liczbową wyznacznika detA stopnia n można obliczyć, stosując poniższe wzory:
det A = ai1Ai1 + ai2Ai2 + ... + ainAin (1 ? i ? n)
[2.10]
det A = a1j A1j + a2j A2j + ... + anj Anj (1 ? j ? n)
[2.11]
gdzie A1k (k = 1, 2, ..., n) oznacza wyznacznik stopnia n - l, który powstaje przez skreślenie z macierzy A pierwszego wiersza i k-tej kolumny.
Własności wyznaczników:
- wyznacznik macierzy trójkątnej jest równy iloczynowi elementów głównej przekątnej det(L) = l11 - l22, ..., - lnn, det(U) = u11- u22, ..., - unn,
- transponowanie macierzy kwadratowej nie zmienia wyznacznika tej macierzy detAT = detA,
- zamiana dwóch wierszy (kolumn) macierzy kwadratowej zmienia wartość wyznacznika tej macierzy na przeciwną,
- jeśli w macierzy dwa wiersze (dwie kolumny) są identyczne , to wyznacznik tej macierzy wynosi zero,
- jeżeli macierz B powstaje z macierzy A przez pomnożenie wszystkich elementów pewnego wiersza (kolumny) przez liczbę ?, to detB = ?detA,
- dodanie do wiersza (kolumny) wielokrotności innego wiersza (kolumny) nie zmienia wyznacznika tej macierzy,
- wyznacznik iloczynu dwu macierzy kwadratowych jest równy iloczynowi wyznaczników tych macierzy det(A - B) = detA detB.
Obliczanie wyznaczników
1. Macierz drugiego stopnia
2. Macierz trzeciego stopnia:
- wzór Laplace'a
- schemat Sarrusa
- reguła Chio
przy założeniu, że a11 ? 0
Odwracanie macierzy
Macierzą odwrotną A-1 macierzy kwadratowej stopnia n nazywamy macierz spełniającą następujący warunek: AA-1 = A-1A = I. W przypadku iloczynu macierzy obowiązuje zależność (AB)-1 = A-1B-1.
Macierz odwrotna A-1 istnieje wtedy i tylko wtedy, gdy A jest macierzą nieosobliwą, tzn. jeśli det(A) ? 0. W przeciwnym wypadku macierz określa się jako macierz osobliwą.
Macierz odwrotną obliczamy zgodnie z zależnością
gdzie AD - macierz dołączona.
Dla każdego elementu aij (i, j = 1, 2, ..., n) macierzy A możemy obliczyć dopełnienie algebraiczne Aij = (-1)i+jMij, gdzie Mij jest minorem, który powstał przez skreślenie i-tego wiersza oraz j-tej kolumny w macierzy A.
Aby wyznaczyć macierz dołączoną AD macierzy A, należy:
- wyznaczyć macierz dopełnień algebraicznych dla danej macierzy kwadratowej A, tzn. zastąpić każdy element aij macierzy A odpowiadającym mu dopełnieniem algebraicznym,
- wyznaczyć macierz transponowaną macierzy dopełnień algebraicznych.
Przykłady
Aby wyznaczyć macierz odwrotną macierzy , obliczamy: detA = (3 - 1 - 2 + 2 - 1 - 1 + 2 - 0 - 0) - (1 - 1 - 0 + 0 - 1 - 3 + 2 - 2 - 2) = 0, detA = 0, co oznacza, że macierz A-1 nie istnieje.
Aby wyznaczyć macierz odwrotną macierzy , sprawdzamy, że detB = 8, co oznacza, że B-1 istnieje. Obliczamy macierz dopełnień algebraicznych dla macierzy B.
Ostatecznie
Otrzymane rozwiązanie jest poprawne ponieważ BB-1 = I.
Przestrzeń liniowa wektorowa
Wektorem x w przestrzeni n-wymiarowej nazywamy uporządkowany zbiór składający się z n liczb rzeczywistych lub zespolonych: . Liczby x1, x2, ..., xn nazywamy współrzędnymi wektora.
Przykłady wektorów:
- wektor na płaszczyźnie jest wektorem przestrzeni dwuwymiarowej,
- w macierzy A o wymiarze m × n każdy jej wiersz jest wektorem n-wymiarowym, a każda kolumna wektorem m-wymiarowym.
Dwa wektory x oraz y w przestrzeni m-wymiarowej, , są równe, jeśli xi = yi dla i = 1, 2, ..., m.
Sumą/różnicą dwóch wektorów nazywamy wektor przestrzeni m-wymiarowej o współrzędnych odpowiednio xi + yi lub xi - yi (i = 1, 2, ..., m).
W przypadku dodawania wektorów słuszne są następujące równości:
x + y = y + x (przemienność)
(x + y) + z = x + ( y + z) (łączność)
Iloczynem wektora przez liczbę ? nazywamy wektor .
Dowolny zbiór wektorów n-wymiarowych, w którym są określone operacje dodawania wektorów oraz mnożenia wektora przez liczbę, przy czym rezultaty tych działań należą do tego zbioru, nazywamy liniową przestrzenią wektorową i oznaczamy przez En.
Wektor y = c1x1 + c2x2 + ... + ckxk nazywany jest kombinacją liniową wektorów x1, x2, ..., xk. Wektory te są nazywane liniowo niezależnymi, jeśli równość c1x1 + c2x2 + ... + ckxk = 0 zachodzi tylko pod warunkiem, że c1 = c2 = ... = ck = 0. W przeciwnym wypadku wektory te określa się jako liniowo zależne.
Wymiarem liniowej przestrzeni wektorowej nazywamy maksymalną liczbę wektorów liniowo niezależnych. Dowolny zbiór n wektorów liniowo niezależnych n-wymiarowej przestrzeni tworzy bazę tej przestrzeni.
Każdy wektor x ? Rn można więc wyrazić jako
x = a1y1 + a2y2 + ... + anyn
[2.12]
gdzie a1, a2, ..., an nazywa się współrzędnymi tego wektora względem bazy y1, y2, ..., yn.
W przestrzeni liniowej rzeczywistej jest określony iloczyn skalarny, jeśli każdej parze wektorów x oraz y należących do tej przestrzeni przyporządkowana jest liczba rzeczywista oznaczona przez (x, y). Liniowa przestrzeń rzeczywista, w której zdefiniowany jest iloczyn skalarny, nazywa się przestrzenią euklidesową.
W przestrzeni euklidesowej aksjomaty iloczynu skalarnego zdefiniowane są następująco:
(x, x) > 0, jeśli x ? 0,
(x, y) = ( y, x)
(? x, y) = (x, ? y) = ? (x, y) ? - liczba rzeczywista
(x + y, z) = (x, z) + ( y, z), (x, y + z) = (x, y) + (x, z)
Iloczyn skalarny wektorów x oraz y jest liczbą określoną wzorem
a w zapisie macierzowym
Długością (modułem) wektora w przestrzeni euklidesowej nazywamy liczbę
2.2. Wartości własne i wektory własne
Jeśli dla danej liczby ? i wektora x ? 0 spełniona jest równość
Ax = ?x
[2.13]
to ? nazywa się wartością własną macierzy A, a x - wektorem własnym, odpowiadającym wartości ?. Zadanie polega na wyznaczeniu takich wartości parametru ?, przy których powyższe równanie ma rozwiązanie niezerowe.
Równanie (2.13) można zapisać w postaci układu jednorodnych równań algebraicznych liniowych
[2.14]
lub w postaci macierzowo-wektorowej
[2.15]
gdzie A - ?I jest macierzą charakterystyczną dla macierzy A = [aij]n×n.
Wyznacznikiem charakterystycznym macierzy A nazywamy wyznacznik
[2.16]
Równanie (2.16) nazywane jest równaniem charakterystycznym macierzy A, natomiast pierwiastki równania charakterystycznego nazywane są wartościami własnymi macierzy A. Zbiór wartości własnych (?1, ?2, ..., ?n) macierzy kwadratowej A nazywamy widmem tej macierzy. Wielomian charakterystyczny ma postać
Współczynniki ki (i = 1, 2, ..., n) są funkcją elementów macierzy A:
- suma minorów głównych stopnia drugiego macierzy A,
ki - suma minorów głównych stopnia trzeciego itd., kn = det A.
W przypadku przestrzeni R2 macierz A - ?I przyjmuje postać
czyli:
det(A - ?I) = (a11 - ?)(a22 - ?) - a12a21 = ?2 - (a11 + a22)? + (a11a22 - a12a21) = 0
tzn. w tym przypadku ? jest pierwiastkiem wielomianu drugiego stopnia.
Przykład
Niech
czyli
Wielomian charakterystyczny jest postaci
D(?) = det (A - ?I)= (1 - ?)(2 - ?)(3 - ?)
co oznacza, że ?1 = 1, ?2 = 2, ?3 = 3.
Niech
x - wektor własny odpowiada ?1,
y - wektor własny odpowiada ?2,
z - wektor własny odpowiada ?3.
Aby znaleźć x, rozwiązujemy układ równań
czyli
Dla x3 = 0, przyjmując x1 = 1, otrzymujemy x2 = -2, tzn.
Wartości wektora y obliczamy z układu równań
czyli
Ponieważ y3 = 0, y1 = 0, stąd y2 =1, tzn.
Po rozwiązaniu poniższego układu równań otrzymujemy
Twierdzenie Cayleya-Hamiltona
Ważnym twierdzeniem w teorii macierzy jest twierdzenie Cayleya-Hamiltona, które mówi, że każda kwadratowa macierz A spełnia swoje równanie charakterystyczne, tzn. |A -?I| = 0. Jeżeli D(?) jest wielomianem charakterystycznym macierzy A = |a|n×n, to wówczas D(A) = 0.
Jeśli wielomian charakterystyczny macierzy A jest postaci
to wówczas jest spełnione równanie macierzowe
[2.17]
Przykład
Niech
Wtedy
Twierdzenie Cayleye'a-Hamiltona można wykorzystać do znalezienia macierzy odwrotnej. Mnożąc obustronnie równanie (2.17) przez A-1, otrzymujemy
Ostatecznie
Przykład
Aby znaleźć macierz odwrotną do macierzy
obliczamy
Ostatecznie
