Metody numeryczne
Bartłomiej Bułat, Michał Nowak, Konrad Malawski
12 grudnia 2009
Spis treści
1 Reprezentacja stało- i zmiennopozycyjna 2
1.1 Reprezentacja zmiennopozycyjna . . . . . . . . . . . . . . . . . . 2
1.2 BÅ‚Ä…d reprezentacji zmiennoprzecinkowej . . . . . . . . . . . . . . 2
1.3 Reprezentacja stałopozycyjna . . . . . . . . . . . . . . . . . . . . 2
2 Uwarunkowanie zadania 2
3 Stabilność i poprawność numeryczna 3
3.1 Stabliność numeryczna algorytmów . . . . . . . . . . . . . . . . . 3
3.2 Poprawność numeryczna algorytmów . . . . . . . . . . . . . . . . 3
3.3 Złożoność obliczeniowa algorytmów . . . . . . . . . . . . . . . . . 3
4 Miejsca zerowe funkcji nieliniowych 3
4.1 Metoda bisekcji (połowienia przedziałów) . . . . . . . . . . . . . 3
4.2 Metoda Newtona (stycznych) . . . . . . . . . . . . . . . . . . . . 4
4.2.1 Szybkość zbieżności metody . . . . . . . . . . . . . . . . . 4
4.3 Metoda siecznych . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
4.4 Metoda regula falsi (tylko na ćwiczeniach) . . . . . . . . . . . . . 4
5 Interpolacja 5
5.1 Interpolacja Lagrange a . . . . . . . . . . . . . . . . . . . . . . . 5
5.2 Interpolacja Newtona . . . . . . . . . . . . . . . . . . . . . . . . . 6
5.3 Interpolacja Hermite a . . . . . . . . . . . . . . . . . . . . . . . . 6
5.4 Efekt Rungego . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
6 Interpolacja funkcjami sklejanymi 7
6.1 Interpolacja funkcjami sklejanymi 1 stopnia . . . . . . . . . . . . 8
6.2 Interpolacja funkcjami sklejanymi 3 stopnia . . . . . . . . . . . . 8
7 Licencja 9
1
1 Reprezentacja stało- i zmiennopozycyjna
1.1 Reprezentacja zmiennopozycyjna
Znormalizowana reprezentacja dziesiętna (x " R, cecha: C " Z, mantysa: 0,1
|m| < 1, znak: s " { 1, 1}):
C
x = s m 10
1
Znormalizowana reprezentacja dwójkowa ( m < 1):
2
C
x = s m 2
1.2 BÅ‚Ä…d reprezentacji zmiennoprzecinkowej
C
Reprezentacja liczby rzeczywistej: rd(x) = s m 2, gdzie t to ilość bitów
t
C
mantysy, zatem x = s m jest dokładną liczbą x (m o nieskończonym
2
rozmiarze). Błąd wzgldędny wynosi:
x rd(x) m2C m2C 1
t
1 t t
= 2 = 2
x m2C 2
rd(x) = x(1 µ), |µ| 2t
Liczba cyfr w mantysie określa dokładność reprezentacji, a liczba cyfr cechy -
zakres. Cmax = Cmin = 2d, gdzie d to liczba cyfr cechy, zatem reprezentoeane
sÄ… liczby:
1
Cmax
min
2C |x| 2
2
Błąd niedomiaru pojawia się gdy C < Cmin, a błąd nadmiaru gdy C > Cmax.
1.3 Reprezentacja stałopozycyjna
Liczba reprezentowana jest na stałej ilości bitów, w wartość jej reprezentacji
wynosi w systemie dziesiętnym:
d
1
i
a = s ei 10
i=0
gdzie s " { 1, 1}, d to ilość cyfr i e wartość i-tej cyfry. W systemie dwójkowym
i
wartość reprezentacji wynosi:
d
1
a = s ei i2
i=0
2 Uwarunkowanie zadania
Uwarunkowanie zadania nie zależy od metod jego rozwiązania. Jeśli mamy
dane: d0, d1, . . . , dn " R i zamiast ich wartości dysponujemy ich reprezenta-
t
cjami rd(di) = di(1 + µi), gdzie |µi| 2 , to ta niewielka zmiana może spo-
wodować duże zmiany względne rozwiązania. Jeżeli niewielkie względne zmiany
danych powodują duże względne zmiany jego rozwiązania to zadanie nazywamy
z
le uwarunkowanym. Wielkości charakteryzujące wpływ zaburzeń danych na
zaburzenia rozwiÄ…zania nazywamy wskaz
nikami uwarunkowania zadania.
2
3 Stabilność i poprawność numeryczna
3.1 Stabliność numeryczna algorytmów
Definicja 1. Stabilność numeryczna algorytmów
Mamy dane:
D - zbiór danych
a = (a1, a2, . . . , ak) - wektor danych
W = (W1, W2, . . . , Wk) - wektor wyniku
W = W (a)  algorytm idealny
2 2
W = W N(a, µ)  algorytm numeryczny (W - wynik numeryczny)
DN(µ) - zbiór danych, dla których okreÅ›lony jest algorytm W N(a, µ).
Mówimy, że algorytm obliczeniowy jest numerycznie stabilny, jeśli:
"a0 " D "µ0 : "µ < µ0 a0 " DN(µ) '" lim W N(a0, µ) = W (a0)
µ0
czyli algorytm jest numerycznie stabilny wtedy, gdy zwiększając dokładność obli-
czeń można wyznaczyć (z dowolną dokładnością) dowolne istniejące rozwiązanie
zadania.
3.2 Poprawność numeryczna algorytmów
Algorytmy poprawne numerycznie są podzbiorem algorytmów numerycznie sta-
bilnych.
Definicja 2. Algorytm numerycznie poprawny.
Za numerycznie najwyższej jakości uznajemy takie algorytmy, dla których ob-
liczone rozwiązanie jest nieco zaburzonym rozwiązaniem (dokładnym) zadania
o nieco zaburzonych danych. Algorytmy spełniające ten postulat nazywamy nu-
merycznie poprawnymi.
3.3 Złożoność obliczeniowa algorytmów
Definicja 3. Złożoność obliczeniowa algorytmu.
Złożoność obliczeniowa to funkcja określająca maksymalną ilość operacji, które
należy wykonać, aby obliczyć wynik za pomocą danego algorytmu.
4 Miejsca zerowe funkcji nieliniowych
4.1 Metoda bisekcji (połowienia przedziałów)
Warunkami zbieżności tej metody jest to, aby funkcja posiadała dokładnie jeden
pierwiastek w przedziale izolacji, była w num ciągła, a na jego końcach miała
przeciwne znaki.
Algorytm wyznaczania kolejnego przybliżenia:
an 1+bn 1
1. Wyznaczamy środek przedziału: mn = , dla n = 1, 2, . . ..
2
2. Jeśli mn jest wystarczającym przybliżeniem, koniec.
3. Jeśli f(an 1) f(m < 0 to an = an 1, a bn = mn.
)
n
4. Jeśli f(bn 1) f(m < 0 to an = mn, a bn = bn 1.
)
n
3
4.2 Metoda Newtona (stycznych)
Warunkami zbieżności tej metody jest to, aby funkcja posiadała dokładnie jeden
pierwiastek w przedziale izolacji, była w num ciągła, a na jego końcach miała
przeciwne znaki. Funkcja musi być klasy C2 oraz pierwsza i druga pochodna
muszą mieć stały znak w przedziale.
Algorytm wyznaczania kolejnego przybliżenia:
2 2
1. Zerowe przybliżenie wyznaczene jest z warunku: f(a) f(a) > 0. Jeśli
jest on prawdziwy x0 = a, jeśli nie to x0 = b.
f(xn 1
2. Kolejne przybliżenie wyznaczamy ze wzoru: xn = xn 1 f2 (xn ) .
)
1
3. Jeśli przybliżenie nas satysfakcjonuje kończymy obliczenia.
4.2.1 Szybkość zbieżności metody
Definicja 4. Szybkość zbieżności metody Newtona.
Niech x0, x1, x2, . . . bÄ™dzie ciÄ…giem zbieżnym do ; µ = xn .
n
Jeśli istnieją takie dwie liczby p, c ,gdzie c = 0, że:

µn+1
lim sup = c
ninf µp
n
p nazywamy wykładnikiem zbieżności ciągu
c  stałą asymptotyczną błędu.
Dla p = 1, 2, 3 mówimy o zbieżności odpowiednio, liniowej, kwadratowej i sze-
ściennej.
4.3 Metoda siecznych
Warunki zbieżności takie jak w metodzie Newtona.
Algorytm wyznaczania kolejnych przyliżeń:
1. Pierwsze dwa przybliżenia wyzbaczamy z krańców przedziału: x0 = a i
x1 = b.
1 n 1 2
2. Kolejne przybliżenie wyznaczamy ze wzoru: xn = xn 1 f(xn ) (x xn ) ,
f(xn 1) f(xn 2)
dla n = 2, 3, 4, . . ..
3. Jeśli przybliżenie nas satysfakcjonuje kończymy obliczenia.
4.4 Metoda regula falsi (tylko na ćwiczeniach)
Warunki zbieżności takie jak w metodzie Newtona.
Algorytm wyznaczania kolejnych przyliżeń:
a f(b) b f(a)
1. Pierwsze dwa przybliżenia wyzbaczamy ze wzoru: x0 =
f(b) f(a)
2. Wyznaczamy punkt stały metody z warunku: f(a)f(x) < 0, jeśli jest on
prawdziwy to c = a, w przeciwnym przypadku c = b.
xn 1f(c) cf(xn 1)
3. Kolejne przybliżenie wyznaczamy ze wzoru: xn = , dla
f(c) f(xn 1)
n = 1, 2, 3, 4, . . ..
4. Jeśli przybliżenie nas satysfakcjonuje kończymy obliczenia.
4
5 Interpolacja
W przedziale [a, b] danych jest n + 1 różnych punktów x0, x1, ..., xn(węzły in-
terpolacji) oraz wartości funkcji y = f(x) w tych punktach f(x0) = y0, f(x1) =
y1, ..., f(xn) = yn. Zadaniem interpolacji jest znalezienie F (x) , która w węzłach
interpolacyjnych ma te same wartości co f(x) i przybliża f(x) w punktach po-
średnich.
Twierdzenie 1. O istnieniu dokładnie jednego wielomianu interpolacyjnego.
Istnieje dokładnie jeden wielomian interpolacyjny stopnia co najwyżej n, (n >
0), który w punktach x0, x1, ..., xn przyjmuje wartości y0, y1, ..., yn.
Dowód.
Wn = a0 + a1x + a2x2 + ... + anxn (1)
korzystając z n + 1 (par) wartości (xi, y(xi)) mamy układ równań:
Å„Å‚
n
a0 + a1x0 + a2x2 + x+ = y0
a
ôÅ‚
n
0 0
ôÅ‚
òÅ‚
n
a0 + a1x1 + a2x2 + x+ = y1
a
n
1 1
(2)
ôÅ‚
ôÅ‚
ół
n
a0 + a1xn + a2xn + x+ = yn
a
n
2 n
z n + 1 niewiadomymi a0, a1, ..., an
1 x0 x2 ... xn
0 0
1 x1 x2 ... xn
1 1
D = detA = (3)
... ... ... ... ...
1 xn x2 ... xn
n n
Wyznacznik Vandermande a = 0, gdy xi = xj dla i = j

n

1
ai = f(xj)Dij (4)
D
j=0
4 dopełnienie algebraiczne elementów i tej kolumny, j tego wiersza
5.1 Interpolacja Lagrange a
Wzory ostateczne:
(x x1) (x x (x x0) (x x
) )
2 2
Ln(x) = y0 (x0 x1) (x x2) ... (xn)x + y1 (x1 x0) (x x2) ... (xn)x + . . .
... (xxn) ... (xxn)
0 0 1 1
(x x0) (x x
)
1 1
+yn (xn x0) (x x1) ... (xn x )
... (x xn 1)
n n
Oznaczamy: Wn+1(x) = (x x) (x )x . . . (xn) x
0 1

n+1czynnikow
n

Wn+1(x)
Ln(x) = yj
Wn1 (x)
j=0
(x x)
j
(x xj)
x=xj
jest to wzór interpolacyjny Lagrange a.
5
5.2 Interpolacja Newtona
Aby biegle korzystaż z wzoru interpolacyjnego Newtona, przypomn3
my sobie
jak wygląda iloraz różnicowy:
f[xi+1, xi+1] f[x, xi]
i
f[xi, xi+1] =
xi+1 x
i
W przypadku ogólnym korzystać będziemy z jego zfediniowanej rekurencyj-
nie postaci:
f[xi+1, . . . , xi+j+1] f[x, . . . , xi+j]
i
f[xi, . . . , xi+j+1] =
xi+j+1 x
i
Wzór interpolacyjny Newtona:
n i
1
w(x) = ai (x x) (5)
j
i=0 j=0
gdzie ai oznacza ity iloraz różnicowy (patrz poniższy algorytm).
Algorytm "wizualnego"ułatwienia sobie liczenia ilorazów różnicowych wyż-
szych rzędów:
1. Stwórz tabelę o 2 kolumnach i tylu wierszach ile masz punktów interpolacji
2. Pierwszą kolumnę wypełn3
wartościami x tych punktów a drugą f(x)
3. Każdą kolejną kolumnę utwórz poprzez obliczenie ilorazu różnicowego z
 2 poprzednich oraz odpowiadających im wartości x z pierwszej kolumny
4. InteresujÄ…ce nas do wzoru interpolujÄ…cego newtona ilorazy sÄ… pierwszymi
elementami w poszczególnych kolumnach
Poniżej przykżadowa tabelka ilorazów różnicowych:
x f(x) a1 a2 a3 a4
2 2 2
0 1 1
3 3 9
2
2 3 1 2
3
2
3 2 3
3
4 5 1
6 7
Rozpisany wzór Newtona dla funkcji przedstawionej w powyższej tabeli:
2 2 2
w(x) = 1 + 1(x 0) (x 0)(x 2) + (x 0)(x 2)(x 3) (x 0)(x 2)(x 3)(x 4) =
3 3 9
2 8 88 35
x4 + x3 x2 + x + 1
9 3 9 3
5.3 Interpolacja Hermite a
Interpolacja Hermita poza wartościami funkcji korzysta również z pierwszej po-
chodnej interpolowanej funkcji, umożliwia zatem uzyskanie dokładniejszego od-
wzorowania.
n n

H2n+1(x) = f(xj)Hn,jx + f2 (xj)Hdn,jx
j=0 j=0
6
gdzie
Hn,j(x) = (1 2(x x (xj))L2 (x)
)L2
j
n,j n,j
oraz
Hdn,j(x) = (x x)L2 (x)
j
n,j
a Ln,j jest oczywiście jtym współczynnikiem wielomianu Lagrange a stopnia
n.
Twierdzenie 2. Zadanie interpolacyjne Hermite a ma jednoznaczne rozwiÄ…za-
nie.
5.4 Efekt Rungego
Przy równoodległych węzłach w środku przedziału są małe błędy, ale na krań-
cach gwałtownie wzrastają i funkcja się  rozjeżdża .
Do optymalnego doboru węzłów interpolacji słyży wielomian Czybyszewa:
Tn(x) = cos(arccos(x))
1. Rekurencyjne obliczanie wielomianu:
T0(x) = 1
T1(x) = x = x T
(x)
0
Tn+1(x) = 2 x T T dla n 1
n n 1
2) Tn(x) jest równoważny pewnemu wielomianowi algebraicznemu stopnia
n określonego w przedziale [ 1, 1] 3) Współczynnik wiodący (przy najwyższej
potędze) Tn wynosi:
2n 1 dlan 1
1 dlan = 1
2k+1
4) tn ma n zer w [ 1, 1] jednokrotnych, rzeczywistych: x = cos ,
k
n 2
2k 1
dla k = 0, 1, ..., n 1 lub x = cos , dla k = 1, 2, ..., n
k
n 2
Przeskalowanie z [ 1, 1] do [a, b]: Dowolną f(t) dla t " [a, b] można prze-
1
skalować, aby była określona w [ 1, 1] i odwrotnie: t =1 (a + b) + (b a) x
2 2
t " [a, b] Ô! x " [ 1, 1] przeskalowanie nie zmienia bÅ‚Ä™du interpolacji.
1 2k+1
Optymalne węzły w [a, b]: tk = (b a) cos + (b + a) dla
2 n 2
k = 0, 1, ..., n 1
6 Interpolacja funkcjami sklejanymi
Funkcje sklejane brzmiÄ… strasznie ale sÄ… bardzo przyjemne, inna nazwa na nie
to z ang. spline, a spalszczamy je (kto wie czy poprawnie) jako śplajny"...
Ogólny pomysł jest następujący: Zamiast 1 funkcji interpolującej, chcemy
uzyskać tyle funkcji ile jest przedziałów interpolacji. Na przykład dla 3 węzłów
interpolujących, uzyskamy 2 splajny, które gdy je śkleimy"utworzą naszą funkcję
interpolujÄ…cÄ….
Zalety splajnów:
Szybciej się liczy spline niż zwyczajne interpolacje wielomianowe, powo-
dem jest stopień wielomianów na których pracujemy (przy spline jest za-
zwyczaj bardzo niski, popularna jest wartość 3)
SÄ… odporne na efekt Rungego (?!?!)
7
6.1 Interpolacja funkcjami sklejanymi 1 stopnia
Inaczej nazywama: liniowa. Liniowa dlatego że spline 1 stopnia, to nic innego
niż funkcja liniowa. A całą interpolację można przedstawić jako "połącz punkty
interpolacji prostymi odcinkami". Funkcja otrzymana w wyniku tej interpolacji
będzie bardzo kanciasta, ale również bardzo szybko ją obliczymy. Wzór na kty
spline wygląda następująco:
yk+1 y
k
sk(x) = yk + (x x )
k
x
k+1 k
Funkcję rzeczywistą nazywamy spline stopnia m z węzłami s = a = x0, x1 N = x
b jeśli
1. w każdym przedziale (xi 1, xi)dlai = 0, 1, , N + 1 s jest wielomianem
stopnia nie wyższego niż m
2. s i jej pochodne rzędu 1, 2, , m 1 są ciągłe na całej osi rzeczywistej;
innymi słowy s " Cm 1
Definicja 5. Funkcję sklejaną stopnia 2m 1 nazywamy naturalną jeśli w
przedziałach ( ", x), (xN , ") jest wielomianem stopnia m 1
0
Twierdzenie 3. Interpolacja naturalną funkcją sklejaną jest ńajgładszą"funkcją
interpolujÄ…cÄ… punktu (x, yi)
6.2 Interpolacja funkcjami sklejanymi 3 stopnia
Lorem ipsum dolor sit amet, consectetur adipiscing elit. Suspendisse aliquam
augue vitae leo ultricies imperdiet. Morbi at nibh nunc, quis suscipit est. Class
aptent taciti sociosqu ad litora torquent per conubia nostra, per inceptos hime-
naeos. Donec sit amet gravida lectus. Nam tincidunt urna velit. Cras sagittis
enim gravida turpis volutpat vel tempus nisi imperdiet. Donec risus diam,
pellentesque id vulputate eu, dignissim non leo. Nullam aliquet odio eu enim
egestas vehicula. Integer commodo porta tortor sit amet accumsan. Curabitur
ac enim id metus malesuada tempor eu id turpis. Aenean magna nibh, ullam-
corper in ornare at, pellentesque nec sapien. Proin tincidunt sapien scelerisque
nulla viverra molestie. Aliquam et massa in sem auctor molestie eget sit amet
mauris. Morbi mi velit, adipiscing sed condimentum at, iaculis ut neque. In
non ultrices diam.
8
7 Licencja
Ten utwór objęty jest licencją Creative Commons:
Uznanie autorstwa
Użycie niekomercyjne
Na tych samych warunkach
3.0 Polska.
Aby zobaczyć kopię niniejszej licencji przejdz
na stronÄ™ http://creativecommons.org/licenses/by-
nc-sa/3.0/pl/ lub napisz do Creative Commons, 171 Second Street, Suite 300,
San Francisco, California 94105, USA.
9