Metody numeryczne dla inżynierów, Skrypty z Matematyka

Pobierz Metody numeryczne dla inżynierów i więcej Skrypty w PDF z Matematyka tylko na Docsity!

M ETODY NUMERYCZNE

DLA INŻYNIERÓW

(notatki do wykładu)

Wrocław, maj 2018

(^4) Spis treści

7.4. Metody ekstrapolacyjno-interpolacyjne ......................................................... 61

8. Literatura ................................................................................................. 63

Skorowidz ...................................................................................................... 65

1. Wstęp

Niniejszy skrypt zawiera opis głównych zagadnień prezentowanych na wykładzie Metody numeryczne dla inżynierów , który jest przeznaczony dla studentów kierunku Automatyka i Robotyka na Wydziale Elektrycznym Politechniki Wrocławskiej. Metody numeryczne są podstawowym narzędziem analitycznym w rękach współczesnego inżyniera i stąd też nietrudno znaleźć wyczerpującą literaturę na ten temat o różnym stopniu zaawansowania - niektóre propozycje podane są w końcowej części pracy. Każde jednak ujęcie tego tematu jest przeznaczone dla określonego czytelnika, o odpowiednim stopniu przygotowania i z myślą o specyficznym zastosowaniu prezentowanych metod. Głównym celem niniejszego opracowania jest prezentacja podstawowych metod numerycznych stosowanych w obliczeniach w elektrotechnice. Zakłada się, że Czytelnik zna podstawowy kurs algebry i analizy matematycznej. Wymagana jest również podstawowa znajomość zasad tworzenia algorytmów obliczeniowych. Wykonanie prezentowanych przykładów obliczeniowych wymaga również elementarnej znajomości korzystania z komputerów. Z wykładem związane są ćwiczenia laboratoryjne, w trakcie których są praktycznie ilustrowane zagadnienia przedstawiane na wykładzie. Podstawowym narzędziem programowym, stosowanym do opisu poszczególnych procedur obliczeniowych, jak i do obliczeń w laboratorium komputerowym jest MATLAB. Program ten jest stosowany tu zarówno do formułowania i sprawdzania prostych algorytmów numerycznych, jak i do rozwiązywania bardziej złożonych zagadnień z wykorzystaniem gotowych procedur. Pakiet programowy MATLAB, jak wiele innych tego typu programów przeznaczonych do rozwiązywania zadań inżynierskich, zawiera sporą liczbę gotowych procedur numerycznych, które są dostępne w postaci pojedynczych instrukcji. Można zatem spytać, jaki jest cel dodatkowego wykładu na ten temat, skoro wystarczy się zapoznać z instrukcją obsługi odpowiedniego programu komputerowego. Jednak każdy użytkownik tego typu oprogramowania specjalistycznego natrafia na problemy związane z wyborem odpowiednich procedur (często do rozwiązania tego samego zadania można zastosować różne algorytmy), interpretacji błędów, dokładności wyników, rozwiązywania zagadnień niestandardowych, czy wreszcie rozumienia i interpretacji tekstu instrukcji. Ważna jest także umiejętność formułowania modeli matematycznych analizowanych zjawisk, które pozwalają określić poszukiwane parametry lub zależności między nimi. W takich przypadkach wymagana jest niekiedy pogłębiona znajomość zagadnień analizy numerycznej. W praktyce inżynierskiej metody numeryczne są narzędziem służącym do formułowania i rozwiązywania praktycznych zagadnień obliczeniowych, a także do przekształcenia znanych modeli ciągłych do adekwatnych postaci dyskretnych. Z tego punktu widzenia metody numeryczne są tu traktowane jako wygodny i wydajny sposób rozwiązywania zadań inżynierskich. Przygotowanie i rozwiązanie takiego typu zadań wiąże się zazwyczaj z wykonaniem następujących działań:

Wstęp 7

Złożoność obliczeniowa algorytmu jest związana z liczbą operacji numerycznych, które prowadzą do uzyskania wyniku. Jest zrozumiałe, że spośród różnych algorytmów, które zapewniają poprawne rozwiązanie, należy wybierać te, które charakteryzują się małą złożonością obliczeniową. Jest to szczególnie istotne w układach sterowania, gdzie pełny cykl obliczeń numerycznych musi być wykonany w czasie określonym przez okres pomiędzy kolejnymi pomiarami wielkości wejściowych. Uwarunkowanie zadania jest cechą metod numerycznych, która określa możliwość uzyskania poprawnych wyników przy stosowaniu dowolnych danych wejściowych z odpowiednio zdefiniowanego zbioru. Jeśli analizowany algorytm służy do rozwiązania zadania y  w ( x ), to stopień uwarunkowania zadania można

