METODY NUMERYCZNE
DLA INŻYNIERÓW
(notatki do wykładu)
[e-mail]
Wrocław, marzec 2012
Spis Treści
1. Wstęp ......................................................................................................... 5
2. Liniowe układy równań .......................................................................... 9
2.1. Wprowadzenie ..................................................................................................... 9
2.2. Metoda eliminacji Gaussa................................................................................. 10
2.3. Metoda rozkładu LU ......................................................................................... 13
2.4. Iteracyjne metody rozwiązywania układu równań liniowych ................... 17
3. Rozwiązywanie równań nieliniowych ............................................... 21
3.1. Zagadnienia jednowymiarowe ........................................................................ 21
Metoda prostej iteracji....................................................................................... 21
Metoda połowienia ............................................................................................ 22
Metoda Newtona ............................................................................................... 23
Metoda siecznych .............................................................................................. 23
Metody wielokrokowe: algorytm Aitkena .................................................... 24
3.2. Rozwiązywanie układów równań nieliniowych .......................................... 25
Metoda Newtona-Raphsona ............................................................................ 25
Metoda siecznych .............................................................................................. 27
4. Interpolacja .............................................................................................. 29
4.1. Wprowadzenie ................................................................................................... 29
4.2. Wielomian interpolacyjny Newtona ............................................................... 30
4.3. Numeryczne różniczkowanie funkcji dyskretnej ......................................... 34
5. Aproksymacja ......................................................................................... 35
5.1. Wprowadzenie ................................................................................................... 35
5.2. Aproksymacja średniokwadratowa ................................................................ 36
5.3. Filtr wygładzający .............................................................................................. 40
5.4. Filtr różniczkujący ............................................................................................. 42
5.5. Przykład obliczeniowy...................................................................................... 43
5.6. Metoda Najmniejszych Kwadratów z wykorzystaniem rozkładu
macierzy według wartości szczególnych  SVD ........................................... 44
6. Całkowanie numeryczne ...................................................................... 47
6.1. Wprowadzenie ................................................................................................... 47
6.2. Metoda Simpsona .............................................................................................. 47
7. Numeryczne rozwiązywanie równań różniczkowych
zwyczajnych ............................................................................................ 49
7.1. Wprowadzenie ................................................................................................... 49
7.2. Metody jednokrokowe ...................................................................................... 51
Metoda Eulera .................................................................................................... 51
Metoda trapezów............................................................................................... 53
Metody Rungego-Kutty.................................................................................... 53
Dokładność metody .......................................................................................... 54
Stabilność metody ............................................................................................. 55
7.3. Metody wielokrokowe ...................................................................................... 58
Metody Geara..................................................................................................... 58
Niejawna metoda Rungego-Kutty .................................................................. 59
4
Spis treści
7.4. Metody ekstrapolacyjno-interpolacyjne ......................................................... 59
8. Literatura ................................................................................................. 61
Skorowidz ...................................................................................................... 63
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 znalez
ć 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ń:
6
Wstęp
- określenie modelu matematycznego analizowanego zjawiska lub opis stanu
obserwowanego systemu;
- wybranie (opracowanie) odpowiedniej metody obliczeń numerycznych;
- analiza i weryfikacja poprawności przyjętego modelu oraz wykonanych
obliczeń.
W niniejszym wykładzie będziemy się zajmować głównie drugim z
wymienionych działań. A
ączy się ono z podaniem sposobu (algorytmu)
numerycznego rozwiązania postawionego zadania. W obliczeniach prowadzonych z
zastosowaniem metod numerycznych należy się liczyć ze specyfiką stosowanych
narzędzi. Liczby reprezentowane w komputerze są przedstawiane z ograniczoną
dokładnością, która zależy od liczby bitów użytych do ich zapisu. Wynikające stąd
błędy najczęściej nie mają znaczenia w dalszym wykorzystaniu wyników obliczeń.
Niekiedy jednak wartość błędów powstających w poszczególnych etapach obliczeń
jest tak duża, że kontynuowanie obliczeń staje się niemożliwe (przekroczenie
zakresu) lub uzyskane wyniki zawierają niedopuszczalne błędy.
Można wyróżnić następujące cztery z
ródła błędów, które ograniczają dokładność
końcowych wyników:
1. błędy w danych wejściowych
2. przybliżony model zjawiska
3. błędy aproksymacji modelu
4. błędy zaokrągleń
Błędy danych wejściowych leżą poza procesem obliczeń, jednak stosowanie
odpowiednich procedur może prowadzić do redukcji ich wpływu na wynik (na
przykład, wygładzanie danych pomiarowych). Problem ten łączy się zatem z drugim
z wymienionych z
ródeł błędów. Należy jednak podkreślić, że błędy danych
wejściowych, w ogólnym przypadku, są nieusuwalne.
Błędny lub przybliżony model analizowanego procesu wynika z uproszczeń
przyjmowanych w trakcie formułowania modelu matematycznego zjawiska lub
opisu stanu. Wynika to z potrzeby redukcji złożoności modelu, która jest
przyjmowana w sposób świadomy lub z braku odpowiednich danych, tak że
analizowany proces jest przedstawiany w sposób uproszczony.
Błąd metody jest związany z tym, że poprawny model jest aproksymowany za
pomocą uproszczonych formuł, w których dodatkowo mogą być stosowane
przybliżone dane raz parametry. Typowym przykładem tego typu podejścia jest
aproksymowanie zależności różniczkowych za pomocą funkcji różnicowych.
Błędy zaokrągleń wynikają ze skończonej długości reprezentacji liczb w
komputerze. Jeśli błędy te mają charakter przypadkowy a nie systematyczny, to
sumaryczny błąd statystyczny nawet długiej serii obliczeń jest zazwyczaj mały.
Systematyczne błędy zaokrągleń mogą jednak prowadzić do szybko rosnącej
niedokładności obliczeń. Znanym z
ródłem takich błędów jest operacja odejmowania
bliskich sobie liczb. Jeśli w algorytmie szeregowo powtarzane są takie działania, to
szybko następuje niedopuszczalna kumulacja błędów. Typowy przykład jest
związany z umieszczeniem takiej różnicy w mianowniku jakiegoś wyrażenia.
Poprawność i efektywność algorytmów obliczeniowych jest określana za pomocą
różnych parametrów. Oto niektóre z nich.
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 z
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ń z
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ń.
2. Liniowe układy równań
2.1. Wprowadzenie
Zagadnienie rozwiązywania układów równań liniowych jest podstawowym
problemem w metodach numerycznych. Metod rozwiązywania tego zagadnienia jest
wiele, a wybór tej czy innej metody zależy od rodzaju zadania, oczekiwanej
dokładności i środków technicznych będących w dyspozycji (szybkość procesora
oraz objętość pamięci).
Załóżmy, że mamy układ trzech równań z trzema niewiadomymi:
10x1 - 7x2 = 6
- 3x1 + 2x2 + 6x3 = 4 (1.1)
5x1 - x2 + 5x3 = 3
Równanie to można zapisać w następującej postaci macierzowej:
10
Ą# - 7 0 x1 6
ń#Ą# ń# Ą# ń#
ó# Ą# ó#4Ą#
= (1.2)
ó#- 3 2 6Ą#ó#x2 Ą# ó# Ą#
Ą#ó#
ó# 5 -1 5Ą#ó#x3 Ą# ó#3Ą#
Ł# Ś#Ł# Ś# Ł# Ś#
Przechodząc do postaci ogólnej mamy:
Ax = b (1.3)
gdzie: A - macierz kwadratowa ( n � n ); w tym przypadku n = 3,
x - wektor niewiadomych ( n �1),
b - wektor współczynników prawej strony ( n �1).
Jeśli wyznacznik macierzy det(A) `" 0 , to rozwiązanie można przedstawić w
następującej postaci
x = A-1b (1.4)
Można pokazać, że poszukiwanie rozwiązania równania (1.3) w postaci (1.4)
prowadzi do algorytmu o dużej złożoności obliczeniowej, co jest związane z
odwracaniem macierzy A . Już zastosowanie reguł 'ręcznego' rozwiązywania układu
równań (1.1) redukuje około n razy liczbę niezbędnych mnożeń potrzebnych do
uzyskania wyniku. Poniżej przedstawimy niektóre najczęściej stosowane metody
rozwiązywania równania (1.3).
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 (1.1), zmienna x1 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 (1.1):
10x1 - 7x2 = 6
- 0.1x2 + 6x3 = 5.8 (1.5)
2.5x2 + 5x3 = 0
które ma następującą formę macierzową:
10
Ą# - 7 0 x1 6
ń#Ą# ń# Ą# ń#
ó# ó#5.8Ą#
0 - 0.1 6Ą#ó#x2 Ą# = (1.6)
ó# Ą#ó# Ą# ó# Ą#
ó# 0 2.5 5Ą#ó#x3 Ą# ó# 0 Ą#
Ł# Ś#Ł# Ś# Ł# Ś#
Z ostatnich dwóch równań można najpierw określić x3 przez eliminację zmiennej
x2 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:
10
Ą# - 7 0 x1 6
ń#Ą# ń# Ą# ń#
ó# ó#5.8Ą#
0 - 0.1 6Ą#ó#x2 Ą# = (1.7)
ó# Ą#ó# Ą# ó# Ą#
ó# 0 0 155Ą#ó#x3 Ą# ó#145Ą#
Ł# Ś#Ł# Ś# Ł# Ś#
Zauważmy, że z ostatniego równania (ostatni wiersz) można już bezpośrednio
określić zmienną x3 . 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:
145 29
x1 = =
155 31
1 29 1 179.8 - 174 - 58
�#5.8 - 6 ś# �# ś#
x2 = = =
ś# ź# ś# ź#
- 0.1 31 - 0.1 31 31
�# # �# #
1 29 5.8 1 18.6 - 40.6 - 22
�#6 ś#
x1 = - 0 - 7 = =
ś# ź#
10 31 3.1 # 10 3.1 31
�#
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 (1.3), gdzie:
a11 a12 ... a1n b1 x1
Ą# ń# Ą# ń# Ą# ń#
ó#a a22 ... a2n Ą# ó#b Ą# ó#x Ą#
21 2 2
ó# Ą# ó# Ą# ó# Ą#
A = , b = , x = (1.8)
ó# M M ... M Ą# ó# M Ą# ó# M Ą#
ó#a an2 ... ann Ą# ó#b Ą# ó#x Ą#
Ł# n1 Ś# Ł# n Ś# Ł# n Ś#
które można zapisać w postaci następującego układu równań
a11x1 + a12x2 + ... + a1n xn = b1
a21x1 + a22x2 + ... + a2n xn = b2
(1.9)
M M M M M
+ an1x1 + an2x2 + ... + ann xn = bn
Stosując pierwszy krok eliminacji w odniesieniu do (1.9) otrzymamy układ
równań, w których poczynając od drugiego z nich, wyeliminowana jest zmienna x1 :
(2)
a11x1 + a12x2 + ... + a1n xn = b1
(2 ( (
+ a22)x2 + ... + a22)xn = b22)
n
(1.10)
M M M M M
( (2 (
+ an2)x2 + ... + ann)xn = bn2)
2
gdzie:
a21 (2 a21 ( a21
(2
a22) = a22 - a12 , a23) = a23 - a13 , ..., a22) = a2n - a1n ,
n
a11 a11 a11
a31 (2 a31 ( a31
(2
a32) = a32 - a12 , a33) = a33 - a13 , ..., a32) = a3n - a1n , ...,
n
a11 a11 a11
an1 ( an1 (2 an1
(
an2) = an2 - a12 , an2) = an3 - a13 , ..., ann) = ann - a1n
2
a11 3 a11 a11
oraz
a21 ( a31 ( an1
(
b22) = b2 - b1 , b32) = b3 - b1 , ..., bn2) = bn - b1
a11 a11 a11
W ostatnim kroku tej procedury układ równań ma następującą postać:
a11x1 + a12x2 + ... + a1n xn a1n xn = b1
(2 ( ( (
+ a22)x2 + ... + a22)xn a22)xn = b22)
n n
(1.11)
M M M M M M
( 1 ( 1 ( -
+ ann-n) x2 + ann-n)x2 bnn11)
-1 -1 -1 -
(n (
+ ann)xn = bnn)
12
Liniowe układy równań
(n-1 (n-1
ann-1) ( ( ( - ann-1)
(n (n ( -1
gdzie: ann) = ann-1) - ann1n) , bnn) = bnn-1) - bnn11) .
- ( -1
- ( -1
ann1n) ann1n)
- -1 - -1
W ten sposób, po wykonaniu procedury eliminacji zmiennych pierwotne równanie
przekształca się do postaci z górną trójkątną macierzą U :
Ux = b (1.12)
gdzie:
u11 u12 ... u1n
Ą# ń#
ó#
0 u22 ... u2n Ą#
ó# Ą#
U =
ó# M M ... M Ą#
ó#
0 0 ... unn Ą#
Ł# Ś#
Niewiadomą xn wyznacza się z równania określonego przez ostatni wiersz:
bn
xn = (1.13)
unn
Dalej, znając niewiadome xn , xn-1, xk +1 z k -tego równania obliczamy:
n
bk - x
"u
kj j
j=k +1
xk = (1.14)
ukk
przy czym uwzględniane są odpowiednio przekształcone współczynniki wektora b .
Ostatecznie otrzymujemy następujący algorytm rozwiązywania układów równań
liniowych metodą eliminacji Gaussa.
{eliminacja zmiennych}
for i := 2 to n do
for k := i to n do
begin
ak
ż#
akl
�# - ai-1,l ,i-1 dla l = i, i + 1, ..., n
akl :=
ai-1,i-1
�#
�#
0 dla l = 1, 2, ..., i - 1
�#
ak
bk := bk - bi-1 ,i-1
ai-1,i-1
end;
end;
{odwrotne podstawianie}
Liniowe układy równań 13
bn
xn =
ann
n
bk - x
"akj j
j=k +1
for k := n -1 to 1 step -1 do xk =
akk
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
aii e" aki , i = 1, 2,..., n
"
k =1
k `"i
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 (1.15)
gdzie:
14
Liniowe układy równań
u11 u12 ... u1,n-1 u1n
Ą# ń#
ó# Ą#
1
Ą# ń#
ó#
u22 ... u2,n-1 u2n Ą#
ó# Ą#
l21 1 0
ó# Ą#
ó# Ą#
ó# Ą#
L = ó# M M ... M MĄ# , U = M M ... M M
ó# Ą#
ó#l ln-1,2 ... 1 Ą#
ó# Ą#
n-1,1
ó# Ą#
0 ... un-1,n-1 un-1,n Ą#
ó#
ó# Ą#
ln,1 ln,2 ... ln,n-1 1Ś#
Ł#
ó# Ą#
ó# ... un,n Ś#
Ą#
Ł#
Załóżmy, że znane są macierze L , U dla danej macierzy A . Wówczas równanie
(1.15) można zapisać w następującej formie
LUx = b (1.16)
Wektor x można określić w dwóch etapach, rozwiązując kolejno następujące
równania
Lz = b (1.17)
Ux = z (1.18)
Ze względu na trójkątną strukturę macierzy L oraz U , równania (1.17)- (1.18)
można rozwiązać bezpośrednio przez odwrotne podstawianie, jak w metodzie
Gaussa. Wymaga to wykonania n2 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 (1.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
i-1
for i := 1 to n do zi := bi - pimzm {rozwiązanie równania (1.17)}
"
m=1
n
�# ś#
ś# ź#
for i := n to 1 step -1 do x := z - p xm ź# / p {rozwiązanie równania (1.18) }
j j " jm jj
ś#
m= j+1
�# #
Algorytm rozkładu LU można łatwo wyznaczyć na podstawie związku (1.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 :
a11 a12 a13 u11 u12 u13
Ą# ń# Ą# ń#
ó# Ą# ó# Ą#
ó#a21 ó#l21u11 l21u12 + u22 l21u13 + u23 Ą#
A = a22 a23Ą# =
ó# Ą# ó# Ą#
ó#a a32 a33 Ą# ó#l u11 l31u12 + l32u22 l31u13 + l32u23 + u33 Ą#
Ł# 31 Ś# Ł# 31 Ś#
Z powyższego przedstawienia można określić sposób obliczania elementów
macierzy L oraz U :
Liniowe układy równań 15
1. u11 = a11
l21 = a21 /u11
l31 = a31 / u11
2. u12 = a12
u22 = a22 -l21u12
l32 = (a32 -l31u12 )/u22
3. u13 = a13
u23 = a23 -l21u13
u33 = a33 -l31u13 -l32u23
Widać, że w każdym z trzech kroków ( n = 3) najpierw są obliczane elementy
macierzy U , a następnie elementy macierzy L w danej kolumnie. Dla ogólnego
przypadku można to zapisać w postaci następującego algorytmu
{warunki początkowe - inicjalizacja macierzy: L = 1 ( n � n ), U = 0 ( n � n )}
for k :=1 to n do
begin
k -1
for i := k to n do uki := aki - "
lkmumi ;
m=1
�#a k -1l umk ś#/u
for j := k +1 to n do l := - "
;
ś# ź#
jk jk jm kk
�# m=1 #
end;
Algorytm ten jest nazywany algorytmem Gaussa-Banachiewicza [8].
Podobnie jak w przypadku algorytmu Gaussa, dla poprawienia skuteczności i
dokładności algorytmu LU można stosować wybór maksymalnego elementu
głównego w kolumnie. W tym celu należy porównać ze sobą wyrazy k -tej kolumny
macierzy A leżące na i poniżej głównej przekątnej ( k d" j d" n ):
k -1
p := a - "
l umk k d" j d" n (1.19)
j jk jm
m=1
i wybrać spośród nich największy co do modułu. Odpowiadający mu wiersz należy
przestawić z rozpatrywanym k -tym wierszem macierzy A . Procedura ta nie
prowadzi do znacznego skomplikowania algorytmu, gdyż wyrażenie (1.18) jest
fragmentem głównego algorytmu i tak musi być obliczone.
Załóżmy, że elementy macierzy L oraz U będą zapisane na odpowiednich
miejscach macierzy A (macierz ta nie zostanie zachowana), a elementy wektora
d ={di} określają numery wierszy macierzy A zgodnie z przestawieniem
wynikającym z wyboru maksymalnego elementu głównego. Przeprowadzone
rozważania prowadzą wówczas do następującego algorytmu rozkładu LU z
wyborem maksymalnego elementu głównego w kolumnie.
16
Liniowe układy równań
{warunki początkowe}
err := 0 ;
for i :=1 to n do d := 0 ;
j
{główny algorytm}
for k :=1 to n do
begin
{wybór elementu głównego}
b:= 0 ;
for j := k to n do
begin
k -1
a := a - amk ;
"a
jk jk jm
m=1
if a > b then
jk
begin b:= a ; w:= j end;
jk
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 ; awj := b end;
s := dk ; dk := dw ; dw := s
end;
{obliczenie uki }
k -1
for i := k to n do aki := aki - "
akmami ;
m=1
{obliczenie l }
jk
for j := k +1 to n do a := a /akk ;
jk jk
end;
Jeśli wynik tego algorytmu jest stosowany łącznie z algorytmem rozwiązywania
równania (1.16), to wektor b należy uszeregować zgodnie z indeksami zawartymi w
wektorze przestawień d :
bi = bd , i =1,2,...,n ,
i
gdzie: wektor b ={bi} może być bezpośrednio użyty w algorytmie (1.17).
Liniowe układy równań 17
Algorytm rozwiązywania układu równań liniowych może być stosowany do
odwracania macierzy. Zauważmy, że
1 0 L 0
Ą# ń#
ó#0 1 L 0Ą#
ó# Ą#
AA-1 = (1.20)
ó#M M L MĄ#
ó#0 0 L 1Ą#
Ł# Ś#
zatem
-
A-1 = A-1[1(1) 1(2) L 1(n)]= [a(-1) a(-1) L a(n1)] (1.21)
(1) (2) ( )
gdzie
1(i) jest wektorem kolumnowym ( n �1) z jedynką na i -tej pozycji i zerami w
pozostałych miejscach,
-
a(i)1) jest i -tą kolumną poszukiwanej macierzy A-1 .
(
-
Można zauważyć, że a(i)1) jest rozwiązaniem równania
(
-
Aa(i)1) = 1(i) (1.22)
(
zatem w celu obliczenia macierzy A-1 należy rozwiązać n równań typu (1.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 (1.3) można zapisać w następującej postaci:
�# ś#
n
ś# ź#
1
xi = - xj , i = 1,2,...,n (1.23)
"aij
aii ś#bi j=1 ź#
ś# ź#
j`"i
�# #
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
18
Liniowe układy równań
fi (x1, x2,..., xn ) = 0 , i = 1,2,...,n (1.24)
który jest równoważny (1.3).
Ogólny schemat iteracyjnego rozwiązywania układu n równań można zapisać
następującą zależnością
xij+1 = xij + j�ij , i = 1,2,...,n (1.25)
gdzie j jest numerem kroku iteracji,
j jest wielkością kroku iteracji,
�ij parametrem określającym 'kierunek' iteracji,
przy założonych początkowych wartościach xi0 , i = 1,2,...,n .
W przypadku układu równań liniowych, odpowiednie metody iteracyjne są
tworzone na podstawie przedstawienia równania (1.3) w następującej postaci
x = Cx + d (1.26)
skąd kolejne przybliżenia rozwiązania są określane zgodnie z równaniem
xk +1 = Cxk +d (1.27)
Zgodnie z tym algorytmem, równanie (1.23) można zapisać w następującej formie
iteracyjnej:
�# ś#
n
ś# ź#
1
xik+1 = - xk , i = 1,2,...,n (1.28)
"aij
aii ś#bi j=1 j ź#
ś# ź#
j`"i
�# #
Zależność ta jest znana jako iteracyjna metoda Jakobiego rozwiązywania równań
liniowych.
Poszczególne metody różnią się sposobem wyboru kroku iteracji  oraz
parametru � . Omówimy poniżej pewną modyfikację metody Jakobiego, znaną jako
metoda Gaussa-Seidla.
W metodzie Gaussa-Seidla kolejne przybliżenie rozwiązania równania (1.3)
określa się zgodnie z następującym podstawieniem
k k k k
a11x1 +1 + a12x2 + a13x3 + L+ a1n xn = b1
k k k k
a21x1 +1 + a22x2 +1 + a23x3 + L+ a2n xn = b2
(1.29)
L L L L L L
k k k k
an1x1 +1 + an2x2 +1 + an3x3 +1 + L+ ann xn +1 = bn
co można zapisać w następującej formie macierzowej
A1xk +1 + A2xk = b (1.30)
Liniowe układy równań 19
gdzie
a11
Ą# ń# a12 ... a1,n-1 a1n
Ą# ń#
ó# Ą#
ó# Ą#
ó# Ą#
ó#
a21 a22 0
... a2,n-1 a2n Ą#
ó# Ą#
ó# Ą#
ó# Ą#
ó#M Ą#
A1 = M M ... M M , A2 = M ... M M
ó# Ą#
ó# Ą#
ó# Ą#
ó# Ą#
0 ... an-1,n Ą#
n-1,1
ó#a an-1,2 ... an-1,n-1 Ą#
ó#
ó# Ą#
ó# Ą#
ó# an,1 an,2 ... an,n-1 ann Ś# Ł#
Ą# ...
ó# Ą#
Ś#
Ł#
Algorytm iteracyjnego poszukiwania rozwiązania wynika bezpośrednio z (1.30).
W następujących po sobie krokach określane jest przybliżenie kolejnej zmiennej po
uwzględnieniu uzyskanych przybliżeń poprzednich zmiennych:
k k k k
x1 +1 = -a12x2 / a11 - a13x3 / a11 - L- a1n xn / a11 + b1 / a11
k k k k
x2 +1 = -a21x1 +1 / a22 - a23x3 / a22 - L- a2n xn / a22 + b2 / a22
(1.31)
L L L L L L
k k k k +1
xn +1 = -an1x1 +1 / ann - an2x2 +1 / ann - L- an,n-1xn-1 / ann + bn / ann
co może być zapisane w następującej ogólnej postaci
bi i-1 aij +1 n aij
xik +1 = - "
xk - "
xk (1.32)
aii j=1aii j j=i+1 aii j
Warunki zbieżności procesu iteracyjnego związanego z algorytmem Gaussa-
Seidla mogą być określone na podstawie badania równania uzyskanego z (1.30)
-1 -1
xk +1 = -A1 A2xk + A1 b (1.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 (1.33) jest
-1
określony przez wartości własne macierzy -A1 A2 . 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
aii > aij i =1,2,...,n (1.34)
"
j=1
j`"i
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.
20
Liniowe układy równań
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
aii > aij i`" j , i =1,2,...,n , j =1,2,...,n
2. Przyjąć warunki początkowe
0 0 0
{x} ={x2 x3 L xn}
0
3. Powtarzać proces iteracyjny (1.33) dla k =1,2,L aż spełniony zostanie
warunek
max xik +1 - xik <�
i=1,2,...,n
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.
3. Rozwiązywanie równań nieliniowych
3.1. Zagadnienia jednowymiarowe
Załóżmy, że dana jest funkcja f (x) rzeczywistego argumentu x . Celem naszych
działań jest określenie rozwiązania następującego równania
f (x) = 0 (1.1)
to znaczy, określenie wartości zmiennej x , dla których spełniona jest zależność (1.1).
Należy zauważyć, że w ogólnym przypadku zadanie to nie jest proste, gdyż ze
względu na nieliniowość funkcji f (x) nie jest nawet wiadomo ilu rozwiązań można
oczekiwać. Nie ma ogólnych, jednoznacznych metod rozwiązywania takich zadań.
Znane są natomiast metody przybliżone, które opierają się na poszukiwaniu
rozwiązań w drodze kolejnych iteracyjnych przybliżeń.
Metoda prostej iteracji
Zapiszmy równanie (1.1) w następującej postaci
x = g(x) (1.2)
Iteracyjne rozwiązanie równania (1.2) polega na wykonaniu następujących działań
xk +1 = g(xk ) (1.3)
przy warunkach początkowych: x0 = x0 .
Powstaje oczywiście pytanie, czy ciąg wartości xk uzyskany w wyniku
stosowania procedury (1.3) prowadzi do rozwiązania, to znaczy, czy metoda jest
stabilna. Dowodzi się, że warunek zbieżności można zapisać następująco. Dla
dowolnie wybranej zmiennej � zachodzi nierówność
g(x) - g(� ) d" K x - � (1.4)
gdzie K < 1.
Jeśli warunek (1.4) jest spełniony, to algorytm (1.3) nazywa się odwzorowaniem
zawężającym, które prowadzi do rozwiązania. Warunek ten w wielu przypadkach
nie jest spełniony i różne metody iteracyjnego rozwiązywania równania (1.1) biorą
się stąd, żeby tak wyrazić równanie (1.1) w formie (1.2), aby poszerzyć obszar
zbieżności rozwiązania i przyśpieszyć proces tego rozwiązania. W ogólnym
przypadku odwzorowanie (1.3) można zapisać następująco
xk +1 = Ś(xk ) (1.5)
22
Rozwiązywanie równań nieliniowych
przy czym, funkcja Ś , znana jako funkcja iteracyjna, jest tak dobrana, że jeśli x' jest
rozwiązaniem równania (1.1), to Ś(x' ) = x' .
Metoda połowienia
Metoda połowienia (metoda bisekcji) wywodzi się z obserwacji, że jeśli na granicach
przedziału [a, b] funkcja f (x) ma różne znaki, to wewnątrz przedziału znajduje się
przynajmniej jedno miejsce zerowe tej funkcji. Z kolei strategia poszukiwania
kolejnego, bliższego rozwiązania polega na wskazaniu w tym celu punktu, leżącego
w środku tego właśnie przedziału. W ten sposób otrzymujemy następujący
algorytm.
{warunki początkowe}
x := a ; y := b ;
fx := f (x) ; fy := f ( y) ; { fx oraz fy powinny mieć różne znaki }
{pętla iteracyjna}
while abs(x - y) > � do
begin
{połowienie}
z := (x + y) / 2 ;
fz := f (z) ;
if sign( fz) = sign( fx) then
begin
p := x ; x := z ; z := p ;
end;
else
begin
p := y ; y := z ; z := p ;
end;
end;
Można zauważyć, że w przypadku cyfrowej reprezentacji liczb, w każdej iteracji
połowienia dokładność rozwiązania wzrasta o jeden bit. Algorytm jest zatem zbieżny
dosyć wolno, chociaż przy poprawnym wyborze początkowego przedziału, zawsze
prowadzi do rozwiązania. Jest on często stosowany jako procedura, która prowadzi
do rozwiązania w skrajnych sytuacjach, gdy zawodzą inne metody.
Rozwiązywanie równań nieliniowych 23
Metoda Newtona
Znaczne przyspieszenie procesu iteracyjnego można uzyskać, jeśli odpowiednio
dobierze się funkcję iteracyjną Ś w (1.5). W tym celu można zastąpić nieliniową
funkcję f (x) w pobliżu rozwiązania (to jest w pobliżu zera) za pomocą jej
rozwinięcia w szereg Taylora
(� - x0 )2
' ''
f (� ) = 0 = f (x0 ) + f (x0 )(� - x0 ) + f (x0 ) +
2!
(1.6)
(� - x0 )k
(k )
... + f (x0 + � ) + � (� - x0 )
k!
Pozostawiając tylko dwa pierwsze wyrazy rozwinięcia (przybliżenie liniowe)
otrzymujemy
'
0 H" f (x0 ) + f (x0 )(� - x0 ) (1.7)
oraz
f (x0 )
'
� H" x0 - , jeśli f (x) `" 0 ,
'
f (x0 )
co w ogólności prowadzi do następującej procedury iteracyjnej
f (xk )
xk +1 = xk - , (1.8)
'
f (xk )
która jest znana jako metoda Newtona rozwiązywania równań nieliniowych.
Można pokazać, że Metoda Newtona dla pierwiastków jednokrotnych ma
przynajmniej zbieżność kwadratową, co odnosi się do stopnia przybliżenia do
rozwiązania w kolejnych iteracjach.
Metoda siecznych
Jeśli w metodzie Newtona zastąpić różniczkowanie funkcji za pomocą wyrażenia
różnicowego, to otrzymamy przybliżenie metody Newtona, które ze względu na
interpretację graficzną jest znane jako metoda siecznych. Przybliżone
różniczkowanie funkcji f (x) może być określone następująco
f (xk ) - f (xk -1)
'
f (xk ) H" (1.9)
xk - xk -1
co, po podstawieniu do (1.9), prowadzi do następującego algorytmu
f (xk )(xk - xk -1)
xk +1 = xk - (1.10)
f (xk ) - f (xk -1)
24
Rozwiązywanie równań nieliniowych
jeśli tylko f (xk ) - f (xk -1) `" 0 .
Metoda siecznych jest w wielu przypadkach wygodniejsza do stosowania
(szczególnie w tych przypadkach, gdy nie ma możliwości określenia pochodnej
funkcji f (x) ), jednak jest ona słabiej zbieżna.
Zauważmy, że powyższe metody mogą być stosowane jedynie wówczas, gdy
spełniony jest warunek o różnej od zera wartości mianownika odpowiedniego
wyrażenia (1.8) lub (1.10). Poprawnie sformułowany algorytm powinien
uwzględniać to i w przypadku, gdy wartość ta jest odpowiednio mała, powinna być
proponowana inna wersja algorytmu.
Metody wielokrokowe: algorytm Aitkena
Dane jest równanie nieliniowe o postaci
f (x) = 0 (1.11)
Metoda prostej iteracji poszukiwania wartości x , dla której spełnione jest
równanie (1.11) polega na przekształceniu go do postaci
x = g(x) (1.12)
dla której można sformułować następującą regułę iteracyjną
xk +1 = g(xk ) (1.13)
z warunkami początkowymi: x0 = x0 .
Algorytm (1.13) prowadzi do rozwiązania, gdy proces iteracyjny jest zbieżny.
Zbieżność jest zapewniona, gdy spełniony jest następujący warunek. Dla dowolnie
wybranej zmiennej � zachodzi nierówność
g(x) - g(� ) d" K x - � (1.14)
gdzie K < 1.
Aby rozszerzyć obszar zbieżności i przyspieszyć zbieżność procesu iteracyjnego
można stosować jego korekcję według metody Aitkena. Jej idea polega na
zastąpieniu problemu rozwiązania równania (1.11) przez zagadnienie poszukiwania
zer funkcji, utworzonej z kolejnych wyników prostej iteracji:
h(xk)= xk - xk-1 = 0 (1.15)
gdzie zmienne xk oblicza się według (1.13).
Problem sprowadza się zatem do określenia sposobu korekcji metody prostej iteracji
w celu uzyskania rozwiązania procesu (1.15). Ponieważ funkcja h(xk) jest dostępna w
postaci numerycznej, więc rozwiązania (1.15) można poszukiwać za pomocą metody
siecznych:
h(xk) xk - xk-1
xk+1 = xk - = xk - , (1.16)
p
"h(xk) (xk+1 - xk)-(xk - xk-1)
"(xk) (xk - xk-1)
przy czym:
Rozwiązywanie równań nieliniowych 25
"h(xk)= (xk+1 - xk)-(xk - xk-1)= h(xk+1)- h(xk),
"(xk)= (xk - xk-1)= h(xk).
Korekcja jest zatem dokonywana na podstawie trzech kolejnych wartości xk-1, xk,
oraz xk+1, przybliżenia, uzyskanych według metody prostej iteracji zgodnie z
następującą regułą:
2
(xk +1 - xk )
xk +1 = xk - (1.17)
p
xk +2 - 2xk +1 + xk
Wynik tej korekcji przyjmuje się w charakterze kolejnego przybliżenia rozwiązania:
xk +1 = xk +1, po czym następują znów dwa kroki procedury (1.13) do kolejnej korekcji
p
(1.17). W ten sposób uzyskuje się algorytm o następującej postaci.
1. Przyjąć warunki początkowe
x0 = x0 , k = 0 - numer kroku iteracji
2. Wykonać dwa kroki prostej iteracji
yk = g(xk ) , zk = g( yk )
3. Skorygować wynik:
2
(yk - xk )
"k =
zk - 2 yk + xk
xk +1 = xk - "k
4. Jeśli abs("k ) > eps , k = k + 1, przejdz
do 2
3.2. Rozwiązywanie układów równań nieliniowych
Układ równań nieliniowych może być w ogólnym przypadku zapisany
następująco
f1(x1, x2,..., xn )
Ą# ń#
ó#
f2 (x1, x2,..., xn )Ą#
ó# Ą#
f (x) = = 0 (1.18)
ó# ... Ą#
ó#
fn (x1, x2,..., xn )Ą#
Ł# Ś#
T
Rozwiązanie tego układu równań oznacza określenie wektora x = [x1 x2 ... xn] ,
dla którego spełnione jest równanie (1.18).
Metoda Newtona-Raphsona
Metody rozwiązywania tego zagadnienia powstają przez odpowiednie
rozszerzenie metod rozwiązywania pojedynczych równań. W szczególności,
równanie (1.7) dla przypadku wielowymiarowego ma następującą postać
'
0 H" f (x0 ) + f (x0 )(� - x0 ) (1.19)
26
Rozwiązywanie równań nieliniowych
gdzie wektor � przedstawia współrzędne punktu, w którym spełniony jest warunek
(1.18).
'
Macierz określająca pochodną f (x0 ) jest nazywana Jakobianem (macierzą
Jakobiego)
"f1 "f1 "f1
Ą# ń#
ó#"x1 "x2 L Ą#
"xn
ó#"f "f2
"f2 Ą#
2
ó# Ą#
"f (x)
L
'
J( f (x)) = f (x) = = (1.20)
ó#"x1 "x2 "xn Ą#
"x
ó# Ą#
M M L M
ó#"fn "fn
"fn Ą#
ó#"x "x2 L Ą#
"xn Ś#
Ł# 1
Analogicznie do (1.8), rozwinięcie (1.24) prowadzi do następującej iteracyjnej
procedury rozwiązywania układu równań (1.18)
xk +1 = xk - J-1(f (xk ))f (xk ) (1.21)
jeśli det[J(f (xk ))]`" 0 , przy czym
J(f (xk ))= J( f (x))
x=xk
Algorytm (1.21) jest znany jako metoda Newtona-Raphsona iteracyjnego
rozwiązywania układu równań nieliniowych. W programach komputerowych wzór
(1.21) jest realizowany przez następujący algorytm
- oblicz f (xk ) ,
'
- oblicz J(f (xk ))= f (xk ) ,
- rozwiąż układ równań liniowych J(f (xk ))zk = f (xk )
- podstaw xk +1 = xk - zk
W charakterze oceny zbieżności procesu iteracyjnego można przyjąć normę
wektora zk odniesioną do normy wektora xk
zk
< � (1.22)
xk
Ze względu na ograniczoną dokładność obliczania funkcji f (xk ) oraz Jakobianu
J(f (xk )), dokładność całego algorytmu jest ograniczona. Objawia się to tym, że
począwszy od pewnej wartości minimalnej, norma wektora zk zacznie narastać. Jest
to sygnał, że należy skończyć obliczenia. Wynika stąd następujące kryterium
zakończenia obliczeń
zk +1 > � zk (1.23)
gdzie � jest rzędu jedności.
Rozwiązywanie równań nieliniowych 27
Metoda siecznych
Również metoda siecznych może być rozszerzona na przypadek
wielowymiarowy. A
atwo zauważyć, że równanie (1.10) można uogólnić następująco
-1
xk +1 = xk - "Xk("Fk ) f (xk ) (1.24)
przez analogię do rozwinięcia (1.19)
-1
0 H" f (xk ) + "Fk("Xk ) (� - xk ) (1.25)
gdzie
"Xk , "Fk są macierzami n � n o kolumnach, odpowiednio:
j k j
"xkj = x - xk oraz "f = f (x ) - f (xk ) , j = k - n, k - n - 1,..., k - 1.
j
Równania (1.22), (1.24) mają sens wtedy, gdy macierze "Xk , "Fk są nieosobliwe.
Jednakże zbieżność ciągu xk wymaga silnej nieosobliwości wszystkich macierzy
"Xk , co oznacza, że moduł wyznacznika tej macierzy powinien być dostatecznie
duży.
Z równania (1.22) widać, że w każdym kroku metody siecznych dla przypadku
wielowymiarowego wymagana jest znajomość n + 1 wartości wektora x oraz tyluż
wartości funkcji f (x) . Algorytm iteracyjny składa się z następujących kroków
- warunki początkowe: założyć wartości wektorów: x-n , x-n+1 , ..., x0
oraz przyjąć numer kroku iteracji k = 0
- obliczyć macierze "Xk , "Fk
- rozwiązać układ równań liniowych "Fkzk = f (xk )
- obliczyć nową wartość wektora xk +1 = xk - "Xkzk
Należy zauważyć, że ograniczenia warunkujące stosowanie metody siecznych
mogą uniemożliwiać wykonanie kolejnych kroków procesu iteracyjnego. Trzeba
zatem stosować odpowiednie rozwiązania (inne metody pomocnicze), pozwalające
uniknąć zatrzymania obliczeń.
4. Interpolacja
4.1. Wprowadzenie
Zadanie interpolacji odnosi się do działań zmierzających do przedstawienia
funkcji w postaci ciągłej, gdy znana jest ona w postaci dyskretnej. Jest to zatem
zdanie odwrotne do dyskretyzacji lub próbkowania wielkości ciągłej.
Załóżmy, że dla danego zbioru zmiennych niezależnych z przedziału < a;b > :
x1, x2,..., xn+1 znane są przyporządkowane im wartości funkcji: y1, y2,..., yn+1 .
Zależność ta jest zazwyczaj przedstawiana w postaci tabelarycznej:
y1 = f (x1) ,
y2 = f (x2 ) ,
...
yn+1 = f (xn+1) .
Zadaniem interpolacji jest wyznaczenie przybliżonych wartości funkcji dla wartości
zmiennych niezależnych z przedziału < a;b > , lecz nie będących punktami ze zbioru
x1, x2,..., xn+1. Jest to bardzo ogólne sformułowanie zadania i łatwo zauważyć, że
istnieje nieskończenie wiele sposobów jego rozwiązania, jeśli nie jest zadany sposób
przeprowadzenia funkcji interpolacyjnej przez zadane punkty (rys. 1.1).
y
x1 x2 x3 x4 x5 x6 x7 x8 x9 x
Rys. 1.1. Zasada interpolacji funkcji dyskretnej
Najczęściej poszukuje się funkcji interpolacyjnej o ściśle określonej postaci, tak,
aby zachowywała się ona w określony sposób. Są to często wielomiany algebraiczne
lub trygonometryczne.
Podstawowym celem interpolacji jest określenie wartości funkcji danej w postaci
stabelaryzowanej dla zmiennej x mieszczącej się pomiędzy danymi zawartymi w
tablicy. Można w ten sposób zapamiętać w komputerze zależność określoną na
30
Interpolacja
podstawie pomiaru. Typowymi przykładami zastosowania interpolacji jest
obliczanie całek oraz pochodnych funkcji dyskretnych w czasie. W takim przypadku,
w celu poprawienia dokładności obliczenia całki można skorzystać z wartości funkcji
aproksymującej określonej dla dowolnych wartości argumentu.
4.2. Wielomian interpolacyjny Newtona
Załóżmy, że dana jest funkcja f (x) w postaci tablicy, w której punktom
x1, x2,..., xn , zwanym węzłami interpolacji, przyporządkowane są wartości
f (x1), f (x2 ),..., f (xn ) . Zakłada się, że xi `" x dla i `" j . Funkcja interpolacyjna może
j
być określona w postaci wielomianu:
P(x) = b1 + b2x + b3x2 + ... + bn+1xn (1.1)
Jeśli funkcja dyskretna f (x) dana jest w dwóch punktach (n = 2), to funkcja
interpolacyjna w postaci (1.1) redukuje się do prostej (n+1 = 2). Podobnie, przez trzy
punkty (n = 3) można jednoznacznie poprowadzić parabolę, określoną przez
wielomian drugiego stopnia (n+1 = 3). Można dowieść, że w ogólnym przypadku,
dla n+1 punktów węzłowych (x , y ), istnieje tylko jeden wielomian P(x) spełniający
i i
warunek [1], [13]:
P(xi ) = yi , i=1, 2, ..., n+1 (1.2)
W przypadku wzoru interpolacyjnego Newtona, poszukiwany wielomian
interpolacyjny jest zapisywany w postaci:
P(x) = a1 + a2(x - x1) + a3(x - x1)(x - x2) + ... + an+1(x - x1)(x - x2)...(x - xn)
(1.3)
= b1 + b2x + b3x2 + L + bn+1xn
Korzystając z właściwości (1.2), powyższy zapis pozwala napisać następujący układ
równań:
P(x1) = y1 = a1
P(x2 ) = y2 = a1 + a2 (x2 - x1)
(1.4)
M
P(xn+1) = yn+1 = a1 + a2 (xn+1 - x1) + ... + an+1(xn+1 - x1)(xn+1 - x2 )...(xn+1 - xn )
co można zapisać w następującej postaci macierzowej:
"X �" a = y (1.5)
gdzie:
1 a1 y1
Ą# ń# Ą# ń# Ą# ń#
ó# Ą# ó# Ą# ó# Ą#
0 a2 y2
ó#1 x2 - x1 Ą# ó# Ą# ó# Ą#
"X = , a = , y =
ó#M M Ą# ó# Ą# ó# Ą#
M M M M M
ó# Ą# ó# Ą# ó# Ą#
ó#1 xn+1 - x1 xn+1 - x2 L xn+1 - xn+1Ą# ó#a Ą# ó#
yn+1Ą#
Ł# Ś# Ł# n+1Ś# Ł# Ś#
Interpolacja 31
Jak widać, macierz "X ma dogodną trójkątną postać, co pozwala bezpośrednio podać
rozwiązanie równania (1.5), wykonując podobne działania, jak w procedurze
odwrotnego podstawiania algorytmu Gaussa:
a1 = y1
y2 - a1
a2 =
x2 - x1
(1.6)
y3 - a1 - a2 (x3 - x1)
a3 =
(x3 - x1)(x3 - x2 )
...
Przykład. Interpolowana funkcja dana jest dla: x = 1, x = 2, x = 4, przy czym
1 2 3
wartości funkcji przyjmują wartości: f(x ) = 1, f(x ) = 4, f(x ) = 0. Określić
1 2 3
funkcję interpolacyjną.
Poszukiwana funkcja ma postać jak w (1.3), przy czym współczynniki są obliczane zgodnie z
(1.6):
a = f(x ) = 1,
1 1
f (x2) - a1 4 -1
a2 = = = 3 ,
x2 - x1 2 -1
f (x3) - a1 - a2(x3 - x1) 0 -1- 3(4 -1) - 5
a3 = = = .
(x3 - x1)(x3 - x2) (4 -1)(4 - 2) 3
Na podstawie (1.3), współczynniki funkcji interpolacyjnej (1.1) przyjmą następujące wartości:
b = a  a x + a x x =  16/3, b = a  a (x + x ) = 8, b = a = -5/3.
1 1 2 1 3 1 2 2 2 3 1 2 3 3
Zatem, funkcja interpolacyjna ma następującą postać:
-1
P(x) = (16 - 24x + 5x2).
3
Na rys. 1.2 pokazany jest przebieg uzyskanej funkcji interpolacyjnej z zaznaczonymi
wartościami danej funkcji dyskretnej.
6
4
P(x)
2
0
-2
-4
-6
-8
0 1 2 3 4 5
Rys. 1.2. Przebieg funkcji interpolacyjnej
32
Interpolacja
Powyższe zależności wygodnie jest zapisać wprowadzając pojęcie ilorazów
różnicowych. Oznaczmy i - tą różnicę
hi = xi+1 - xi (1.7)
Wyrażenia
f (xi+1) - f (xi ) f (xi+1) - f (xi )
"i = = (1.8)
hi xi+1 - xi
nazywa się ilorazami różnicowymi pierwszego rzędu. Odpowiednio także:
"i+1 - "i (xi+1 - xi )( f (xi+2 ) - f (xi+1))- (xi+2 - xi+1)( f (xi+1) - f (xi ))
"(i2) = = (1.9)
hi+1 + hi (xi+2 - xi+1)(xi+1 - xi )(xi+2 - xi )
- iloraz drugiego rzędu
"(i2) - "(i2) "(i2) - "(i2)
+1 +1
"(i3) = = - iloraz trzeciego rzędu (1.10)
hi+2 + hi+1 + hi xi+3 - xi
i ogólnie:
"(ik -1) - "(ik -1)
+1
"(ik ) = dla k = 1, 2,..., n - 1, i = 1, 2,..., n - k (1.11)
xi+k - xi
Dla zbioru n par (xi , f (xi )) można utworzyć tablicę ilorazów różnicowych, zwaną
schematem ilorazów różnicowych (patrz Tablica).
Wielomian interpolacyjny Newtona ma następującą postać:
P(x) = f (x1) + "1(x - x1)+ "(2)(x - x1)(x - x2)+ ... + "(n-1)(x - x1)(x - x2)...(x - xn-1) (1.12)
1 1
Widać, że w postaci wielomianu P(x) występują współczynniki z górnej przekątnej
schematu ilorazów różnicowych. Można sprawdzić, że
P(xi ) = f (xi ) dla i = 1, 2,..., n .
Algorytm interpolacyjny Newtona sprowadza się więc do obliczenia ilorazów różni-
cowych oraz określenia wartości wielomianu dla konkretnej wartości zmiennej x .
Ważne znaczenie ma przypadek, gdy wszystkie punkty stałe (węzły) są jednakowo
od siebie oddalone. Wówczas mamy: h = hi = xi+1 - xi = const oraz
Interpolacja 33
xi f (xi ) Ilorazy różnicowe
I rząd II rząd III rząd IV rząd V rząd
x1 f (x1)
"1
"(2)
x2 f (x2 )
1
"(3)
"2
1
"(4)
x3 f (x3) "(2) 1
2
"(5)
"3 "(3) 1
2
"(4)
x4 f (x4 ) "(2) 2
3
"4 "(3)
3
x5 f (x5)
"(2)
4
"5
x6 f (x6 )
f (xi+1) - f (xi ) "i+1 - "i f (xi+2 ) - 2 f (xi+1) + f (xi )
"i = , "(i2) = = ,
h 2h 2h2
"(ik -1) - "(ik -1)
+1
"(ik ) = dla k = 1, 2,..., n - 1, i = 1, 2,..., n - k
kh
dzięki czemu upraszcza się reprezentacja funkcji interpolacyjnej (1.12), gdyż:
x2 = x1 + h , xk = x1 + (k - 1)h , k = 1, 2,..., n (1.13)
Wprowadzając nową zmienną
q = x - x1 ,
a zatem: q - h = x - x2 = x - x1 - h , q - (n - 2)h = x - xn-1 = x - x1 - (n - 2)h
otrzymamy z (4.6)
P(q) = f (x1) + q"1 + q(q - h)"(2) + ... + q(q - h)...(q - (n - 2)h)"(n-1) (1.14)
1 1
Nie ma więc potrzeby obliczania współczynników wielomianu (1.3).
Przykład. Określić ogólną postać wielomianu interpolacyjnego Newtona dla funkcji
dyskretnej reprezentowanej przez trzy kolejne równooddalone punkty.
W tym przypadku wielomian interpolacyjny (1.14) jest ograniczony do trzech wyrazów:
P(q) = f (x1) + q"1 + q(q - h)"(2) ,
1
co, po podstawieniu odpowiednich wartości, przyjmuje następującą postać:
f (x2 ) - f (x1) f (x3) - 2 f (x2 ) + f (x1)
P(q) = f (x1) + q + q(q - h)
h 2h2
Podstawiając d = q / h ( d jest w ten sposób względną odległością od początku przedziału),
otrzymamy:
34
Interpolacja
f (x3) - 2 f (x2 ) + f (x1)
P(d ) = f (x1) + d( f (x2 ) - f (x1))+ d(d -1)
2
Po uporządkowaniu otrzymamy:
1
2
P(d ) = (2 f (x1) - d(3 f (x1) - 4 f (x2 ) + f (x3))+ d ( f (x1) - 2 f (x2 ) + f (x3)))
2
W ten sposób, na przykład, wartość funkcji w środku drugiego ( d =1,5 ) odcinka można
oszacować następująco:
1 3 9
�#2
P(1.5) = f (x1) - (3 f (x1) - 4 f (x2 ) + f (x3))+ ( f (x1) - 2 f (x2 ) + f (x3))ś#
ś# ź#
2 2 4
�# #
1
= (3 f (x3) + 6 f (x2 ) - f (x1))
8
4.3. Numeryczne różniczkowanie funkcji dyskretnej
Funkcję interpolującą można wykorzystać do określenia algorytmu numerycznego
różniczkowania funkcji dyskretnej. Odpowiednią formułę uzyskuje się przez
różniczkowanie funkcji interpolującej:
d
D(x) = P(x)
dx
Na przykład, dla aproksymacji 3-punktowej (jak w Przykładzie 4.1), otrzymamy:
d 1 f (x3) - 2 f (x2 ) + f (x1) f (x3) - 2 f (x2 ) + f (x1)
�# ś#
P(d ) = f (x2 ) - f (x1) - + 2d
ś# ź#
dd h 2 2h
�# #
- f (x3) + 4 f (x2 ) - 3 f (x1) d( f (x3) - 2 f (x2 ) + f (x1))
= +
2h h
Dla różniczkowania na końcu przedziału ( d = 2 ) otrzymamy:
3 f (x3) - 4 f (x2 ) + f (x1)
D(2) =
2h
f (x3) - f (x1)
Podobnie, w środku przedziału: D(1) = .
2h
5. Aproksymacja
5.1. Wprowadzenie
Zadanie aproksymacji polega na przybliżeniu funkcji Y (x) za pomocą innej funkcji
f (x) , która odnosi się do tego samego obszaru. W praktycznych zastosowaniach
inżynierskich Y (x) jest najczęściej funkcją dyskretną Y (x) = Y (xi ) = yi ,
i = 0,1,...,M -1, reprezentującą dane pomiarowe znane dla M różnych wartości
zmiennej niezależnej, a f (x) przedstawia model (wzorzec) obserwowanego procesu.
Zależność Y (x) jest nazywana funkcją aproksymowaną, natomiast f (x) jest funkcją
aproksymującą. Zakłada się, że dostępne próbki yi obarczone są błędami, zatem
aproksymacja ma na celu najlepsze przybliżenie danych za pomocą funkcji
aproksymującej zgodnie z przyjętym kryterium.
Funkcję f (x) wygodnie jest przedstawić w następującej postaci:
f (x) = a0�0(x) + a1�1(x) + ... + aN -1�N -1(x) (1.1)
gdzie �i (x) , i = 0,1,..., N -1 są funkcjami bazowymi, natomiast ai , i = 0,1,..., N -1
przedstawiają poszukiwane parametry funkcji. Można zauważyć, że funkcja (1) jest
liniowa względem nieznanych parametrów
Problem ten ilustruje rys. 1.1.
1,5
funkcja aproksymująca f(x)
1
funkcja aproksymowana Y(x)
0,5
0
-0,5
-1
0 5 10 x
Rys. 1.1. Ilustracja zasady aproksymacji funkcji
Parametry funkcji aproksymującej f (x) są określane na podstawie przyjętego
kryterium. Na przykład, żąda się, aby spełnione było minimum normy różnicy obu
36
Aproksymacja
funkcji: min(Y (x) - f (x) ). W przypadku funkcji dyskretnej warunek ten jest
równoważny kryterium najmniejszych kwadratów:
M -1
2
S2 = Y (x) - f (x) = (xi) - f (xi)) (1.2)
"(Y
i=0
Inne kryterium aproksymacji jest sformułowane za pomocą zależności:
S1 = sup f (x) - Y (x) (1.3)
x")#a,b*#
co oznacza, że poszukiwana funkcja f (x) powinna dać najmniejsze maksimum
różnicy pomiędzy daną funkcją i jej aproksymacją. Jest ono znane jako kryterium
aproksymacji jednostajnej.
Najszersze zastosowanie praktyczne znalazła aproksymacja według metody
najmniejszych kwadratów (MNK). Jest to związane z istnieniem bardzo efektywnych
obliczeniowo algorytmów, które wywodzą się z kryterium minimalizacji funkcji
(1.2).
5.2. Aproksymacja średniokwadratowa
Załóżmy, że znana jest funkcja Y (xi ) = yi na zbiorze dyskretnym xi ,i = 0,1,..M -1 w
przedziale < a, b > . Chcemy, aby wartości funkcji yi były przybliżone przez funkcję
f (x) o postaci (1.1) w tym samym przedziale.
Jeśli odwołać się do interpretacji pomiarowej, to yi , i = 0,1,..M -1 przedstawia
zbiór M pomiarów, w stosunku do których zakładamy, że spełnione są następujące
zależności:
y0 = a0�0(x0 ) + a1�1(x0 ) L aN -1�N -1(x0 ) + v0
y1 = a0�0(x1) + a1�1(x1) L aN -1�N -1(x1) + v1
(1.4)
M M M M M M
yM -1 = a0�0(xM -1) + a1�1(xM -1) L aN -1�N -1(xM -1) + vM -1
gdzie wielkości vi , i = 0,1,.., M -1 przedstawiają odchyłki (błędy) pomiędzy
postulowaną wartością funkcji aproksymującej f (xi ) , a dyskretnymi wartościami
funkcji aproksymowanej yi (zmierzonymi wartościami), N jest liczbą składników
funkcji aproksymującej (a zatem, liczbą nieznanych współczynników a , i=0,1,..,N-1).
i
Dalsze rozważania wygodnie jest prowadzić, korzystając z zapisu macierzowego.
Układ równań (1.4) przyjmie następującą formę:
y0 �0(x0 ) �1(x0 ) L �N -1(x0 ) a0 v0
Ą# ń# Ą# ń#Ą# ń# Ą# ń#
ó# Ą# ó# Ą#ó#
y1 �0(x1) �1(x1) L � (x1) a1 Ą# ó# v1 Ą#
N -1
ó# Ą# ó# Ą#ó# Ą# ó# Ą#
= + (1.5)
ó# M Ą# ó# M M M M Ą#ó# M Ą# ó# M Ą#
ó# Ą# ó#� (xM ) �1(xM ) L �N (xM )Ą#ó#aN Ą# ó#v Ą#
yM
Ł# -1Ś# Ł# 0 -1 -1 -1 -1 Ś#Ł# -1Ś# Ł# M -1Ś#
co w ogólnym zapisie wygląda następująco:
Aproksymacja 37
y = Ha + v (1.6)
gdzie:
y = [ y0 y1K yM -1]T - wektor reprezentujący aproksymowaną funkcję dyskretną,
H =[h(0) h(1) K h( M -1)]T - macierz modelu określonego przez funkcje
bazowe:
h(i) = [�0(xi ) �1(xi ) L �N -1(xi )] , i = 0,1,..M -1 ,
v = [v0 v1KvM -1]T - wektor błędów pomiarowych,
a = [a0 a1 L aN -1]T - wektor estymowanych parametrów.
Można zauważyć, że odchyłki pomiędzy danymi wartościami aproksymowanej
funkcji Y (xi ) = yi , a wartościami postulowanej funkcji aproksymującej f (xi ) można
określić następująco:
vi = F(xi )-f (xi) = yi - h(i) �" a, i = 0,..,M - 1 (1.7)
Funkcja kryterialna S2 = S2(a) (1.2) jest określona, jako suma kwadratów odchyłek
(błędów) (1.7), co można zapisać w następującej postaci:
M -1 M -1
2
2
S2(a) = = )-f (xi)) = vT v , (1.8)
"vi "(F(xi
i=0 i=0
przy czym, na podstawie (6):
v = y - Ha
Funkcja (1.8) osiąga minimum, gdy jej pochodne względem parametrów,
określonych przez wektor współczynników a , przyjmują zerowe wartości:
M -1
"S2(a) "
2
= )-f (xi)) = 0
"(F(xi
"a0 "a0 i=0
M -1
"S2(a) "
2
= )-f (xi)) = 0
"(F(xi
"a1 "a1 i=0 (1.9)
L
M -1
"S2(a) "
2
= )-f (xi)) = 0
"(F(xi
"aN -1 "aN -1 i=0
co można zapisać w postaci wektorowej:
"S2(a) "
T
= ((y - Ha) (y - Ha))= 2HTHa - 2HTy = 0 (1.10)
"a "a
Korzysta się tu z zależności:
(y - Ha)T (y - Ha)=(yT - aT HT )(y - Ha)= yT y - yT Ha - aT HT y + aT HT Ha
(1.11)
= yT y - 2aT HT y + aT HT Ha
38
Aproksymacja
Ostatnia równość wynika z faktu, że wielkość S2 (a) jest skalarem, a więc także:
yT Ha = aT HT y = p ,
a więc:
" "
(yT Ha - aT HT y)= 2 (aT HT y)= 2HT y
"a "a
p  wielkość skalarna.
Podobnie, różniczkując ostatni składnik w (1.11), otrzymamy:
"
(aT HT Ha)= HT Ha + aT HT H = 2HT Ha .
"a
Zatem, z (1.10) uzyskuje się następującą równość:
HT Ha = HT y (1.12)
Zauważmy, że macierz HT H jest kwadratowa o wymiarze N , wyrażenie z prawej
strony (1.12) jest wektorem o długości N , a N -elementowy wektor a zawiera
poszukiwane współczynniki aproksymującej funkcji (1.1). Równanie (1.12)
przedstawia zatem klasyczny liniowy układ równań z N niewiadomymi. Można go
rozwiązać jedną ze znanych metod.
Równanie w postaci (1.12) jest nazywane równaniem normalnym. Formalnie, dla
warunku M e" N , jego rozwiązanie można zapisać w postaci:
-1
a = (HTH) HTy (1.13)
-1
gdzie: macierz prostokątna H+ = (HTH) HT jest nazywana macierzą pseudo-
odwrotną (macierzą Moore a-Penrose a). W przypadku, gdy M = N , macierz
pseudo-odwrotna H+ w (1.13) odpowiada macierzy odwrotnej H-1 (jeśli macierz H
jest nieosobliwa).
Niektóre właściwości macierzy pseudo-odwrotnej:
T T
H+H = (H+H) , HH+ = (HH+) ,
H+HH+ = H+ , HH+H = H .
Przykład
Dane są cztery punkty na płaszczyz
nie (x,y): (-2, 20), (1, 2), (2, 7), (3, 12)  rys. 5.2. Określić
parabolę, która najlepiej, w sensie kryterium najmniejszych kwadratów, aproksymuje
podaną funkcję dyskretną.
Funkcja aproksymująca jest określona w postaci wielomianu drugiego stopnia:
f (x) = ax2 + bx + c
Podstawiając kolejne punkty do powyższego równania, otrzymamy następujący
nadokreślony układ równań:
4a - 2b + c = 20
a + 2b + c = 2
4a + 2b + c = 7
9a + 3b + c = 12
którego postać macierzowa jest następująca:
Aproksymacja 39
4
Ą# - 2 1 20
ń# Ą# ń#
Ą#
ó#1 1 1Ą#ó#ań# ó# Ą#
2
ó# Ą# ó# Ą#
bĄ# =
ó# Ą#ó# Ą# ó# Ą#
4 2 1 7
ó# Ą#ó#cĄ# ó# Ą#
Ś#
Ł#9 3 1Ś#Ł# Ł#12Ś#
Stosując zależność (1.13) otrzymamy:
-1
�# 4 20
Ą# - 2 1 ś#
ń#ź# Ą# ń#
ś#
a 4 1 4 9 4 1 4 9
Ą# ń# Ą# ń#ó# Ą# ń#ó# Ą#
1 1 1Ą#ź# ó# 2
ś#
ó#bĄ# ó#
Ą# Ą#
=
ś#
ó# Ą# ó#- 2 1 2 3Ą#ó# 2 1Ą#ź# ó#- 2 1 2 3Ą#ó# 7 Ą#
Ą#ó#4 Ą#ó#
ś# ź#
ó#cĄ# ó# 1 1 1 1Ą#ó# ó# 1 1 1 1Ą#ó# Ą#
Ł# Ś# Ł# Ś# Ł# Ś#
ś#
Ł#9 3 1Ą#ź# Ł#12Ś#
Ś#
�# #
20
Ą# ń#
92 -196 - 68 172 759
Ą# ń#ó# Ą# Ą# ń# 2,097
Ą# ń#
1
ó# Ą#ó# 2 Ą# 1 ó# Ą# ó#
= = 3,569Ą#
ó#- Ą#
ó#- 376 140 152 84 Ą#ó# 7 Ą# = 362 ó#-1292Ą#
1448
ó# Ą#
4,384
ó# Ą#
ó# 324 1104 516 - 496Ą#ó# Ą# Ł# 1587
Ł# Ś#
Ł# Ś# Ś#
Ł#12Ś#
Przebieg uzyskanej funkcji aproksymującej jest pokazany na rys. 1.2.
y
25
20
f(x)
15
10
x
-3 -2 -1 0 1 2 3 4
Rys. 1.2. Aproksymacja funkcji za pomocą paraboli
W ogólnym przypadku, niektóre pomiary reprezentowane przez aproksymowaną
funkcję mogą być bardziej wiarygodne od innych. Wówczas ich wpływ na przebieg
obliczanej funkcji aproksymującej powinien być większy. Można to uwzględnić
przez wprowadzenie współczynników wagowych wi do funkcji kryterialnej (1.8):
M -1 M -1
2
S2(a) = vi2 = (1.14)
i i
"w "w(x )(F(xi )-f (xi )) = vT Wv ,
i =0 i =0
gdzie: W jest wagową macierzą kwadratową diagonalną o wymiarach M � M ; na jej
przekątnej leżą współczynniki w(xi ) , których wartości są zwykle normalizowane:
40
Aproksymacja
0 d" w(xi) d" 1. Przy braku informacji o wspomnianej wiarygodności pomiarów,
współczynniki wagowe przyjmują wartość 1 ( W jest wówczas macierzą
jednostkową).
Uwzględnienie macierzy wagowej w (1.14) prowadzi do następującej postaci
równania (1.12):
HTWHa = HTWy (1.15)
oraz, odpowiednio:
-1
a = (HTWH) HTWy (1.16)
Metoda najmniejszych kwadratów, zgodnie z którą współczynniki funkcji
aproksymującej są określane według (1.13) lub (1.16), jest bardzo użytecznym i
praktycznym narzędziem w wielu zastosowaniach. Niektóre z nich są rozpatrywane
poniżej.
5.3. Filtr wygładzający
Problem wygładzania danych pomiarowych polega na przybliżeniu
obserwowanych (mierzonych) parametrów procesu za pomocą przyjętych
zależności, które stanowią model tego procesu. Odchylenia od modelu są traktowane
jako zakłócenia. Aby uzyskać poprawne wygładzanie danych pomiarowych,
przyjęty model powinien być zbliżony do przebiegu obserwowanego procesu
(funkcja aproksymująca powinna mieć dostatecznie dużo stopni swobody), a
jednocześnie nie powinien on być ściśle dopasowany do danych empirycznych
(zakłócenia nie powinny być odzwierciedlone w modelu).
Załóżmy, że funkcja Y (x) przedstawia proces, który jest obserwowany poprzez
próbkowanie określonego parametru i dostępnych jest M ostatnich jego próbek:
Y (xi ) = yi , i = 0,1,..., M -1. Proces ten jest reprezentowany (modelowany) za pomocą
funkcji aproksymującej w postaci wielomianu stopnia N -1 < M :
f (xi ) = a0 + a1xi + ... + aN -1xiN -1 (1.17)
Dla uproszczenia załóżmy, że próbkowanie odbywa się ze stałym okresem T .
Ponieważ dostępne są próbki w oknie pomiarowym o szerokości M , więc
aproksymacja jest równoważna filtracji nierekursywnej w tym właśnie oknie.
Odpowiedni filtr jest określony następującym równaniem:
M -1
gx (k) = h(i)yk -M +i+1 = hy(k) (1.18)
"
m
i=0
gdzie xm oznacza wartość zmiennej niezależnej względem której określona jest k -ta
próbka odpowiedzi filtru, 0 d" xm d" (M -1)T ; h = [h(0) h(1) ... h(M -1)] - wektor
współczynników filtru; y(k) = [yk -M +1 yk -M +2 ... yk ]T - wektor zawierający
ostatnie M próbek sygnału wejściowego.
Aproksymacja 41
Współczynniki filtru (1.18) należy tak dobrać, aby funkcja gx (k) aproksymowała
m
przebieg dyskretny określony przez wektor y(k) w punkcie xm , licząc od początku
przedziału, zgodnie z modelem f (x) (1.17) według kryterium najmniejszych
kwadratów.
Zgodnie z przedstawionym opisem funkcja f (x) aproksymuje mierzony
dyskretny przebieg według następującej relacji:
f (xi ) = yi + vi - w punktach odpowiadających kolejnym próbkom. (1.19)
xi =iT
Zakłada się zatem, że funkcja będzie aproksymowana tylko w węzłach
odpowiadających punktom próbkowania.
Stosowne zależności dla metody najmniejszych kwadratów można napisać przy
założeniu, że punkt odpowiadający zmiennej xm znajduje się w początku układu
współrzędnych ( xm = 0 ). Wówczas dla jednego zbioru M próbek otrzymamy:
a0 + a1(-xm) + a2(-xm)2 + ... + aN -1(-xm)N -1 = y0
a0 + a1(-xm + T ) + a2(-xm + T )2 + ... + aN -1(-xm + T )N -1 = y1
...
a0 + a1(-xm + (m -1)T ) + a2(-xm + (m -1)T )2 + ... + aN -1(-xm + (m -1)T )N -1 = ym-1 (1.20)
a0 + a1(0) + a2(0)2 + ... + aN -1(0)N -1 = ym
...
a0 + a1((M -1)T - xm)+ a2((M -1)T - xm)2 + ...+ aN -1((M -1)T - xm)N -1 = yM -1
W punkcie odniesienia równanie modelu ma zatem następującą postać: a0 = ym .
Jest to wynikiem odpowiedniego przesunięcia osi czasu, jednak takie założenie jest
pomocne dla uproszczenia kolejnych kroków procedury syntezy filtru.
Równania (1.20) można zapisać w postaci macierzowej:
Aa = y (1.21)
gdzie, jeśli dla uproszczenia przyjąć T = 1:
Ą#1 - m y0
Ą# ń#
(- m)2 ... (- m)N -1 ń# Ą# a0 ń#
ó# Ą#
ó# Ą#
-1
y1
ó# Ą#
ó#1 (1- m) (1- m)2 ... (1- m)N Ą#
ó# Ą#
a1
ó# Ą#
ó#.. ... Ą#
ó# Ą#
M
... ... ...
ó# Ą#
A = ó# Ą# , a = a2 , y =
ó# Ą#
ym
1 0 0 ... 0 ó# Ą#
ó# Ą#
ó# Ą#
M
ó# Ą#
ó#.. ... Ą#
ó# Ą#
... ... ...
ó#aN -1Ś# ó# M Ą#
Ą#
ó# Ą#
Ł#
ó# - m -1) (N - m -1)2 ... (N - m -1)N -1Ś#
Ą# ó# Ą#
Ł#yM -1Ś#
Ł#1 (N
Zgodnie z przyjętymi założeniami, aproksymacja odbywa się w odniesieniu do m -tej
próbki w oknie pomiarowym ( m =1...M ), co oznacza, że z lewej strony tej próbki
znajduje się nL = m -1 próbek, a z prawej: nP = M - m -1 próbek. Wielkość nP
określa opóz
nienie odpowiedzi algorytmu na sygnał wejściowy i jest nazywane
opóz
nieniem grupowym [12].
Wektor poszukiwanych współczynników wielomianu jest określony następująco:
42
Aproksymacja
-1
a = (AT A) AT y . (1.22)
Macierz A nie zależy od pomiarów, zatem część powyższego równania może być
określona przed rozpoczęciem obliczeń. A
atwo zauważyć, że [11]:
nP nP
i+ j
{aij}= {AT A}= Aki Akj = (1.23)
" "k , i, j = 0,1,..., N -1.
k =nL k =nL
Wracając do problemu syntezy filtru (1.18) można zauważyć, że sygnał wyjściowy
gx (k) = gm (k) jest estymatą próbki y(m) H" a0 (co wynika ze struktury równania
m
(1.6)). Zatem, równanie (1.22) przedstawia filtr (1.18), jeśli obliczać w nim tylko
współczynnik a0 . Kolejne współczynniki filtru są utworzone przez pierwszy wiersz
-1
macierzy (AT A) AT . Można je określić zgodnie z następującymi zależnościami:
-1
-1 -1
h(0) = (AT A) AT v(0) , h(1) = (AT A) AT v(1) , h(M -1) =(AT A) AT v(M -1) (1.24)
T T
gdzie: v(0) = [1 0 L 0] , v(1) = [0 1 L 0] , v(M -1) = [0 L 0 1]T .
Obliczanie współczynników filtru (1.24) można uprościć, jeśli zauważyć, że [11]:
N -1
-1 -1
ż# ż# �#
h(l) = (AT A) AT v(l)�# = (AT A) li (1.25)
"
�# Ź# �# Ź#
�# �#0 i=0 �# �#0i
Utworzony w ten sposób filtr nierekursywny (1.18) aproksymuje zbiór M kolejnych
danych pomiarowych, uzyskanych w równych odstępach czasu, według
aproksymującej funkcji wykładniczej (1.17).
5.4. Filtr różniczkujący
Filtr służący do określania pierwszej pochodnej funkcji dyskretnej może być łatwo
utworzony na podstawie przedstawionego powyżej algorytmu. Można zauważyć, że
pochodna funkcji aproksymującej (1.17) ma następującą postać:
df (xi )
= a1 + 2a2xi + ... + (N -1)aN -1xiN -2 (1.26)
dx
Współczynniki poszukiwanego filtru różniczkującego:
M -1
dxm (k) = d(i)yk -M +i+1 = dy(k) (1.27)
"
i=0
są zatem określone przez zbiór współczynników a1 funkcji (1.26). Można je obliczyć
-1
zgodnie z (1.25), przy czym, należy wziąć drugi wiersz macierzy (AT A) AT
(oznaczony indeksem 1):
N -1
-1 -1
ż# ż# �#
d(l) = (AT A) AT v(l)�# = (AT A) li (1.28)
"
�# Ź# �# Ź#
�# �#1 i=0 �# �#1i
Aproksymacja 43
W celu uzyskania większej dokładności filtracji (mniejsza wariancja wyników), w
charakterze punktu odniesienia należy brać punkt, leżący najbliżej środka okna
pomiarowego.
Należy zauważyć, że zarówno filtr wygładzający, jak i różniczkujący, pozwalają
aproksymować wejściową funkcję dyskretną (lub jej pochodną) tylko w punktach,
odpowiadających momentom próbkowania. Jeśli funkcja aproksymująca powinna
być określona dla dowolnych wartości argumentu, to należy stosować postać (1.17), a
w związku z tym, powinny być wyznaczone wszystkie współczynniki funkcji ai ,
i = 0,1,..., N -1.
5.5. Przykład obliczeniowy
Rozważmy problem wygładzania danych pomiarowych i określania pochodnej
funkcji reprezentowanej tymi danymi za pomocą odpowiednich filtrów
rekursywnych. Zarejestrowany przebieg jest przedstawiony na rys. 5.3a, krzywa 1.
Próbkowanie odbywa się ze stałą częstotliwością.
Do wygładzania tego przebiegu zastosowano filtr (1.18), który został
zaprojektowany na podstawie funkcji aproksymacyjnej (1.17) 2-go rzędu. Założono,
że w oknie przetwarzania filtru znajduje się M = 9 próbek. Filtracja jest prowadzona
w odniesieniu do środkowej próbki w oknie pomiarowym, a więc m = 5 .
Po zastosowaniu przedstawionej powyżej procedury uzyskuje się następującą
funkcję impulsową filtru wygładzającego i różniczkującego:
15
a)
3
10
1
5
2
0
0 10 20 30 40 xk
1
b)
0.5
0
-0.50 10 20 30 40 xk
Rys. 1.3.
1
h = [- 21 14 39 54 59 54 39 14 - 21]
231
1
d = [- 4 - 3 - 2 -1 0 1 2 3 4]
60
W wyniku filtracji (1.18) oraz (1.27) otrzymuje się przebiegi wyjściowe (rys. 5.3):
wygładzonych danych (krzywa 2) oraz estymaty pochodnej (rys. 5.3b). W celu
lepszego porównania przebiegu oryginalnego z wygładzonym (aproksymowanym),
44
Aproksymacja
ten ostatni został przesunięty o liczbę próbek stanowiących opóz
nienie grupowe (4
próbki, krzywa 3).
Na zakończenie kilka uwag praktycznych dotyczących problemu aproksymacji.
- Przy wyborze funkcji bazowych należy się kierować podobieństwem
pomiędzy aproksymowaną funkcją i jej reprezentacją  wygładzoną ; w
charakterze funkcji bazowych najczęściej wybiera się szereg potęgowy lub
trygonometryczny.
- Liczba wyrazów w funkcji aproksymacyjnej decyduje o dokładności
odwzorowania oryginalnego przebiegu; niski rząd tej funkcji powoduje
zgrubne przybliżenie, dzięki czemu efekt filtracji jest bardziej wyrazisty.
- Kwestia ta łączy się również z szerokością okna przetwarzania (liczba M
jednocześnie branych pod uwagę próbek funkcji oryginalnej): im szersze jest
okno pomiarowe, tym lepsze jest wygładzanie danych. Przy szerokim oknie
pomiarowym można także stosować funkcje aproksymujące wyższego rzędu z
zachowaniem wierności odtwarzania. Niestety, zwiększenie długości okna
pomiarowego prowadzi do zwiększenia opóz
nienia grupowego, co ma istotne
znaczenie wówczas, gdy aproksymacja odbywa się bezpośrednio w czasie
pomiarów. Związana z tym zwłoka czasowa sprawia, że informacja o stanie
nadzorowanego procesu jest dostępna z określonym opóz
nieniem.
W przypadku filtracji sygnałów, model użyty do projektowania filtrów wygodnie
jest tworzyć na bazie funkcji trygonometrycznych sinus/kosinus. Można wówczas
łatwo uzyskać procedury do pomiaru amplitudy i fazy mierzonych sygnałów lub ich
harmonicznych [12].
5.6. Metoda Najmniejszych Kwadratów z wykorzystaniem rozkładu
macierzy według wartości szczególnych  SVD
Dowolna macierz A ( m � n ) może być przedstawiona w następującej postaci:
A = UWVT (1.29)
gdzie:
U - macierz ( m � m ), której kolumny spełniają następującą zależność:
m
uij = 1, 1d"kd"n, 1d"jd"n, to znaczy, są wzajemnie ortogonalne;
"uik
i=1
Ł 0
Ą# ń#
W =
ó#0 0Ą# - macierz ( m � n ) wartości szczególnych, przy czym:
Ł# Ś#
�1
Ą# ń#
ó# Ą#
�
2
ó# Ą# *
Ł = , �1 e" � e" L e" � e" 0 , r d" n  rząd macierzy A ; (1.30)
2 r
ó# L Ą#
ó# Ą#
�
Ł# r Ś#
V - macierz kwadratowa ( n � n ), której kolumny spełniają następującą zależność:
*
Rząd macierzy r jest największą liczbą niezależnych wierszy (lub kolumn) macierzy.
Aproksymacja 45
n
vij = 1, 1d"kd"n, 1d"jd"n, a więc są również wzajemnie ortogonalne.
"vik
i=1
Z powyższego opisu wynika, że:
UT U = VT V = VVT = 1 (1.31)
Kolumny macierzy U są wektorami własnymi macierzy kwadratowej AAT ,
natomiast kolumny macierzy V są wektorami własnymi macierzy AT A . Wartości
szczególne rozkładu (1.29)  a więc elementy diagonalne macierzy Ł - są natomiast
pierwiastkami kwadratowymi wartości własnych macierzy AAT lub AT A .
Zależność (1.29) jest nazywana rozkładem macierzy A według wartości
szczególnych (ang. SVD  Singular Value Decomposition). Wynika stąd sposób
obliczania macierzy rozkładu (1.29).
Wektor własny x macierzy H spełnia następujące równanie:
Hx = x (1.32)
przy czym  jest wartością własną macierzy H . W takim przypadku mówi się, że
wektor x jest skojarzony w wartością własną .
W celu określenia wartości własnych  oraz odpowiadających im wektorów
własnych x należy rozwiązać równanie (1.32), co jest równoważne następującej
zależności:
(H - 1)x = 0 (1.33)
Jednoznaczne rozwiązanie (1.33) uzyskuje się wówczas, gdy spełniony jest warunek:
h11 -  h12 L h1m
h21 h22 -  L h2m
det(H - 1) = = 0 (1.34)
M M L M
hm1 hm2 L hwmm - 
Rozwiązanie równania (1.34) jest równoważne znalezieniu pierwiastków
wielomianu m-tego stopnia względem  . W ogólnym przypadku jest zatem m
wartości własnych: 1,2 ,..m . Podstawiając kolejno te wartości własne do (1.33)
otrzymuje się m równań:
h11
Ą# - i h12 L h1m x1i
ń#Ą# ń#
ó#
h21 h22 - i L h2m Ą#ó# x2i Ą#
ó# Ą#ó# Ą#
= 0 , i = 1,2,..,m , (1.35)
ó# Ą#ó# Ą#
M M L M M
ó# Ą#ó# Ą#
hm1 hm2 L hmm - i Ś#Ł#xmi Ś#
Ł#
T
po rozwiązaniu których uzyskuje się m wektorów własnych: xi = [x1i x2i L xmi ] ,
i = 1,2,..,m .
W rozpatrywanym przypadku macierz H = AAT (jeśli obliczana jest macierz U )
lub H = AT A (jeśli obliczana jest macierz V ). Macierz Ł (1.30) powstaje przez
uporządkowanie pierwiastków z wartości własnych macierzy H . W standardowych
programach do rozkładu macierzy według wartości szczególnych stosowane są
46
Aproksymacja
zazwyczaj inne, bardziej efektywne algorytmy w stosunku do przedstawionego
powyżej.
Właściwości rozkładu SVD:
A = UWVT , AT = VWT UT (1.36)
AT A = VWT UT UWVT = VWT WVT , AAT = UWWT UT (1.37)
-1
A# = (AT A) AT = VW-1UT - macierz pseudo-odwrotna (1.38)
Ą# ń#
Ł-1 0
gdzie W-1 = W# = ; elementy diagonalne macierzy Ł-1 można obliczyć
ó# Ą#
0 0Ś#(n�m)
Ł#
jako odwrotności odpowiednich elementów diagonalnych macierzy Ł .
Ostatnią zależność można bezpośrednio wykorzystać do rozwiązywania zadań
MNK:
x = H#y = VW#UT y (1.39)
gdzie: y - wektor pomiarów, x - wektor poszukiwanych parametrów funkcji
aproksymującej.
Zastosowanie z użyciem programu MATLAB:
[U,W,V]=svd(A); rozkład macierzy A według wartości szczególnych,
[C,D]=eig(A*A'); wartości własne macierzy A*A'.
6. Całkowanie numeryczne
6.1. Wprowadzenie
Załóżmy, że należy obliczyć przybliżoną wartość całki oznaczonej:
b
I ( f ) = f (x)dx (1.1)
+"
a
na odcinku < a, b > .
W celu wyprowadzenia odpowiedniej formuły obliczeniowej przyjmijmy, że odcinek
< a, b > jest podzielony na n przedziałów elementarnych < xi , xi+1 > , i = 1, 2,..., n ,
przy czym x1 = a , xn+1 = b oraz x1 < x2 < ... < xn+1. Długość i -tego odcinka
oznaczymy przez hi :
hi = xi+1 - xi (1.2)
Wartość całki (1.1) można przedstawić w postaci sumy całek na elementarnych
odcinkach:
n
I( f ) = (1.3)
"I
i
i=1
przy czym
xi+1
Ii = Ii ( f ) = f (x)dx
+"
xi
Numeryczne przybliżanie wartości całek na danym odcinku nazywa się
kwadraturą.
6.2. Metoda Simpsona
Najprostsza kwadratura powstaje przez przyjęcie, że całka na elementarnym
odcinku < xi, xi+1 > jest równa polu trójkąta wyznaczonego przez boki: hi oraz
wartości funkcji w środku przedziału:
xi + xi+1
yi =
2
Zatem:
xi + xi+1
�# ś#
Ii H" hi f ( yi ) = hi f (1.4)
ś# ź#
2
�# #
48
Całkowanie numeryczne
lub
hi
�# ś#
Ii H" hi f xi + (1.5)
ś# ź#
2
�# #
Jeśli funkcja podcałkowa dana jest tylko w węzłach xi , to wówczas należy
stosować formułę:
IPi H" hi f (xi ) (1.6)
Zwiększenie dokładności obliczeń można uzyskać przez prostą operację zamiany
prostokąta do obliczania pola pod krzywą na trapez. Uzyskuje się wówczas
następującą zależność:
f (xi ) + f (xi+1)
ITi H" hi (1.7)
2
Metodą kombinacji liniowej obu metod otrzymuje się znacznie bardziej dokładną
metodę Simpsona:
2 1 1 �# xi + xi+1 ś#
�# ś#
ISi = IPi + ITi = hi ś# f (xi ) + 4 f + f (xi+1)ź# (1.8)
ś# ź#
3 3 3 2
�# #
�# #
Dla przypadku, gdy funkcja jest dostępna tylko w węzłach xi powyższa formuła
przyjmuje następującą postać:
1
ISi = hi( f (xi ) + 4 f (xi+1) + f (xi+2 )) (1.9)
3
7. Numeryczne rozwiązywanie równań różniczkowych
zwyczajnych
7.1. Wprowadzenie
Numeryczne metody rozwiązywania równań różniczkowych są stosowane w takich
przypadkach, gdy rozwiązania w postaci analitycznej nie są znane. Ponieważ jednak
tylko dla nielicznych równań można podać analityczne rozwiązanie, więc metody
numeryczne są w takich przypadkach bardzo pomocnym narzędziem poszukiwania
rozwiązania. Ważnym obszarem zastosowania tych metod są numeryczne modele
układów i zjawisk dynamicznych, które służą do ich komputerowej symulacji.
W odróżnieniu jednak od metod analitycznych, w metodach numerycznych
zakłada się, że zmienna niezależna dostępna jest tylko w wybranych, dyskretnych
wartościach. Konsekwencją tego jest nieuchronne przybliżenie rozwiązania.
Rozpatrzmy równanie różniczkowe pierwszego rzędu:
dy(t)
y' = = f ( y,t) (1.1)
dt
Równanie tego typu ma, w ogólnym przypadku, rodzinę rozwiązań. Na przykład,
równanie:
y' (t) = - y(t)
ma następujące rozwiązanie ogólne
y(t) = Ce-t ,
gdzie C jest dowolną stałą.
Konkretne rozwiązanie jest związane z tą trajektorią, na której znajduje się jakieś
rozwiązanie, spełniające określone wymagania, na przykład, warunki początkowe:
y(0) = y0 (rys. 1.1). Można wówczas wyznaczyć stałą C :
y0 = Ce-t = C , zatem: C = y0 .
t=0
W wielu przypadkach obiekt (zjawisko) jest opisywane za pomocą większej liczby
równań różniczkowych. Wówczas do jednoznacznego określenia trajektorii
rozwiązania potrzebne są odpowiednie warunki dla każdego z równań.
Jednoznaczne rozwiązania można uzyskać, jeśli warunki te są dane dla tej samej
wartości zmiennej niezależnej t. Na przykład, dla układu dwóch równań
różniczkowych:
x' = f (x, z,t)
z' = g(x, z,t)
50
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych
potrzebne są dwa warunki początkowe: x(t0 ) , z(t0 ) .
y(t)
t
Rys. 1.1.
Układ równań różniczkowych można zapisać w postaci macierzowej:
dy
= F(y,t) (1.2)
dt
gdzie dla przypadku dwóch równań:
x(t) f (x, z,t)
Ą# ń# Ą# ń#
y = , F(y,t) =
ó#z(t)Ą# ó#g(x, z,t)Ą#
Ł# Ś# Ł# Ś#
A
atwo można pokazać, że równanie n-tego rzędu można sprowadzić do n równań
pierwszego rzędu. Na przykład
2
d u
równanie: = g(u,u',t) można zapisać w postaci następujących równań:
dt2
z' = g(u, z,t),
u' = z.
Podstawowy sposób numerycznego rozwiązywania równań różniczkowych polega
na zastosowanie jakiejś numerycznej procedury całkowania obu stron tego równania.
W tym celu, dla równania o postaci (1.1), funkcja wymuszająca f ( y,t) (prawa strona
równania) powinna zostać poddana dyskretyzacji, co oznacza, że zmienna
niezależna t przyjmuje wartości dyskretne t0 , t1, t2 , .... W zależności od sposobu
dyskretnej reprezentacji funkcji wymuszającej (może tu być zastosowana interpolacja
lub aproksymacja), powstają różne metody rozwiązywania równań różniczkowych.
Ogólnie, metody te dzielą się na jednokrokowe i wielokrokowe w zależności od tego, czy
do obliczania bieżącej próbki rozwiązania korzysta się z wartości funkcji i
rozwiązania odległych o jeden krok do tyłu (metody jednokrokowe), czy też z
historii odległej o większą liczbę kroków (metody wielokrokowe).
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych 51
f(y,t)
tn
t0 t1 t2 tn-2 tn-1 tn tn+1
Rys. 1.2.
7.2. Metody jednokrokowe
Metoda Eulera
Rozpatrzmy równanie (1.1) dla warunków początkowych y0 = y(t0 ) . Załóżmy, że
poszukiwane jest rozwiązanie tego równania dla t = t1 , przy czym: T = t1 - t0 .
Poszukiwane rozwiązanie można przedstawić, korzystając z jej rozwinięcia w szereg
Taylora w pobliżu punktu początkowego t0 :
2
T
y1 = y(t0 + T ) = y(t0 ) + Ty' (t0 ) + y'' (t0 ) + ... (1.3)
2!
Ponieważ, w danym przypadku, dostępna jest tylko pierwsza pochodna
poszukiwanej funkcji, więc do przybliżenia wartości y(t0 + T ) można wykorzystać
pierwsze dwa składniki rozwinięcia (1.3):
y1 H" y(t0 ) + Ty' (t0 ) (1.4)
Zgodnie z (1.1), w miejsce pochodnej można podstawić prawą stronę równania
różniczkowego (funkcję wymuszającą), co prowadzi do następującej zależności:
y1 H" y(t0 ) + Tf (y(t0 ),t0) (1.5)
Powtarzając ten wywód dla kolejnych kroków otrzymamy ogólną zależność:
yn+1 = yn + Tf ( yn ) (1.6)
która jest znana jako jawna metoda Eulera.
52
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych
Określenie  jawna bierze się stąd, że do określenia kolejnej wartości rozwiązania
wykorzystuje się znane wartości z poprzedniego okresu próbkowania
(dyskretyzacji). Ilustracja tej metody jest pokazana na rys. 1.3. Widać tam błąd
wynikający z obcięcia szeregu Taylora: rzeczywista wartość funkcji dla t = tn+1 jest
równa f ( yn+1) , podczas gdy jej aproksymacja według pierwszej pochodnej wynosi
fa ( yn+1) .
f(y,t)
fa(yn+1)
f(yn+1)
f(yn)
T tn
tn tn+1
Rys. 1.3.
Zależność (1.6) można także uzyskać korzystając ze wspomnianej metody
obustronnego całkowania wyjściowego równania (1.1). Wartość rozwiązania dla
t = tn+1 można wówczas określić następująco:
tn+1 tn tn+1
yn+1 = f (y(t))dt = f (y(t))dt + f (y(t))dt (1.7)
+" +" +"
t0 t0 tn
Ponieważ rozwiązanie yn jest znane (gdyż wcześniejsze kroki: t0 , t1,.., tn , zostały już
wykonane, a wartość początkowa f (y0) także jest znana), to:
tn
f (y(t))dt = yn
+"
t0
Pozostaje więc wyznaczyć drugą całkę w (1.7). Przybliżając ją metodą prostokątów
(w którym jeden bok jest utworzony przez odcinek T, a drugi  przez odcinek
f ( yn ) ), otrzymamy zależność (1.6). Ta interpretacja wyjaśnia drugą nazwę
rozpatrywanego algorytmu: metoda prostokątów (w tym przypadku  jawna,
ekstrapolacyjna).
Wartość całki na odcinku T = tn+1 - tn można także przybliżyć prostokątem o boku
określonym przez wartość funkcji w bieżącym punkcie tn+1: f (yn+1). Wówczas,
zależność (1.6) przyjmie zastępującą postać:
yn+1 = yn + Tf ( yn+1) = yn + Ty' (1.8)
n+1
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych 53
która jest znana jako niejawna metoda Eulera (prostokątów). Jak widać, nazwa jest
uzasadniona tym, że postać funkcji (1.7) jest uwikłana, gdyż z obu stron znaku
równości występuje odwołanie do wartości funkcji (lub jej pochodnej) dla tej samej
wartości zmiennej niezależnej (czasu). Spotykana jest także inna nazwa: formuła
interpolacyjna Eulera (prostokątów).
Metoda trapezów
Można zauważyć, że znaki błędów w obu powyższych metodach (jawnej i niejawnej)
są przeciwne (niezależnie od przebiegu funkcji). Można to wykorzystać w celu
zwiększenia dokładności przybliżenia rozwiązania: uzyskuje się to przez obliczenie
średniej z wyników uzyskanych obiema metodami:
f (yn) + f (yn+1) y' + y'
n n+1
yn+1 = yn + T = yn + T (1.9)
2 2
W ten sposób, całka na odcinku T = tn+1 - tn jest określona przez pole trapezu
wyznaczonego przez boki f ( yn ) oraz f ( yn+1) . Można zauważyć, że metoda
trapezów jest także niejawna.
Metody Rungego-Kutty
W odróżnieniu od przedstawionych powyżej metod jednokrokowych, w metodach
Rungego  Kutty nie ma odwołania do pochodnej funkcji, występującej w równaniu
różniczkowym. W jej miejsce występuje odpowiednia kombinacja wartości samej
funkcji, obliczanej w stosownych miejscach. W przypadku najbardziej znanej metody
Rungego-Kutty czwartego rzędu, kolejne przybliżenie rozwiązania jest określane
według następującego wzoru:
1
yn+1 = yn + (F1 + 2F2 + 2F3 + F4) (1.10)
6
gdzie:
F1 = Tf ( yn,t) , F2 = Tf ( yn + F1 / 2,t +T / 2) ,
F3 = Tf (yn + F2 / 2,t + T / 2) , F4 = Tf (yn + F3,t + T ) .
Jest to metoda 4 rzędu, gdyż błąd przybliżonego wzoru (1.10) wynosi O(T5). Metoda
ma charakter jawny, gdyż obliczana wartość yn+1 nie występuje po prawej stronie
formuły (1.10).
Realizacja tej metody wymaga wykonania w każdym kroku czasowym obliczeń
czterech wartości funkcji prawej strony równania różniczkowego f (y,t) , a więc,
funkcja ta musi być dostępna w jawnej formie. Przy spełnieniu tych wymagań,
omawiana metoda jest prosta w realizacji i zapewnia dużą dokładność rozwiązania.
Jest najczęściej stosowną metodą w zastosowaniach inżynierskich i naukowych. Bez
żadnych dodatkowych ograniczeń może być stosowana do rozwiązywania układów
równań różniczkowych, również nieliniowych. Profesjonalne programy z tą metodą
rozwiązywania równań różniczkowych są najczęściej wyposażone w mechanizm
automatycznego doboru długości kroku całkowania T, co może przyśpieszyć czas
54
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych
obliczeń przy zachowaniu założonej dokładności. Niestety, przy stosowaniu tej
metody mogą wystąpić problemy ze stabilnością w przypadku układów sztywnych
(gdy w modelowanym systemie występują znacznie różniące się stałe czasowe.
Dokładność metody
Omówione powyżej metody można zapisać następująco:
1
y' = (yn - yn-1) jawna metoda prostokątów,
n
T
1
y' = (yn+1 - yn) niejawna metoda prostokątów,
n+1
T
1 1
(y' + y' )= (yn+1 - yn) metoda trapezów (jest to także metoda niejawna).
n+1 n
2 T
Ogólna postać tych związków przybiera następującą formę:
a1yn+1 + a0 yn = T(b1y' + b0 y' ) (1.11)
n+1 n
Współczynniki tego równania określają odpowiedni algorytm:
a1 = 1, a0 = -1 - dla wszystkich metod,
b1 = 0 , b0 = 1 - jawna metoda prostokątów,
b1 = 1, b0 = 0 - niejawna metoda prostokątów,
1
b1 = b0 = - metoda trapezów.
2
Na podstawie równania (1.11) można przeprowadzić dyskusję dokładności
rozpatrywanych metod. Analogicznie do (1.3) napiszemy:
2
T T3
yn+1 = y(tn + T ) = y(tn) +Ty'(tn) + y''(tn) + y'''(tn) + ... (1.12)
2! 3!
Po zróżniczkowaniu:
2
T T3
y' = y'(tn +T ) = y'(tn) + Ty''(tn) + y'''(tn) + y''''(tn) + ... (1.13)
n+1
2! 3!
Podstawiając powyższe zależności do formuły (1.11) otrzymamy:
2 2
�# ś# �# ś#
T T3 T
ś# ź# ś# ź#
a1ś# yn + Ty' + y'' + y''' +...ź# + a0 yn = Tb1ś# y' + Ty'' + y''' +...ź# + Tb0 y' (1.14)
n n n n n n n
2 6 2
�# # �# #
skąd:
1 1 1
2 3
yn(a1 + a0)+ Ty' (a1 - b1 - b0) = -T y'' �# a1 - b1 ś# -T y''' �# a1 - b1 ś# -... (1.15)
n n n
ś# ź# ś# ź#
2 6 2
�# # �# #
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych 55
Równanie (1.15) jest spełnione dla dowolnych wartości funkcji i jej pochodnych
wtedy, gdy współczynniki określone przez a0 , a1, b0 , b1 będą równe zero, a więc:
a1 + a0 = 0
a1 - b1 - b0 = 0
1
a1 - b1 = 0
2
1 1
a1 - b1 = 0
6 2
Dla metody prostokątów ( a1 = 1, a0 = -1 oraz b1 = 0 , b0 = 1 lub b1 = 1, b0 = 0 ),
pierwszy niezerowy współczynnik stoi przy drugiej pochodnej funkcji:
1
2 2
T y'' �# �"1-1ś# = C2T y'' ,
n n
ś# ź#
2
�# #
1 1
skąd: C2 = lub C2 = - .
2 2
1
Podobnie, dla metody trapezów ( a1 = 1, a0 = -1 oraz b1 = b0 = ) pierwszy niezerowy
2
współczynnik stoi przy trzeciej pochodnej funkcji:
1 1 1 1
ś#
3 3 3
-T y''' �# �"1- �" = T y''' + ... = C3T y''' +...
n n n
ś# ź#
6 2 2 12
�# #
1
skąd: C3 = .
12
Powyższe związki charakteryzują dokładność metody zgodnie z następującymi
regułami.
- Rząd metody p określa się rzędem pochodnej (p+1), dla której pierwszy
współczynnik jest różny od zera. Zatem: dla metody prostokątów p=1; dla
metody trapezów p=2;
- Błąd odcięcia jest związany z wartościami wyrazów, które nie spełniają
równości (1.15). Współczynnik stojący przy najniższej pochodnej, który nie
spełnia tego równania jest właśnie błędem odcięcia. Dla metody trapezów jest
on równy: C3 = Cp+1 = 1/12 .
Stabilność metody
Do badania stabilności określonego algorytmu numerycznego rozwiązywania
równań różniczkowych wygodnie jest posługiwać się wybranym równaniem
modelowym. Rozpatrzmy równanie o postaci:
y'(t) = y(t) (1.16)
z warunkiem początkowym y(0) = y0 , przy czym, współczynnik  jest w ogólnym
przypadku zespolony:  = u + jw ,
jwt
Rozwiązaniem równania (1.16) jest funkcja wykładnicza: y(t) = et = eute .
56
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych
1. Zastosujmy do rozwiązania równania (1.16) jawną metodę prostokątów. Załóżmy,
że dyskretna postać zmiennej niezależnej t dostępna jest ze stałym przedziałem T.
Rozwiązanie w kolejnych krokach przybiera następujące wartości:
y1 = y0 + Ty' = y0 + Ty0 = (1+ T)y0 ,
0
2
y2 = y1 + Ty' = (1+ T)y0 + Ty1 = (1+ T)y0 + T(1+ T)y0 = (1+ T) y0 ,...
1
i ogólnie:
n
yn = (1+ T) y0 . (1.17)
Jeśli wyjściowe równanie ciągłe jest stabilne (u = Re() < 0 ), to uzyskana
aproksymacja dyskretna rozwiązania także powinna być stabilna. W odniesieniu do
(1.17) jest to równoważne warunkowi:
1+ T d" 1 (1.18)
Zatem:
2 2
1+ Tu + jTw d" 1, czyli: (1+ Tu) + (Tw) d" 1.
Równanie powyższe przedstawia okrąg na płaszczyz
nie zespolonej o promieniu
jednostkowym i o środku leżącym w punkcie (-1, 0) (jeśli współrzędne są
przeskalowane: (Tu, Tw)). Obszar stabilności leży wewnątrz okręgu (rys. 1.4). Można
zauważyć, że przy danej wartości  obszar stabilności zależy od długości przedziału
T: im większa jest wartość T tym mniejszy obszar stabilności (granica stabilności:
Tu=Tw=1).
Re(T)
Tw
Re(T)
-1
Tu
Rys. 1.4.
2. Powtórzmy powyższe rozważania dla niejawnej metody prostokątów. Tym razem,
rozwiązanie w kolejnych krokach jest określone następująco:
1
y1 = y0 +Ty' . Stąd: y1 = y0
1
1-T
2
1 1
�# ś#
y2 = y1 + Ty' = y0 + Ty2 . Stąd: y2 = y0
2
ś# ź#
1-T
�#1-T #
i ogólnie:
n
1
�# ś#
yn = y0 . (1.19)
ś# ź#
�#1- T #
Ciąg ten jest ograniczony dla:
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych 57
1
d" 1, (1.20)
1-T
co jest równoważne :
2 2
(1-Tu) + (Tw) e" 1
Równanie to jest spełnione dla całej płaszczyzny zespolonej za wyjątkiem wnętrza
okręgu (przy przeskalowanych osiach Tu, Tw) o środku w punkcie (1,0)  rys. 1.5. A
zatem, przy danej wartości  , procedura mogłaby być niestabilna dla małych
wartości T. Jednak, w rzeczywistych warunkach, zawsze u d" 0 (warunek stabilności
równania ciągłego), co oznacza, że niejawna procedura Eulera numerycznego
rozwiązywania równania (1.16) jest globalnie stabilna, jeśli tylko wyjściowe
równanie jest stabilne.
Re(T)
Tw
Re(T)
1
Tu
Rys. 1.5.
3. Dla formuły trapezów mamy:
T 2 + T T 2 + T
y1 = y0 + (y0 + y1) = y0 + y1 , skąd: y1 = y0 .
2 2 2 2 -T
Ogólnie:
n
2 +T
�# ś#
yn = y0 (1.21)
ś# ź#
2 -T
�# #
Ciąg ten jest stabilny dla:
2 + T
d" 1, (1.22)
2 -T
co jest równoważne relacji:
2 2
2 + Tu + jTw (2 + Tu) + (Tw)
d" 1 oraz: d" 1
2 2
2 -Tu + jTw
(2 -Tu) + (Tw)
Powyższy związek jest słuszny dla: u d" 0 , co oznacza, że procedura jest stabilna dla
całej ujemnej części płaszczyzny zespolonej współczynnika  (rys. 1.6) - a więc
procedura jest również stabilna globalnie, jeśli tylko wyjściowe równanie ciągłe jest
stabilne.
58
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych
Re(T)
Tw
Re(T)
Tu
Rys. 1.6.
7.3. Metody wielokrokowe
Metody Geara
Spośród metod wielokrokowych, szczególnie ważne są metody całkowania,
projektowane z myślą o zastosowaniu w odniesieniu do tzw. sztywnych systemów.
Sztywność układu równań różniczkowych oznacza szybkie zanikanie (dążenie do
zera) rozwiązania. Stosowane w takim przypadku metody nie mogą być wysokiego
rzędu, które należą do grupy metod niejawnych. Są to przede wszystkim metody
Geara i niejawne metody Rungego Kutty (R K)1, które odznaczają się dużą
stabilnością [5].
Metody Geara wywodzą się bezpośrednio z numerycznej reprezentacji pochodnej
w podstawowym równaniu (1.1), Na przykład, przy najprostszym zapisie
pochodnej:
dy(t) yn+1 - yn
y' = H"
dt T
oraz przy założeniu, że funkcja prawej strony równania będzie liczona według
zasady ekstrapolacji: f (yn+1,tn+1) , to otrzymamy:
yn+1 = yn + Tf (yn+1,tn+1) , (1.23)
co pokrywa się z niejawną metodą prostokątów.
W celu podwyższenia rzędu metody, można zastosować dokładniejszą zależność
na obliczanie pochodnej. W tym celu można zastosować wielomian interpolacyjny
Newtona, w którym pochodna na końcu ostatniego przedziału jest obliczana na
podstawie trzech punktów, w których znana jest różniczkowana funkcja. Prowadzi
to do następującej zależności (p. 4.3):
dy(t) 3yn+1 - 4yn + yn-1
y' = H"
dt 2T
1
Należy zauważyć, że powszechnie są stosowane jawne metody R K (p. 7.2), które nie mają
wymaganych tu właściwości w odniesieniu do układów sztywnych, np. metoda R K IV rzędu
[5], [14].
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych 59
Podstawienie tej zależności do (1.1) prowadzi do następującej formuły Geara 2-go
rzędu:
4 1 2
yn+1 = yn - yn-1 + Tf (yn+1,tn+1) (1.24)
3 3 3
Metody wyższych rzędów nie są stosowane, gdyż występują również problemy z
ich stabilnością w przypadku sztywnych systemów równań.
Niejawna metoda Rungego-Kutty
Prezentowane w p. 7.2 jednokrokowe algorytmy Rungego-Kutty są metodami
jawnymi i w związku z tym zazwyczaj nie są stosowane do modelowania układów
dynamicznych. Odrębną grupę stanowią wielokrokowe niejawne metody Rungego-
Kutty [2]. Spośród tych ostatnich zachęcający jest, zwłaszcza w algorytmach
modelowania numerycznego, 2-stopniowy algorytm R K II rzędu. W literaturze
anglojęzycznej jest on oznaczany akronimem: 2S-DIRK (ang. 2-stage diagonally
implicite Runge Kutta). Numeryczne przybliżenie rozwiązania równania (1.1) w k-
tym kroku jest określane za pomocą następującego 2-stopniowego algorytmu [5]:
~
~(k)
tk ~(k)
1. y = y(k -1) + Tf (~ , y ) (1.25)
~(k ~(k)
aproksymacja: y -1) = ą y(k -1) + � y (1.26)
~(k ~
2. y(k) = y -1) + Tf (tk , y(k)) (1.27)
gdzie zmienne oznaczone tyldą są wielkościami pomocniczymi, odnoszącymi się do
punktu pośredniego w czasie, pomiędzy tk 1 i tk:
~ ~ 1
�#1+ ś#T
~
tk = tk-1 +T , T = (1.28)
ś# ź#
ą
�# #
ą = - 2 , ą + � =1.
Można zauważyć, że w obu stopniach algorytmu jest realizowana niejawna
~(k
metoda Eulera z początkowymi punktami w y(k -1) i y -1) , odpowiednio, z
~
krokiem T = (1 2 / 2 )T. A
atwo więc wyznaczyć odpowiednie zależności odnoszące
się do konkretnego zastosowania metody.
7.4. Metody ekstrapolacyjno-interpolacyjne
Powyższe rozwiązania pokazują, że metody niejawne rozwiązywania równań
różniczkowych mają szersze zastosowanie, gdyż wykazują większą stabilność.
Jednak formuła niejawna nie zawsze może być stosowana, gdyż w ogólnym
przypadku, wymaga rozwiązania równania uwikłanego (niewiadoma występuje po
obu stronach równania), co niekiedy może być kłopotliwe. W takim przypadku
można stosować dwie różne formuły w każdym kroku rozwiązania równania,
według następującego schematu:
~n+1
1. Obliczyć przybliżone rozwiązanie y według procedury jawnej.
60
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych
2. Poprawić rozwiązanie yn+1 według procedury niejawnej, z wykorzystaniem
rozwiązania przybliżonego.
Pierwszy krok jest prognozą (ekstrapolacją) rozwiązania, a drugi: korekcją
(interpolacją), skąd schemat ten jest nazywany metodą prognozy i korekcji.
W pierwszym kroku stosuje się metodę niejawną o jeden rząd mniej dokładną, niż
w drugim kroku z metodą niejawną. Na przykład, niejawna metoda prostokątów (I
rząd) może być połączona z metodą trapezów (niejawna metoda II rzędu):
~n+1 = yn + Tf
y (yn,tn)
(1.29)
T
yn+1 = yn + f (yn,tn)+ f (~
yn+1,tn+1)
2
Zazwyczaj są w takim przypadku stosowane metody wyższych rzędów.
8. Literatura
[1] AL-KHAFAJI A.W., TOOLEY J.R., Numerical methods in engineering practice,
Holt, Rinehart and Winston, Inc., New York,1986.
[2]  #(/ .., 
" !.$., '8A;5==K5 <5B>4K @5H5=8O
>1K:=>25==KE 48DD5@5=F80;L=KE C@02=5=89 (7040G0 >H8). Dostępne w:
http://www.srcc.msu.su/num_anal/list_wrk/sb3_doc/part6.htm
[3] DRYJA M., JANKOWSCY J. i M., Przegląd metod i algorytmów numerycznych.
Część 2, WNT, Warszawa, 1982.
[4] FORSYTHE G.E., MALCOLM M.A., MOLER C.B., Computer methods for
mathematical computations, Englewood Cliffs, N.J., Prentice-Hall, Inc., 1977.
[5] FORTUNA Z., MACUKOW B., W
SOWSKI J., Metody numeryczne, WNT,
Warszawa, 2003.
[6] JANKOWSCY J. i M., Przegląd metod i algorytmów numerycznych. Część 1,
WNT, Warszawa, 1981.
[7] KACZOREK T., Wektory i macierze w automatyce i elektrotechnice, WNT,
Warszawa, 1998.
[8] KIEA
BASIC
SKI A., SCHWETLICK H., Numeryczna algebra liniowa, WNT,
Warszawa, 1992.
[9] KINCAID D., CHENEY W., Analiza numeryczna. WNT, Warszawa, 2005.
[10] KRUPKA J., MORAWSKI R.Z., OPALSKI L.J., Wstęp do metod numerycznych
dla studentów elektroniki i technik informacyjnych. Oficyna Wydawnicza
Politechniki Warszawskiej, Warszawa 1999.
[11] PRESS W.H., TEUKOLSKY S.A., VETTERLING W.T., FLANNERY B.P.,
Numerical recipes in C. The art of scientific computing. Cambridge University
Press, Cambridge 1992 (oddzielne fragmenty książki dostępne są na stronie
internetowej: http://apps.nrbook.com/empanel/index.html# ).
[12] ROSOA
OWSKI E., Cyfrowe przetwarzanie sygnałów w automatyce
elektroenergetycznej. Akademicka Oficyna Wydawnicza EXIT, Warszawa 2002.
[13] ROSOA
OWSKI E., Komputerowe metody analizy elektromagnetycznych stanów
przejściowych. Oficyna Wydawnicza Politechniki Wrocławskiej, Wrocław
2009.
[14] STOER J., BULRISCH R., Wstęp do analizy numerycznej, PWN, Warszawa,
1987.
[15] WANAT K., Wybór metod numerycznych, Wydawnictwo DIR, Gliwice, 1993.
[16] WEISSTEIN E.W., Singular Value Decomposition. From MathWorld--A Wolfram
Web Resource:
http://mathworld.wolfram.com/SingularValueDecomposition.html
Podstawowe idee związane z aproksymacją średniokwadratową zostały
sformułowane niezależnie przez dwóch matematyków:
Johann Carl Friedrich Gauss (1777-1855)  Niemiec,
62
Literatura
Adrien-Marie Legendre (1752-1833)  Francuz.
Inne ciekawe strony:
http://www.csit.fsu.edu/~burkardt/
http://public.lanl.gov/mewall/kluwer2002.html
http://www.american.edu/cas/mathstat/People/kalman/pdffiles/svd.pdf
http://www.cs.toronto.edu/NA/index.html - Department of Computer Science,
University of Toronto.
Skorowidz
Metoda eliminacji Gaussa ................... 10 Uwarunkowanie zadania .................... 7
Stabilność numeryczna algorytmu .... 7 Złożoność obliczeniowa ...................... 7