Table Of ContentTHE COLLEGE OF COMPUTER SCIENCE
E dw ard K ącki, A ndrzej M ałolepszy
A licja R om anow icz
M etody num eryczne dla
inżynierów
W yższa Szkoła In fo rm atyki
i v Łodzi
Spis treści
1. Informacje wstępne 7
1.1. Reprezentacje liczb 7
1.2. Rachunek zmiennopozycyjny 9
1.3. Algorytmy numeryczne 10
1.4. Numeryczna stabilność algorytmów 14
2. Obliczanie wartości funkcji 16
2.1. Obliczanie wartości wielomianu 16
2.2. Obliczanie wartości wielomianu trygonometrycznego 19
3. Aproksymacja 23
3.1. Sformułowanie zagadnienia 23
3.2. Wielomiany ortogonalne 26
3.3. Wielomiany ortogonalne na zbiorze dyskretnym 29
3.4. Aproksymacja w przestrzeniach Hilberta 30
3.5. Aproksymacja na zbiorze dyskretnym 32
3.6. Przykłady 34
4. Interpolacja 37
4.1. Interpolacja wielomianowa 37
4.2. Wielomian interpolacyjny Lagrange’a 38
4.3. Wzór interpolacyjny Lagrange’a dla węzłów
równoodległych 40
4.4. Metoda Aitkena 41
4.5. Wzór interpolacyjny Newtona 44
4.6. Wzór interpolacyjny Hermite’a 46
4.7. Przykłady 48
5. Rozwiązywanie układów algebraicznych równań liniowych 49
5.1. Metody dokładne 49
5.1.1. Metoda eliminacji Gaussa 49
5.1.2. Oszacowanie błędu 56
5.1.3. Przykłady 59
5.2. Metody iteracyjne 63
5.2.1. Metoda iteracji prostej 64
3
5.2.2. Metoda Seidla 69
6 . Rozwiązywanie równań nieliniowych oraz ich układów 70
6.1. Metody wyznaczania pierwiastków w zadanym
przedziale 70
6.1.1. Metoda połowienia przedziału 70
6.1.2. Metoda cięciw zwana również regułą falsi 73
6.1.3. Metoda stycznych zwana również metodą
Newtona 76
6.1.4. Metoda mieszana 77
6.2. Rozwiązywanie układów równań nieliniowych 80
6.2.1. Metoda Newtona 80
6.2.2. Metoda iteracji 84
7. Całkowanie numeryczne 89
7.1. Kwadratury liniowe 89
7.2. Kwadratury Newtona-Cotesa 91
7.3. Złożone kwadratury Newtona-Cotesa 94
7.4. Kwadratury Gaussa 98
7.5. Obliczanie całek niewłaściwych i całek z
osobliwościami 104
7.6. Obliczanie numeryczne całek wielokrotnych 110
8 . Rozwiązywanie równań różniczkowych zwyczajnych 112
8.1. Metody różnicowe 112
8.1.1. Metody jednokrokowe 112
8.1.2. Liniowe metody wielokrokowe 114
8.1.3. Metoda współczynników nieoznaczonych 117
8.1.4. Metody iteracyjne 119
8.1.5. Obliczanie wartości wstępnych 122
8.2. Metody Rungego-Kutty 123
8.3. Równania n-tego rzędu i układy równań
różniczkowych zwyczajnych 128
8.3.1. Metoda Rungego-Kutty 128
8.3.2. Metody ekstrapolacyjne 133
8.3.3. Metody interpolacyjne 135
8.3.4. Metoda różnicowa zwykła i ulepszona 138
8.4. Metody wariacyjne 142
4
9. Rozwiązywanie równań różniczkowych cząstkowych 151
9.1. Metody różnicowe 151
9.1.1. Ogólne uwagi o dyskretyzacji i funkcjach
siatkowych 152
9.1.2. Równanie Laplace’a 157
9.1.3. Równanie Poissona 163
9.1.4. Równanie Helmholtza 164
9.1.5. Równanie przewodnictwa 165
9.1.6. Równanie fali płaskiej 171
9.1.7. Równanie telegrafistów 173
9.1.8. Ulepszona metoda różnicowa 174
9.2. Stabilność schematów różnicowych 175
9.2.1. Uwagi wstępne 175
9.2.2. Stabilność schematów różnicowych dla
równania typu parabolicznego 178
9.2.3. Stabilność schematu różnicowego dla
równania typu hiperbolicznego 186
9.3. Metody Galerkina 190
9.4. Przykłady 197
9.5. Metoda elementu skończonego 205
10. Rozwiązywanie równań całkowych 212
10.1. Metoda kolejnych przybliżeń 212
10.2. Metoda sum skończonych 215
11. Metody Monte Carlo 217
11.1. Rozwiązywanie równań z jedną niewiadomą 217
11.2. Obliczanie całek oznaczonych 220
5
1. Informacje wstępne
1.1 Reprezentacja liczb
W komputerach liczby są reprezentowane przez skończoną liczbę
cyfr ich rozwinięć pozycyjnych o pewnej podstawie rozwinięcia.
Najczęściej stosowanymi podstawami rozwinięć liczb są podstawy 2, 8
lub 16. Arytmetyka tych rozwinięć jest podobna, dlatego też możemy
ograniczyć się do arytmetyki dwójkowej inaczej zwanej binarną. W
systemie dwójkowym mamy dwie cyfry 0 lub 1. W komputerach na
reprezentację liczby przeznacza się słowo o skończonej długości np. d +
1 bitów (cyfr dwójkowych). Sposób rozwinięcia liczby zależy od jej typu.
Liczby całkowite przedstawiane są w sposób stałopozycyjny.
Dowolną liczbę całkowitą I w systemie dwójkowym możemy przedstawić
w postaci
1 = e' 2‘ ’ (1-1)
i=0
gdzie zn - jest znakiem liczby /, e, - są jej cyframi dwójkowymi, przy
czym jeśli I * 0, to en^Q. W komputerze na reprezentację liczby
stało pozycyjnej jest przeznaczone słowo o długości d + 1 bitów i jeden
bit jest przeznaczony na znak liczby, zatem liczba całkowita I może być
przedstawiona w postaci (1.1) jeśli / e [-2d + 1, 2d - 1]. Oznacza to, że
n <d i mówimy wtedy, że liczba jest reprezentowana dokładnie.
Przykład 1.1. Liczba I = 18 w systemie dwójkowym może być
przedstawiona w postaci / = 1 ■ 24 + 0 • 23 + 0- 22 + 1 • 21 + 0 • 2°. Może
więc być przechowywana w słowie o długości d + 1 > 5 bitów jako
000...010010
--------V---------'
2« ci cyfr rozwinięcia
Liczby rzeczywiste są przedstawione w sposób zmiennopozycyjny.
Dowolną liczbę rzeczywistą x ^ 0 (zero jest na ogół zapamiętane jako
słowo o wszystkich bitach zerowych) możemy przedstawić w postaci
x = zn ■ 2C m, (1.2)
gdzie zn - jest znakiem liczby, c - jest liczbą całkowitą zwaną cechą, m
- jest liczbą rzeczywistą z przedziału [y,l) zwaną mantysą. Liczba
7
rzeczywista jest więc przedstawiona za pomocą dwu grup bitów. Jeśli
liczba x jest przechowywana w jednym słowie o długości d + 1 bitów,
wtedy jeden bit przeznaczony jest na znak liczby, t bitów na
reprezentację mantysy m, a pozostałe d - t bitów na reprezentację
cechy c, przy czym jeden z bitów tej reprezentacji musi być
przeznaczony na jej znak.
Mantysa jako liczba rzeczywista z przedziału [y,l), może być
przedstawiona w systemie dwójkowym w postaci nieskończonego
rozwinięcia
oo
m = Yue-i2 ' > (1.3)
gdzie e_, = 1 oraz e_; równe 0 lub 1 dla /> 1.
Jednak na jej zapamiętanie jest przeznaczonych t bitów, może być
zatem zapamiętane tylko jej t początkowych cyfr dwójkowych, tzn.
mantysa zostaje prawidłowo zaokrąglona do t cyfr w następujący
sposób:
(1.4)
7 = 1
W wyniku zaokrąglenia, a więc zastąpienia m (1.3) przez mt (1.4)
mamy
\m - < -i- • 2 1. (1.5)
Natomiast reprezentacja zmiennopozycyjna liczby rzeczywistej x danej
wzorem (1.2 ) oznaczana jest symbolem rd{x) i jest równa
rd(x) = zn • 2C mt. (1.6)
Ze wzorów (1.2), (1.4) i (1.6) otrzymujemy, że dla x #0
rd(x) - x <
x
co oznacza, że
rd(x) = x(l + £■), gdzie | e | < 2”'. (1.7)
Liczbę 2“' nazywamy dokładnością maszynową.
Przykład 1.2. Liczba x = 18,5 daje się zgodnie ze wzorem (1.2)
przedstawić w postaci
x = 2(1'22+a2l+,-20)(l • T x + 0 • T 1 + 0 ■ 2“3 + 1 ■ 2“4 + 0 ■ 2"5 + 1 • 2“6) .
Może być więc przechowywana w słowie o długości d + 1 >11 bitów,
o t> 6 bitach mantysy jako
0 0 0 0 ...0 1 0 1 1 0 0 1 0 1 0 ...0 .
zn znc /bitów mantysy
d-t bitów cechy
W przedstawieniu zmiennopozycyjnym liczby rzeczywistej liczba cyfr
mantysy decyduje o dokładności przedstawienia liczby, natomiast liczba
cyfr cechy decyduje o zakresie liczb, które przy pomocy takiego
przedstawienia mogą być reprezentowane. Przypadek kiedy cecha jest
za duża nazywany jest nadmiarem cechy i w większości komputerów
taki przypadek będzie sygnalizowany jako błąd, natomiast przypadek
kiedy cecha jest za mała nazywamy niedomiarem cechy i zwykle nie jest
on sygnalizowany. Na ogół w przypadku niedomiaru cechy liczba
reprezentowana jest zerem i nie jest wtedy spełniony warunek (1.7),
ponieważ błąd względny w tym przypadku jest równy 100%.
We współczesnych komputerach nadmiary i niedomiary cechy są
stosunkowo rzadkie, ponieważ na przedstawienie cechy rezerwowana
jest wystarczająco duża ilość cyfr, a także można ich unikać poprzez
odpowiednie przeskalowanie danych lub też modyfikacje algorytmów
obliczeniowych.
1.2 Rachunek zmiennopozycyjny
Działania arytmetyczne w komputerze wykonywane są na ogół na ich
reprezentacjach zmiennopozycyjnych. Realizację działań w arytmetyce
zmiennopozycyjnej oznacza się skrótem fi Wynik operacji
arytmetycznych nie musi być reprezentowalny w słowie o długości d + 1
cyfr pomimo, że argumenty tej operacji są liczbami reprezentowalnymi.
Można oczekiwać, że zamiast operacji arytmetycznych
realizowane są operacje zastępcze tzw. operacje zmiennopozycyjne,
które możemy oznaczyć następująco:
Jl(x + y) = rd(x + y),
fl(x-y) = rd{x-y), (1.8)
fl(x *y) = rcl(x * y),
fl(x t y) = rd(x / y),
gdzie fl(W) oznacza wynik zmiennopozycyjnego działania na
argumentach x i y równych swoim reprezentacjom x = rd(x), y = rd(y).
Nadmienić należy, że dla arytmetyki zmiennopozycyjnej fałszywe są
prawa łączności i rozdzielczości.
Przykład 1.3. Rozpatrzmy dwie sumy
jl(a + (b + c)) = (a + (b + c)(l + £■, ))(1 + s2) =
= (a + b + c)f 1 + - —-](1 + s2) = (a + b + ć){\ + 5,),
V a + b + cJ
fl((a + b) + c) = ((a + 6)(1 + e3) + c)( l + s4) =
b)s3
= (a + b + c)I 1 + - 1(1 + s4) = (a + b + c)( 1 + ó2). ’
a + b + c
gdzie
«. , (b + c)ex(\ + e2) e t (a + 6)Ą (l + f4)
o, — 1 + f, + , a, — 1 + + ,
a+6+c " a+ó+c
a zatem przy różnej kolejności sumowania wytworzone błędy c), i S2 są
wyrażone różnymi wzorami.
1.3 Algorytmy numeryczne
Załóżmy, że mamy skończoną liczbę danych rzeczywistych
x = (x,,x2,... ,xn) i na ich podstawie chcemy obliczyć skończenie wiele
wyników rzeczywistych y = (y1,y2,... ,ym) . Będziemy więc chcieli
określić wartość y według odwzorowania
y = <p(x), (1.9)
gdzie (p\D -» R'" jest ciągłym odwzorowaniem, i D ci?".
Algorytmem będziemy nazywać jednoznaczny przepis obliczania
wartości odwzorowania cp(x), przy czym przepis ten składa się ze
skończonej sekwencji odwzorowań elementarnych.
Jeśli takie zadanie rozwiązujemy numerycznie, to zamiast z
dokładnymi danymi x,,x2,... ,xn, (x,. <= R), mamy do czynienia z
reprezentacjami rd(x]), rd(x2) , ... , rd(xn) określonymi wzorem
(1.7). Operacje elementarne nie mogą więc być realizowane dokładnie
tylko poprzez odwzorowania zastępcze fl. Ta niewielka zmiana danych
10
powoduje, że różne algorytmy rozwiązania tego samego zadania dają
na ogół niejednakowe wyniki. Dużą rolę odgrywa tutaj przenoszenie się
błędów zaokrągleń, co widać w przykładzie 1.3. W zależności od tego,
czy mniejsza jest wartość \a + Ą czy też \b + cj, korzystniej jest
utworzyć sumę a + b + c według wzoru (a + b) + c lub też wzoru
a + (b + ć).
Często też niewielkie zmiany danych, powodują duże względne
zmiany rozwiązania zadania. Takie zadanie nazywamy źle
uwarunkowanym. Wielkości charakteryzujące wpływ zaburzeń danych
na zaburzenia rozwiązania nazywamy wskaźnikami uwarunkowania
zadania.
Przykład 1.4. (J.H. Wilkinson) Rozpatrzmy wielomian
20
w(x) = a20x20 + a]9x'9 H------b a{x + aQ = ]^[ (x - k).
k=\
Zerami tego wielomianu są liczby naturalne 1,2, ... ,20. Współczynnik
np. a]g jest równy a19 = -210. Zaburzmy ten współczynnik
wprowadzając nowy współczynnik al9 = -(2 1 0 + 2~23) , a więc
al9(e) = a19(l - £•), gdzie s - (210- 223)”1 < 2~30. Otrzymamy wtedy
wielomian w£(x) = w(x) - 2~2V 9, który posiada już pierwiastki
zespolone i np. najbliższe pierwiastkowi 15 wielomianu w(x), są
pierwiastki 13.992358137 ±2.518830070/. A więc niewielka zmiana tylko
jednego współczynnika wielomianu, spowodowała duże zmiany jego zer.
Z przykładu 1.4. wynika, że ważnym zagadnieniem jest badanie
uwarunkowania zadania numerycznego, a więc badanie jaki wpływ na
otrzymane wyniki y = (yx,y2, ,ym) ma zaburzenie danych
x = (x15x2, ... ,x j , szacując względną zmianę wyniku.
Metodę badania przenoszenia się błędów , można rozbudować do
różniczkowej analizy błędów algorytmu (1.9). Niech 5xi = rd(xi) - x j,
dla i = 1, 2 n, wtedy
" drp. (x)
fy = (p(rd(x)) - ę(x) = ^ 5xt , (1.10)
/=i
dla j = 1, 2, ... , m. Wzór ten otrzymujemy ze wzoru Taylora po
odrzuceniu wielkości wyższego rzędu. We wzorze (1.10)
11
współczynnikiem określającym wrażliwość y, na bezwzględną zmianę
d(pix)
8xi wartości x,. jest — -—
ćki
Analogiczny wzór można wyprowadzić dla przenoszenia się błędów
względnych. Jeśli y . * 0 dlay = 1, 2 m i xt, ^ 0, dla /' = 1, 2, ... , n,
to
x. dę ix)
We wzorze (1.11) współczynnik — j——i— nazywamy wskaźnikiem
(Pj{x) dxt
uwarunkowania. Wskaźnik uwarunkowania określa w jakim stopniu błąd
względny wartości x,. wpływa na błąd względny wartości y .. Jeśli
wskaźniki uwarunkowania są co do wartości bezwzględnej duże, to
zadanie jest źle uwarunkowane.
Przykład 1.5. Zbadajmy uwarunkowanie zadania obliczania sumy
y - a + b + c, gdzie a,b,c eR. Oszacujmy względną zmianę wyniku
dla zaburzonych danych
a( 1 + £■,) + b{ 1 + £,) + c(l + £j) - a - b - c a ■ £j + b ■ £, + c ■ £3
<
a + b + c a + b + c
a
< £ + +
a + b + c a + b + c a + b + c
Wskaźnikiem uwarunkowania tego zadania jest
a b c
+ +
a + b + c a + b + c a + b + c
Wnioskujemy skąd, że zadanie jest dobrze uwarunkowane, jeśli każdy
składnik sumy jest mały w stosunku do sumy y = a + b+c. W
przypadku, gdy wszystkie składniki sumy są tego samego znaku,
wskaźnik uwarunkowania jest równy 1. Maksymalne zaburzenie danych
może się przenieść na zaburzenie względne wyniku co najwyżej z takim
mnożnikiem, a więc zadanie jest bardzo dobrze uwarunkowane, tzn.
mało wrażliwe na zaburzenie danych.
12