mierzyć za pomocą ilorazu w ( x   x ) x ( x )/  x. Niekiedy używa się też terminu

czułość zadania. Mówi się, zatem, że zadanie jest dobrze uwarunkowane lub źle uwarunkowane. W pierwszym przypadku zadanie jest stabilne względem danych wejściowych, co oznacza, że rozwiązanie w sposób ciągły zależy od dokładności

danych wejściowych tak, że dla  x  0 jest  y  0. W przypadku złego

uwarunkowania zadania, możliwość uzyskania poprawnego rozwiązania zależy od wartości danych wejściowych. Cecha ta jest wykorzystywana do odpowiedniej korekcji zadań źle uwarunkowanych, które nie mogą być inaczej rozwiązane. W wielu przypadkach algorytm zastosowany do rozwiązania zadania dobrze uwarunkowanego (stabilnego) może być niestabilny. Stabilność numeryczna algorytmu odnosi się do możliwości uzyskania określonej dokładności obliczeń. Algorytm jest stabilny numerycznie, gdy zwiększając dokładność obliczeń można z dowolną dokładnością określić dowolne z istniejących rozwiązań.

(^10) Liniowe układy równań

2.2. Metoda eliminacji Gaussa

Powyższy przykład z układem trzech równań liniowych można rozwiązać stosując metodę, która jest zbliżona do tradycyjnej metody 'szkolnej'. Polega ona na kolejnej eliminacji zmiennych. Korzysta się przy tym z prostych działań, takich jak: mnożenie obu stron równania przez stałą wartość lub dodawanie równań stronami.

W rozważanym przypadku (2.1), zmienna x 1 może być wyeliminowana z

drugiego równania przez odjęcie od niego równania pierwszego pomnożonego przez współczynnik 3 /10  0,3. Podobnie można postąpić z trzecim równaniem: w tym przypadku pierwsze równanie przed odjęciem go od równania trzeciego należy pomnożyć przez współczynnik 5 /10  0,5. Po wykonaniu tych operacji

otrzymamy następującą postać równania (2.1):

2 3

2 3

1 2

x x

x x

x x (2.5)

które ma następującą formę macierzową:

3

2

1

x

x

x (2.6)

Z ostatnich dwóch równań można najpierw określić x (^) 3 przez eliminację zmiennej x 2 z trzeciego równania. Można to uzyskać przez dodanie drugiego równania po

jego pomnożeniu przez współczynnik 2,5 / 0,1  25. Ostatecznie otrzymamy:

3

2

1

x

x

x (2.7)

Zauważmy, że z ostatniego równania (ostatni wiersz) można już bezpośrednio określić zmienną x 3. Ten etap obliczeń nazywa się etapem eliminacji zmiennych.

Poczynając teraz od ostatniego równania (ostatniego wiersza w zapisie macierzowym) można otrzymać kolejne rozwiązania. Jest to postępowanie odwrotne. Zatem, w celu uzyskania wartości wszystkich niewiadomych wykonujemy następujące działania:

x 1  

2

^ 

x 

75.^8

1

  ^ 

x  ^  

Liniowe układy równań 11

Powyższe operacje można zapisać dla ogólnego przypadku. W tym celu rozpatrzmy ogólną postać równania (2.3), gdzie:

11 12 1 21 22 2

1 2

n n

n n nn

a a a a a a

a a a

 ^ 

A

1 2

n

b b

b

 ^ 

b 

1 2

n

x x

x

 ^ 

x 

które można zapisać w postaci następującego układu równań

11 1 12 2 1 1 21 1 22 2 2 2

1 1 2 2

n n n n

n n nn n n

a x a x a x b a x a x a x b

a x a x a x b

Stosując pierwszy krok eliminacji w odniesieniu do (2.9) otrzymamy układ równań, w których poczynając od drugiego z nich, wyeliminowana jest zmienna x 1 :

(2) 11 1 12 2 1 1 (2) (2) (2) 22 2 2 2

(2) (2) (2) 2 2

n n n n

n nn n n

a x a x a x b a x a x b

a x a x b

gdzie:

11

22 (^2 ) 22 12 21 a

a  a  a a , 11

23 (^2 ) 23 13 21 a

a  a  a a , ..., 11

2 (^2 ) 2 1 21 a

a a a a n ^ n  n ,

11

32 (^2 ) 32 12 31 a

a  a  a a , 11

33 (^2 ) 33 13 31 a

a  a  a a , ..., 11

3 ( 2 ) 3 1 31 a

a a a a n ^ n  n , ...,

11

1 2 12

( 2 ) (^2) a a a a an n ^ n  , 11

1 3 13

( 2 ) (^3) a a a a an n ^ n  , ..., 11

1 1

( 2 ) a

a a a an nn  nn  n

oraz

11

2 (^2 ) 2 1 21 a

b  b  b a , 11

3 ( 2 ) 3 1 31 a

b  b  b a , ..., 11

( 2 ) 1 1 a

b b b an n  n 

W ostatnim kroku tej procedury układ równań ma następującą postać:

11 1 12 2 1 1 1 (2) (2) (2) (2) 22 2 2 2 2

( 1) ( 1) ( 1) 1 1 2 1 2 1 ( ) ( )

n n n n n n n n

n n n n n n n n n n nn n n

a x a x a x a x b a x a x a x b

a x a x b a x b

      

Liniowe układy równań 13

nn

n n (^) a

b x 

for k :  n  1 to 1 step  1 do kk

n

jk

k kj j k (^) a

b a x x

 

W powyższym algorytmie zakłada się, że w pierwszym etapie nie jest tworzona nowa macierz U , natomiast tworzona macierz trójkątna jest zapisywana na miejscu macierzy A.

Jak widać, w operacjach arytmetycznych ważną rolę odgrywają elementy leżące na przekątnej macierzy współczynników równania. Przez nie są dzielone odpowiednie równania w pierwszym etapie eliminacji zmiennych. Także w wyniku dzielenia uzyskuje się kolejne rozwiązania na etapie podstawiania zmiennych. Rozwiązanie staje się nieosiągalne, gdy któryś z tych elementów diagonalnych jest równy zero (wówczas macierz parametrów jest osobliwa). Również przy małych wartościach elementów diagonalnych można spodziewać się dużych błędów (gdyż występuje dzielenie przez małą liczbę, która - z racji reprezentacji dyskretnej - może być przedstawiona niedokładnie. Aby tego uniknąć stosuje się modyfikację metody, która polega na tak zwanym częściowym wyborze elementu wiodącego. W tym celu, przed eliminacją kolejnej zmiennej (etap wprzód), spośród równań pozostających do rozpatrzenia (poniżej danego wiersza) wybiera się to, które ma w redukowanej kolumnie (w pierwszej niezerowej) największą wartość i zamienia się go z danym równaniem. Odpowiedni algorytm zostanie pokazany w następnym rozdziale.

Optymalne metody rozwiązywania układów równań liniowych powinny przewidywać takie uporządkowanie równania, aby macierz A była diagonalnie dominującą. Oznacza to, że moduły elementów na przekątnej są nie mniejsze od sumy modułów pozostałych elementów w tym samym wierszu (wówczas jest to macierz diagonalnie dominująca kolumnowo), co można zapisać następująco

 



n

ki

ii (^) k ki a a 1

, i  1 , 2 ,..., n

2.3. Metoda rozkładu LU

Załóżmy, że kwadratowa macierz współczynników równania A zostanie przedstawiona w postaci iloczynu dwóch macierzy trójkątnych:

A  LU (2.15)

gdzie:

(^14) Liniowe układy równań

21

1,1 1, ,1 ,2 , 1

n n n n n n

l

l l l l l

  

 ^ 

L ^ ^ ^ ^ ,

11 12 1, 1 1 22 2, 1 2

1, 1 1, ,

n n n n

n n n n n n

u u u u u u u

u u u

 

  

 ^ 

U ^ ^ ^ 

Załóżmy, że znane są macierze L , (^) U dla danej macierzy A. Wówczas równanie (2.15) można zapisać w następującej formie

LUx  b (2.16)

Wektor x można określić w dwóch etapach, rozwiązując kolejno następujące równania

Lz  b (2.17)

Ux  z (2.18)

Ze względu na trójkątną strukturę macierzy L oraz U , równania (2.17)- (2.18) można rozwiązać bezpośrednio przez odwrotne podstawianie, jak w metodzie Gaussa. Wymaga to wykonania n^2 operacji mnożenia i dzielenia, a więc tyle, ile potrzeba na pomnożenia macierzy przez wektor. Dużą oszczędność uzyskuje się wówczas, gdy równanie (2.16) trzeba rozwiązać dla różnych wartości wektora b. Należy zauważyć, że macierze L oraz U mogą być zapisane w jednej macierzy P   L \ U , gdyż elementy diagonalne macierzy L są zawsze równe 1, więc nie

muszą być pamiętane.

Wektor x można określić za pomocą następującego algorytmu

for i : 1 to n do (^) 





1

1

i

m

zi bi pimzm {rozwiązanie równania (2.17)}

for i : n to 1 step  1 do (^) jj

n

m j

x (^) j : zj pjmxm / p 1 

   

{rozwiązanie równania (2.18) }

Algorytm rozkładu LU można łatwo wyznaczyć na podstawie związku (2.15). Na przykład, dla n  3 macierz A wyraża się w następujący sposób za pomocą współczynników macierzy L oraz (^) U :

31 11 31 12 32 22 31 13 32 23 33

21 11 21 12 22 21 13 23

11 12 13

31 32 33

21 22 23

11 12 13

l u l u l u l u l u u

l u l u u l u u

u u u

a a a

a a a

a a a

A

Z powyższego przedstawienia można określić sposób obliczania elementów macierzy L oraz U :

1. u 11 (^)  a 11

l (^) 21  a 21 / u 11

(^16) Liniowe układy równań

for i : 1 to n do d (^) j : 0 ;

{ główny algorytm }

for k : 1 to n do

begin { wybór elementu głównego } b : 0 ; for j : k to n do begin







1

1

k

m

a (^) jk ajk ajmamk ;

if a (^) jk  b then

begin b : ajk ; w : j end;

end; if b  0 then begin err : 1 ; halt end; { brak rozwiązania } { przestawienie wierszy } if w  k then begin for j : 1 to n do

begin b : akj ; akj : awj ; a (^) wj : b end;

s : d k ; d (^) k : dw ; d (^) w : s

end; {obliczenie u (^) ki }

for i : k to n do  

 

1 1

k ki ki m km mi a a a a ;

{obliczenie l (^) jk } for j :  k  1 to n do a (^) jk : ajk / akk ;

end ;

Jeśli wynik tego algorytmu jest stosowany łącznie z algorytmem rozwiązywania równania (2.16), to wektor b należy uszeregować zgodnie z indeksami zawartymi w wektorze przestawień d :

bi  bd i , i  1 , 2 ,..., n ,

gdzie: wektor b  b (^) i może być bezpośrednio użyty w algorytmie (2.17).

Algorytm rozwiązywania układu równań liniowych może być stosowany do odwracania macierzy. Zauważmy, że

Liniowe układy równań 17

1



 ^ 

AA

zatem

1 1 ( 1) ( 1) ( 1) (1) (2) ( n ) (1) (2) ( n )

A A (^)  1 1  (^1)   a a  a  (2.21)

gdzie

(^1) ( (^) i ) jest wektorem kolumnowym ( n  1 ) z jedynką na i -tej pozycji i zerami w

pozostałych miejscach,

( 1 ) ()

 a (^) i jest i -tą kolumną poszukiwanej macierzy A ^1.

Można zauważyć, że a ((^  i )^1 ) jest rozwiązaniem równania

( )

( 1 ) Aa (^) ( i )  (^1) i

zatem, w celu obliczenia macierzy A ^1 należy rozwiązać n równań typu (2.22). W przedstawionych metodach wymaga to tylko jednokrotnego rozkładu macierzy A (na macierz trójkątna lub na macierze LU). Złożoność obliczeniowa takiego algorytmu jest z grubsza równa trzykrotnej złożoności rozwiązania pojedynczego układu równań liniowych.

2.4. Iteracyjne metody rozwiązywania układu równań liniowych

Przedstawione powyżej metody eliminacji nie uwzględniają różnych właściwości macierzy współczynników, które w metodach iteracyjnych mogą prowadzić do uproszczenia obliczeń, co jest szczególnie ważne w zadaniach o dużych rozmiarach. Ma to miejsce, na przykład, w przypadku macierzy o silnie dominującej przekątnej, gdy wiele elementów leżących poza przekątną ma małą wartość lub są to elementy zerowe. Można w takim przypadku założyć, że wszystkie elementy leżące na przekątnej macierzy współczynników równania są różne od zera. W taki przypadku równanie (2.3) można zapisać w następującej postaci:





n

ji

j

i ij j ii

i b a x

a

x

1

, i  1 , 2 ,..., n (2.23)

Przy zadanych wartościach początkowych poszukiwanych niewiadomych zdefiniowanych przez wektor x , kolejne przybliżenia można uzyskać zgodnie z algorytmem iteracyjnym. Metody iteracyjne sprowadzają się do poszukiwania rozwiązania układu równań o postaci

f (^) i ( x 1 , x 2 ,..., xn ) 0 , i  1 , 2 ,..., n (2.24)

Liniowe układy równań 19

11 21 22 1 1,1 1,2 1, 1 ,1 ,2 , 1

n n n n n n n n nn

a a a

a a a a a a a

    

 ^ 

A ^ ^ ^ ^ ,

12 1, 1 1 2, 1 2 2 1,

n n n n

n n

a a a a a

a

 



A    

Algorytm iteracyjnego poszukiwania rozwiązania wynika bezpośrednio z (2.30). W następujących po sobie krokach określane jest przybliżenie kolejnej zmiennej po uwzględnieniu uzyskanych przybliżeń poprzednich zmiennych:

1 1 12 2 11 13 3 11 1 11 1 11 1 1 2 21 1 22 23 3 22 2 22 2 22

1 1 1 1 1 1 2 2 , 1 1

k k k k n n k k k k n n

k k k k n n nn n nn n n n nn n nn

x a x a a x a a x a b a x a x a a x a a x a b a

x a x a a x a a x a b a

  

     

co może być zapisane w następującej ogólnej postaci



 

  n ji

k j ii

i ij j

k j ii

ij ii

k i i (^) a x

a x a

a a

b x 1

1 1

Warunki zbieżności procesu iteracyjnego związanego z algorytmem Gaussa- Seidla mogą być określone na podstawie badania równania uzyskanego z (2.30)

x k^ ^1   A 1 ^1 A 2 x k  A 1 ^1 b (2.33)

Można pokazać (patrz rozdział dotyczący iteracyjnego rozwiązywania układów równań nieliniowych), że warunek zbieżności procesu określonego przez (2.33) jest określony przez wartości własne macierzy  A 1  1 A 2. Dostatecznym i wystarczającym

warunkiem zbieżności metody jest to aby moduły wszystkich wartości własnych tej macierzy były mniejsze od jedności. Jest to równoważne następującemu warunkowi odnoszącemu się do współczynników macierzy A

  



n

ji

ii (^) j ij a a 1

i  1 , 2 ,..., n (2.34)

co oznacza, że rozwiązanie iteracyjne jest możliwe, jeśli moduły elementów diagonalnych są większe od sumy modułów wszystkich pozostałych elementów w wierszu macierzy. Większość zagadnień spotykanych w technice spełnia ten warunek. Niekiedy należy wcześniej odpowiednio przekształcić wyjściowy układ równań.

Ostatecznie, metoda Gaussa-Seidla iteracyjnego rozwiązywania układów równań liniowych przybiera formę następującego algorytmu.

1. Uporządkować wyjściowy układ n równań tak, aby w macierzy współczynników A największe co do modułu elementy znalazły się na przekątnej, co jest określone następującym warunkiem a (^) ii  aiji  j , i  1 , 2 ,..., n , j  1 , 2 ,..., n

(^20) Liniowe układy równań

2. Przyjąć warunki początkowe   x (^) 0  x 20 x^03  xn^0 
3. Powtarzać proces iteracyjny (2.33) dla k  1 , 2 , aż spełniony zostanie warunek

i max 1 , 2 ,..., nxi^ k ^1  xi^ k ^ 

gdzie  - założona dokładność obliczeń. Metody iteracyjne stosowane są zazwyczaj do rozwiązywania dużych układów równań, w których wiele współczynników ma wartość zerową (są to tak zwane równania z macierzami rzadkimi). Wówczas można oczekiwać mniejszej złożoności obliczeniowej takiego podejścia niż stosowanie metod skończonych. Metody iteracyjne są także stosowane do poprawiania (zwiększania dokładności) wyników uzyskanych w rezultacie stosowania metod skończonych.