Metody numeryczne I (wykład do wyboru) Krzysztof Tabisz Instytut Matematyczny Uniwersytet Wrocławski 2002/2003 2002 by Krzysztof Tabisz 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Spis treści Wstęp 1 1 Wykład 1 2 2 Wykład 2 i 3 5 3 Wykład 4 9 4 Wykład 5 13 5 Wykład 6 18 6 Wykład 7 24 7 Wykład 8 i 9 29 8 Wykład 10 39 9 Wykład 11 43 10 Wykład 12 53 11 Wykład 13 57 12 Wykład 14 61 13 Metoda gradientów 63 A Uzupełenienia 64 A.1 Analiza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 A.2 Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 A.2.1 Struktura form dwuliniowych . . . . . . . . . . . . . . . . 64 A.3 Analiza funkcjonalna . . . . . . . . . . . . . . . . . . . . . . . . . 68 Literatura 69 Indeks 70 iii Krzysztof Tabisz 25 marca 2003, 8:39am iv 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Wstęp Notatki do wykładu, stan na: 31.I.2003r 0.Met1Int 1 Krzysztof Tabisz 25 marca 2003, 8:39am 1 Wykład 1 Wykład 1 Postawowe fakty z ANALIZY: " A " R jest ograniczony, gdy "M "a"A |a| d"M; " Niech A " R ograniczony i A =. Wtedy ą =sup A, gdy
"a"A a d" ą, ("a"A a d" ) =! ą d" , " aksjomat ciągłości R: dla każdego niepustego ograniczonego A " R istnieje sup A; " aksjomat ciągłości R jest równoważny faktowi: każdy ciąg monotoniczny ograniczony jest zbieżny; " każdy odcinek [a, b] " R ograniczony jest zwarty, tzn. z każdego ciągu jego elementów można wybrać podciąg zbieżny (granica wybranego podciągu ma należeć do [a, b]!); " f : R - R ma granicę w punkcie x0, tzn. lim f(x) =g gdy xx0 (Cauchy) ">0"()>0"x 0 < |x - x0| <() =!|f(x) - f(x0)| <,
n n+ " n " " f : R - R jest ciągła w x0, gdy lim f(x) =f(x0); xx0 " trzy podstawowe twierdzenia o funkcjach f : [a, b] - R ciągłych: y = f(x) jest ograniczona na [a, b], y = f(x) osiąga swoje kresy na [a, b], przyjmuje wszystkie wartości pośrednie między minx"[a,b] f(x) a maxx"[a,b] f(x). Innymi słowy f([a, b]) = [minx"[a,b] f(x), maxx"[a,b] f(x)]. f(x)-f(x0) " określenie pochodnej y = f(x) w x0: f (x0) = lim ; xx0 x-x0 " warunek konieczny ekstremum y = f(x) przyjętego w x0: f (x0) =0; " trzy twierdzenia o wartości średniej dla funkcji f : [a, b] - R ciągłych i różniczkowalnych na (a, b): (Rolle) jeżeli dodatkowo f(a) =f(b) =0, to ""(a,b)f () =0, (Lagrange) ""(a,b) f(x)-f(x0) = f (), x-x0 (wzór Taylor a) jeżeli dodatkowo f : (a, b) - R jest (n +1) razy różniczkowalna, to wtedy n
f(k)(x0) "x,x "(a,b)""(x ,x) f(x) = (x - x0)k + En(x), 0 0 k! k=0 gdzie f(n+1)() En(x) = (x - x0)n+1 (n +1)! 21.L1-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Przykład strzelanie z armaty na odległość d: " model matematyczny:
y(0) = 0, y (t) =-g, , y (0) = V0 sin 2V0 sin " rozwiązanie teoretyczne: T = , g 2V0 sin " czas lotu pocisku: T = , g 2 2V0 cos sin " przebyta droga przez pocisk: = d, g " określenie nieliniowego równania: 2V02 cos sin f() a" - d =0, g " dyskusja warunków istnienia rozwiązania, " założenie o funkcji y = f() ciągłości, " algorytm znajdywania rozwiązania nieliniowego równania przez połowienie przedziału: while(abs(b-a)>eps) { c = (b+a)/2; if(c==a||c==b) return; fc = f(c); if(fc==0) { a = b = c; fa = fb = fc; return; } if(sign(fc)!=sign(fb)) { a = c; fa = fc; } else { b = c; fb = fc; } } Zadania na ćwiczenia: 1.L1-met 3 Krzysztof Tabisz 25 marca 2003, 8:39am 1. Pokazć, że |x2 - 4| < dla 0 < |x - 2| <(5 + )-1 oraz korzystając tylko z definicji udowodnić: lim x2 =4. x2 2. Pokazć, że lim (4x +2) =6 wprost z -- definicji. x1 3. Dla f(x) =3 - 2x + x2 i odcinka [a, b] =[1, 3] wyznaczyć z twierdzenia o wartości średniej. 4. Wyznaczyć szereg Taylor a funkcji f(x) =cosh x dla x0 =0. 1 5. Pokazać, że funkcja f(x) =x sin z f(0) = 0 jest ciągła w 0 ale nie jest x różniczkowalna w 0. Zadania na Laboratoria: 1. Napisać w Pascal u procedurę rozwiązywania równania f(x) =0 metodą przez połowienie przedziałui. 2. Każdy student dla innego równania podanego przez prowadzących powi- nien zastosować procedurę do: " empirycznego określenia maksymalnęj sensownęj dokładność względ- nęj i bezwzględnęj jaką można uzyskać przy wykonywaniu obliczeń przy zastosowaniu typów rzeczywistychReal, Single,Double, Extended,Comp; " wyliczyć ilość iteracji jaką trzeba wykonać aby uzyskać zadaną do- kładność; " wykonać obliczenia dla różnych danych wejściowych; " uzyskane wyniki oddać w formie sprawozdania; " termin: do 25.X.00r. Literatura: 1. G.W. Stewart, Afternotes on Numerical Analysis, SIAM, Philadelphia 1996; 2. D. Kincaid & W. Cheney, Numerical Analysis, Brooks 1996; 3. Anthony Ralston, Wstęp do analizy numerycznej, PWN, Warszawa 1983. 41.L1-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 2 Wykład 2 i 3 Uzupełnienie wykładu 2 Twierdzenie (Metoda Newton a) Niech f : R - R będzie dwukrotnie różnicz- kowalna w sposób ciągły oraz r będzie jednokrotnym zerem równania f(x) =0. Wtedy istnieje otoczenie r i stała C takie, że jeżeli przyjmiemy punkt startowy x0 z tego otoczenia dla ciągu iteracji f(xn) xn+1 = xn - , n e" 0 f (xn) to lim xn = r oraz n " |xn+1 - r| d"C(xn - r)2, n e" 0. Dowód: r jest jednokrotnym zerem funkcji y = f(x), gdy f(r) =0 oraz f (r) =0. Niech
en = xn - r, oznacza błąd przybliżenia. Wtedy mamy f(xn) en+1 = xn+1 - r = xn - - r (2.1) f (xn) f(xn) enf (xn) - f(xn) = en - = f (xn) f (xn) Ze wzoru Taylor a mamy 1 0 =f(r) =f(xn - en) =f(xn) - enf (xn) + e2 f (n), gdzie n " (r, xn). n 2 Skąd otrzymujemy 1 enf (xn) - f(xn) = f (n)e2 . n 2 Co wstawiając do (1) otrzymujemy 1 f (n) 1 f (r) en+1 = e2 H" e2 = Ce2 . (2.2) n n n 2 f (xn) 2 f (r) Dla C H" 1 i en H" 10-4 z (2) mamy: en+1 H" 10-8, en+2 H" 10-16 . . . . W ten sposób błąd en+1 jest w przybliżeniu równy e2 co oznacza, iż mamy do n czynienia z tzw. kwadratową zbieżnością. Przejdzmy teraz do formalnego dowodu zbieżności metody. Niech max |f (x)| 1 |x-r|d" def c() = , > 0. (2.3) 2 min |f (x)| |x-r|d" Ponieważ f (r) =0 oraz y = f (x) jest funkcją ciągłą, to lim c() =0. Więc
0 istnieje taka , iż = c() < 1. Teraz, jeżeli zaczniemy proces iteracji z punktu 2.L3-met 5 Krzysztof Tabisz 25 marca 2003, 8:39am x0, spełniającego |x0 - r| d", to wtedy |e0| d" i |0 - r| d" i z określenia c() otrzymujemy
f (0) 1
d" c().
2 f (x0) Tym samym z (3) mamy |x1 - r| = |e1| d"e2c() =|e0||e0|c() d"|e0|c() =|e0| <|e0| d". 0 To pokazuje, że kolejny punkt x1 jest w otoczeniu r i iteracja może być kon- tynuowana |e1| d" |e0| |e2| d" |e1| d" 2|e0| |e3| d" |e2| d" 3|e0| . . . Ogólnie mamy |en| d" n|e0|. Poniewż 0 d" <1, to mamy lim n =0 i lim en =0. n " n "
Jako ćwiczenie udowodnić: Twierdzenie Jeśli f " C2(R) jest malejąca, wypukła i istnieje r rozwiązanie równania f(x) =0, to rozwiązanie to jest jedyne i ciąg iteracji Newton a star- tujący z dowolnego x0 " R jest do niego zbieżny. Wykład 3 Rozwiązywanie równania f(x) =0 metodą siecznych. W metodzie Newton a f(xn) xn+1 = xn - , n =0, 1, 2, . . . f (xn) przybliżamy f(xn) - f(xn-1) f (xn) H" , n =1, 2, . . . xn - xn-1 otrzymuje wtedy (metoda siecznych)
f(u)-f(v) (u, v) = f(u) # u - , dla u = v f (u) to wtedy ciąg iteracyjny z metody siecznych możemy zapisać w postaci xn+1 = (xn, xn-1), gdzie określa iteracje. Szukamy punktu stałego r, r = (r, r) 62.L3-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 tej iteracji. Ponieważ prawa strona zależy od dwóch kolejnych punktów to me- todę tą nazywa się też dwupunktową. Analiza zbieżności Będziemy wzorować się na dowodzie Twierdzenia w metodzie Newton a. Niech, tak jak poprzednio en = xn - r. Wtedy f(xn)xn-1 - f(xn-1)xn f(xn)en-1 - f(xn-1)en en+1 = xn+1 - r = - r = f(xn) - f(xn-1) f(xn) - f(xn-1) Dzieląc przez enen-1 i wstawiajćxn-xn-1 =1 otrzymujemy xn-xn-1 f(xn) f(xn-1) xn - xn-1 en - en-1 en+1 = enen-1 (2.4) f(xn) - f(xn-1) xn - xn-1 Z wzoru Taylor a mamy 1 f(xn) =f(r + en) =f(r) +enf (r) + e2 f (r) +O(e3 ). n n 2 Ponieważ f(r) =0, to f(xn) 1 = f (r) + enf (r) +O(e2 ) n en 2 oraz dla wskaznika n - 1 f(xn-1) 1 = f (r) + en-1f (r) +O(e2 ). n-1 en-1 2 Ostatecznie f(xn) f(xn-1) 1 - = (en - en-1)f (r) +O(e2 ). n-1 en en-1 2 czywiście xn - xn-1 = en - en-1. Możemy teraz f(xn) f(xn-1) - 1 en en-1 H" f (r). xn - xn-1 2 Pierwszy czynnik w wyrażeniu (4) xn - xn-1 1 H" . f(xn) - f(xn-1) f (r) Pokazaliśmy w ten sposób, iż 1 f (r) en+1 H" enen-1 = Cenen-1. (2.5) 2 f (r) Otrzymane równanie (5) jest analogiczne do równania (2) z metody Newton a. Aby oszacować rząd zbieżności metody załóżmy, że |en+1| <"A|en|ą, gdzie A>0 stała. 2.L3-met 7 Krzysztof Tabisz 25 marca 2003, 8:39am Tzn. ą rząd zbieżności metody oznacza, iż |en+1| lim =1. n " A|en|ą Wstawiając 1 ą |en| <"A|en-1|ą i |en-1| <"(A-1|en|) do (5), otrzymujemy 1 1 ą ą A1+ |C|-1 <"|en|1-ą+ . 1 Ponieważ prawa strona zależy od en, lim en = 1, to 1 - ą + = 0, więc ą n " " 1+ 5 ą = H" 1.62 . W ten sposób otrzymaliśmy, iż szybkość zbieżności metody 2 1 siecznych jest superliniowa tzn. lepsza niż liniowa. Z równania 1 + = ą ą otrzymujemy 0.62 1
1 f (r) 1 1+
ą A = |C| ą = |C| = |C|ą-1 = |C|0.62 = . 2f (r)
Z tak wyliczonym A otrzymujemy " 1+ 5 2 |en+1| H"A|en| . Zadania na Ćwiczenia 1 1. Dla równania =0 znalezć najmniejszy dodatni punkt startowy me- tan x tody Newton a, dla którego metoda ta jest jeszcze zbieżna. 2. Metoda Steffensen a. Niech f(xn) f(x + f(x)) - f(x) xn+1 = xn - , gdzie g(x) = . g(xn) f(x) Przy odpowiednich założeniach, wykazać kwadratową zbieżność metody. 3. W metodzie Halley a rozwiązywania równaia f(x) =0 ciąg iteracji okre- śla się w następujący sposób:
fnfn xn+1 = xn - , gdzie fn = f(xn).
fnfn
(fn)2 - 2 f " Przez podstawienie f - wyprowadzić powyższy wzór z metody New- f ton a. Zadania na Laboratoria 1. Metodą Newton a znalezć rozwiązania równania x =tan x najblizsze punk- tów 4.5 i 7.7. tan x 2. Znalezć najmniejsze dodatnie x0 dla którego funkcja f(x) = osiąga x-2 swoje minium. Zastosować metodę Newton a do f . 3. Równanie 2x4 +24x3 +61x2 -16x+1 = 0 ma dwa pierwiastki w otoczeniu 0.1. Posługując się metodą Newton a znalezć je. 82.L3-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 3 Wykład 4 Wykład 4 Metody przez połowienie przedziału , Newton a (stycznych) oraz siecznych wy- znaczenia pierwiastka r " R będącego rozwiązaniem równania f(r) =0 (3.6) są przykładami metod iteracyjnych w których r jest ganicą ciągu xn+k = (xn, xn+1, . . . , xn+(k-1)), n, k " N, k ustalone. (3.7) Dla metody Newton a mamy: f(u) (u) =u - , f (u) a dla siecznych ż# vf(u)-uf (v)
# , dla u = v f(u)-f(v) (u, v) = # f(u) u - , dla u = v f (u) Pierwiastek r równania (3.6) jest granicą ciągu r = lim xn, n " i dodatkowo jest spełniona równość wynikająca z (3.7) r = (r, . . . , r). Definicja1 rząd metody (iteracyjnej) Niech en = xn - r, będzie błędem n tej iteracji. Załóżmy dodatkowo, że ciąg przybliżeń (3.7) rów- nania (3.6) jest zbieżny lim xn = r. n " Mówimy, że w punkcie r metoda jest rzędu p, jeśli istnieje liczba rzeczywista p e" 1 taka, że |xn+1 - r| |en+1| lim = lim = C =0.
n " - r|p |en|p n " |xn Stałą C =0 nazywamy stałą asymptotyczną błędu, jest ona zależna od funkcji
y = f(x). Dla metody Newton a wykazaliśmy, iż p =2 a przy metodzie siecznych otrzy- " 1+ 5 maliśmy p = =1.618 . . . . 2 1 A. Ralston, Wstęp do analizy numerycznej, PWN, Warszawa 1983, ss. 298 299. 3.L4-met 9 Krzysztof Tabisz 25 marca 2003, 8:39am Arytmetyka na liczbach zmiennopozycyjnych. 1. Przedstawienie liczb rzeczywistych w komputerze: " typy rzeczywiste w TP ver. 6.0 dokładność zajętość typ zakres wcyfrach wbajtach Real 2.9 10-39 . . . 1.7 1038 11 - 12 6 Single 1.5 10-45 . . . 3.4 1038 7 - 84 Double 5.0 10-324 . . . 1.7 10308 15 - 16 8 Extended 3.4 10-4932 . . . 1.1 104932 19 - 20 10 Comp -263 +1 . . . 263 - 119 - 20 8 " typy rzeczywiste w TC ver. 3.1 dokładność zajętość typ zakres wcyfrach wbajtach float 3.4 10-38 . . . 3.4 1038 74 double 1.7 10-308 . . . 1.7 10308 15 8 long double 3.4 10-4932 . . . 1.1 104932 19 10 " sposób przedstawienia liczby rzeczywistej x = m d gdzie m mantysa, podstawa i d cecha liczby; " postać znormalizowana liczby, tzn. w przedstawieniu x = m d jest 0.1 d"|m| < 1 dla x =0.
" standard IEEE przedstawienia liczb zmniennopozycyjnych single precision 32 bitowe z 1 bitem na znak, 8 bitową cechą oraz 23 bitową mantysą (ponieważ zapis liczby jest znormalizo- wany, to pierwszą cyfrą jest 1 bo bazą jest 2, która nie jest pamiętana), double precision 64 bitowe z 1 bitem na znak, 11 bitową cechą i 52 bitową mantysą. 2. Nadmiar i niedomiar przedstawienia liczb rzeczywistych " przedstawienie zera: -0, +0, " wartości: - ", + ", 0 " wynik nieokreślony (np. dla ): NaN. 0 3. Rodzaje zaokrągleń 10 3.L4-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 " obcięcie ( chopping , truncation ), " zaokrąglenie ( rounding ), " błąd zaokrąglenia dla a = 0 liczby o podstawie 2 ( unit roundoff
error lub machine epsilon ):
|b - a| 2-n, zaokrąglenie = gdzie n liczba cyfr w mantysie, |a| 2-n+1, obcięcie b - a def fl(a) = a(1 + ), gdzie = , a =0,
a " wyliczenie (problematyczne !) maksymalnego błędy zaokrągleń: x = 1; while(1+x!=1) x = x/2; 4. Działania na liczbach przybliżonych " arytmetyka fl(x(y + z)) = [fl(y + z)](1 + 1) = [x(y + z)(1 + 2)](1 + 1) = x(y + z)(1 + 2 + 1 + 21) H" x(y + z)(1 + 1 + 2) = x(y + z)(1 + 3) 5. Analiza błędów zaokrągleń (przykład) Twierdzenie Niech x0, x1, . . . , xn będą dodatnimi liczbami rzeczywistymi, które są przedstawione w komputerze z roundoff error . Wtedy błąd n
względny przedstawienia xi wynosi (1 + )n - 1 H" n. i=0 Zadania na Ćwiczenia do metody siecznych 1. W metodzie siecznych wykazać, że jeżeli xn q gdy n "i f (q) =0,
to wtedy f(q) =0. 2. Korzystając ze wzoru Taylor a uzasadnić następujące przybliżenie na f (x): k2f(x + h) - h2f(x + k) +(h2 - k2)f(x) f (x) H" . (k - h)kh 3. Metodę siecznych zastosować do funkcji f(x) =x2 - 2 z x0 =0 i x1 =1. Wyliczyć x2. 4. Jakie jest x2 z metodysiecznych, gdyx0 =1, x1 =2, f(x0) =2 oraz f(x1) = 1.5 ? 3.L4-met 11 Krzysztof Tabisz 25 marca 2003, 8:39am 5. Równość asymptotyczną dwóch ciągów określamy w następujący sposób: xn xn <" yn !! lim =1. n " yn Udowodnij, że jeżeli xn <" yn, un <" vn i c =0, to wtedy
c a) cxn <" cyn, b) xc <" yn, c) xnun <" ynvn, n d) yn <" un =! xn <" vn, e) yn <" xn. Zadania na Laboratoria 1. Napisać procedurę rozwiązywania równania f(x) =0 metodą siecznych oraz przetestować ją na funkcjach: x a) sin - 1; b) ex - tan x; c) x3 - 12x2 +3x +1. 2 2. Która z metod: połowienie przedziału , Newton a, siecznych jest naj- właściwsza (podać kryterium !) do znalezienia rozwiązania równania f(x) =0: a) x20 - 1 na [0, 10]; b) tan x - 30x na [1, 1.57]; c) x2 - (1 - x)10 na [0, 1]; d) x3 +10-4 na [-0.75, 0.5]; e) x19 +10-4 na [-0.75, 0.5]; f) x5 na [-1, 10]; 2 g) x9 na [-1, 10]; h) xe-x na [-1, 4]. Zadania na Ćwiczenia z przybliżeń komputerowych Przyjmujemy, że liczby są 32 bitowe ( single precision ) i mają przedstawienie zgodne ze standardem IEEE. 3 1. Przedstawić w postaci 12.a1a2 . . . a24. Jaki jest roundoff error tego 5 przedstawienia?
2 2. Liczbę 1 - 2-24 zapisać w postaci single precision . 3 3. Ile jest liczb single precision w komputerze? 4. Niech x1, x2, . . . , xn będą dodatnimi liczbami single precision oraz niech " Sn = x1 + x2 + . . . + xn będzie sumą matematyczną a Sn oznacza sumę komputerową. Udowodnić, że jeśli dla każdego i : xi+1 e" 2-24Si, to wtedy " |Sn - Sn| n - 1 d" . Sn 224 5. Uzasadnić nierówność
x - x" 1
d" , x =0
x 1+224 gdzie x" jest zapisem w postaci single precision liczby x. 12 3.L4-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 4 Wykład 5 Wykład 5 Podstawowe fakty z ALGEBRY LINIOWEJ2 Macierzą A " Rnm o n wierszach i m kolumnach nazywamy układ nm liczb: Ą# ń# a11 a12 . . . a1m ó# a21 a22 . . . a2m Ą# ó# Ą# A = ó# Ą# . . . . . . . . Ł# Ś# . . . . an1 an2 . . . anm oraz przyjmujemy, iż aij def (A)ij. Wektory x " Rm, b " Rn zapisujemy w = postaci Ą# ń# Ą# ń# x1 b1 ó# Ą# ó# Ą# x2 b2 ó# Ą# ó# Ą# x = ó# Ą# ó# Ą# . , b = . . . . Ł# Ś# Ł# Ś# . . xm bn Przy tak przyjętych oznaczeniach układ równań ż# a11 x1 + a12 x2 + . . . + a1m xm = b1 # # # a21 x1 + a22 x2 + . . . + a2m xm = b2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . # # # an1 x1 + an2 x2 + . . . + anm xm = bn (4.8) możemy zapisać w postaci: A x = b. Definicja Dwa układy równań A x = b, B x = d są równoważne, gdy mają dokładnie te same rozwiązania. Twierdzenie Niech Ei oznacza i te równanie układu (5.10). Wtedy, jeżeli jeden układ równań można otrzymać z drugiego przez skończony ciąg operacji: " przestawienie dwóch równań Ei ! Ej, " pomnożenie równania przez niezerową liczbę Ei - Ei, " dodanie pomnożone przez liczbę inne równania Ei + Ej - Ei, to oba te układy są równoważne. Definicja: Macierz transponowana AT jest określona wzorem: (AT )ij = aji, gdzie aij =(A)ij. Definicja: Macierz A jest symetryczna, gdy: AT = A. Operacje na macierzach 2 B. Gleichgewicht, Algebra, PWN, Warszawa 1983 4.L5-met 13 Krzysztof Tabisz 25 marca 2003, 8:39am " dodawanie A + B: (A + B)ij =(A)ij +(B)ij, " mnożenie przez liczbę A, gdzie " R: (A)ij = (A)ij, " iloczyn dwóch macierzy A " Rnp i B " Rpm: p
(AB)ij = aik bkj. k=1 Uwaga: Istnieją macierze A, B takie, że AB = BA. Oznacza to, że
iloczyn macierzy nie jest przemienny! Własności operacji macierzowych: " (A + B) +C = A +(B + C) łączność dodawania, " (AB)C = A(BC) łączność mnożenia, " A(B + C) =AB + AC rozdzielność mnożenia względem dodawania, " (A)T = AT , " (A + B)T = AT + BT , " (AB)T = BT AT . Definicja: Macierz identycznościowa (jednostkowa) I " Rnn Ą# ń# 1 0 0 . . . 0 ó# Ą# 0 1 0 . . . 0 ó# Ą# def ó# Ą# 0 0 1 . . . 0 I = ó# Ą# ó# Ą# . . . . . . . . . . Ł# Ś# . . . . . 0 0 0 . . . 1 Własność macierzy identycznościowej AI = IA = A, dla A " Rnn. Definicja: Lewostronna A " Rnm i prawostronna B " Rmn macierz odwrotna jest określona wzorem AB = I. Twierdzenie Jeżeli A, B " Rnn i AB = I, to BA = I. Określenie Jeżeli dla A " Rnn (macierzy kwadratowej) istnieje B " Rnn: AB = I, tomacierz A nazywamy odwracalną (nieosobliwą) oraz ozna- czamy A-1 ozn. B. = Definicja: Przekształcenie f : Rn - Rm nazywamy liniowym, gdy spełnia następujący warunek: f(x + y) =f(x) +f(y), dla każdych , " R & x, y " Rn. 14 4.L5-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Związek przekształceń liniowych z macierzami Jeżeli f : Rn - Rm jest przekształceniem liniowym oraz (e1, . . . , en) jest bazą Rn a (e , . . . , e ) bazą Rm, to macierz 1 m (A)ij =(f(ei))j, tzn. f(ei) =ai1e + . . . + aime , 1 m ma własność: Ax = f(x) dla x " Rn. I odwrotnie, jeżeli A " Rnm jest macierzą, to def f(x) = Ax jest przekształceniem liniowym f : Rn - Rm wyznaczonym przez A. Definicja: Jeżeli A " Rnn oraz Ax = x, gdzie " R, x " Rn i x = 0
to nazywamy wartością własną macierzy A a x wektorem własnym tej ma- cierzy. Definicja: Wyznacznika macierzy A " Rnn:
def det(A) = sgn() a1(1) a2(2) . . . an(n) "Sn gdzie Sn zbiór wszystkich permutacji zbioru {1, 2, . . . , n} a sgn() jest znakiem permutacji . Twierdzenie: (Cauchy ego) Dla A, B " Rnn zachodzi wzór: det(AB) =det(A)det(B) Twierdzenie Rozwinięcie Laplace a wyznacznika macierzy A (może tez być defini- cją indukcyjną wyznacznika): n
det(A) = aijAij, j=1 gdzie Aij jest dopełnieniem algebraicznym wyrazu aij. Twierdzenie (Kroneckera Capellego) Układ równań (5.10) ma rozwiązanie wtedy i tylko wtedy, gdy r(A) =r(Ab) gdzie r(A) jest rzędem macierzy A a macierz Ab powstała z A przez dodanie kolumny wyrazów wolnych b. Rozwiązanie układu (5.10) dla n = m i det(A) =0 dane jest wzorem
(wzory Cramera): det(A ) j xj = , j =1, 2, . . . , n det(A) 4.L5-met 15 Krzysztof Tabisz 25 marca 2003, 8:39am gdzie macierz A jest utworzona z A przez zastąpienie j tej kolumny j kolumną wyrazów wolnych b. Definicja: Macierz A " Rnn nazywamy elementarną, gdy powstała z macierzy identycznościowej w wyniku następujących operacji: " zamiany dwóch wierszy (lub kolumn): As ! At; " pomnożenia wiersza (lub kolumny) przez różną od zera stałą: As ! As; " dodania do wiersza pomnożonego przez stałą innego wiersza: As + At ! As; Twierdzenie Dla A " Rnn następujące warunki są równoważne: " macierz odwrotna A-1 istnieje, tzn. macierz A jest nieosobliwa; " det(A) =0;
" wiersze macierzy A tworzą bazę Rn; " kolumny macierzy A tworzą bazę Rn; " jeżeli Ax = 0, to x = 0; def " odwzorowanie f(x) = Ax jest przekształceniem różnowartościowym f : Rn - Rn; def " odwzorowanie f(x) = Ax jest przekształceniem na f : Rn - Rn; n n " "b"R "!x"R Ax = b; " A jest iloczynem macierzy elementarnych; " 0 nie jest wartością własną macierzy A. Aatwe do rozwiązania układy równań. Systemy równań dolnie trójkąt- nych : Ą# ń# Ą# ń# Ą# ń# a11 0 0 . . . 0 x1 b1 ó# Ą# ó# a21 a22 0 . . . 0 x2 Ą# ó# b2 Ą# ó# Ą# ó# Ą# ó# Ą# ó# Ą# ó# a31 a32 a33 . . . 0 x3 Ą# = ó# b3 Ą# . (4.9) ó# Ą# ó# Ą# ó# Ą# ó# Ą# ó# Ą# ó# Ą# . . . . . . . . . . . . . . Ł# Ś# Ł# Ś# Ł# Ś# . . . . . . . an1 an2 an3 . . . ann xn bn Rozwiązanie (5.11) jest określone rekurenyjnie wzorami (oczywiście zakładamy, iż det(A) =a11a22 . . . ann =0):
b1 x1 = , a11 b2 - a21x1 x2 = , a22 . . . 16 4.L5-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 i-1
bi - aijxj j=1 xi = , aii . . . n-1
bn - anjxj j=1 xn = . ann Co daje następujący algorytm3 rozwiązywania (5.11) trójkątnych dolnie ukła- dów równań liniowych: for(i=0;i{ for(j=0;jb[i] = b[i]-A[i][j]*b[j]; b[i] = b[i]/A[i][i]; } Liczba mnożeń jaka jest wymagana przez algorytm. W instrukcji b[i] = b[i]-A[i][j]*b[j]; jest jedno mnożenie. Biorąc pod uwagę cały algorytm otrzymujemy, że całkowita liczba mnożeń wynosi: n i n
n(n +1) n2 <" 1 = i = = 2 2 i=1 j=1 i=1 Podobna też jest liczba dodawań. Ostatecznie ilość wykonywanych operacji aryt- metycznych przy rozwiązywaniu układu (5.11) jest rzędu n2. Zadania na Laboratoria z przybliżeń komputerowych 1. Znalezć największą liczbę n " N, dla której 1.0+2-n =1.0. Sprawdzić,
czy znalezione n zależy od 1.0 występującego w warunku. 2. Wykonując kolejno dzielenie przez 2 (np. 1.0) i drukując kolejne wyniki jaką najmniejszą liczbę można otrzymać. Uzasadnić otrzymany wynik teo- retycznie. 3. Niechx:=1.0/3;. Dla różnych deklaracjixwydrukować przedstawienie bitowe. Wyjaśnić otrzymany wynik. 4. Dla liczb IEEE double precision pokazać, że: " liczby tego typu są równomiernie rozłożone w odcinku [1, 2] z odstę- pem =2-52 co będzie oznaczało, że x " [1, 2] ma przedstawienie x =1 +k gdzie k =0, 1, 2, . . . , 252; 1 " pokazać, że [ , 1] zawiera tę samą ilość liczb co [1, 2] z odstępem . 2 2 3 W języku C tablice są indeksowane począwszy od 0! 4.L5-met 17 Krzysztof Tabisz 25 marca 2003, 8:39am 5 Wykład 6 Wykład 6 Dla x, y " Rn1 w przyjętym zapisie kolumnowym wektorów Ą# ń# Ą# ń# x1 y1 ó# Ą# ó# Ą# x2 y2 ó# Ą# ó# Ą# ó# Ą# ó# Ą# x3 y = y3 x = ó# Ą# ó# Ą# , ó# Ą# ó# Ą# . . . . Ł# Ś# Ł# Ś# . . xn yn mamy x " Rn1, yT " R1n i xyT " Rnn Ą# ń# Ą# ń# x1 x1y1 x1y2 x1y3 . . . x1yn ó# Ą# ó# x2 x2y1 x2y2 x2y3 . . . x2yn Ą# ó# Ą# ó# Ą# ó# Ą# x3 y1 y2 y3 . . . yn = ó# x3y1 x3y2 x3y3 . . . x3yn Ą# xyT = ó# Ą# ó# Ą# ó# Ą# ó# Ą# . . . . . . . . . . . . Ł# Ś# Ł# Ś# . . . . . . xn xny1 xny2 xny3 . . . xnyn oraz xT y " R11 Ą# ń# y1 ó# Ą# y2 n ó# Ą#
xT y = x1 x2 x3 . . . xn ó# y3 Ą# = xiyi ó# Ą# ó# Ą# . . i=1 Ł# Ś# . yn Aby uprościć zapis wzorów macierzowych wprowadzono zapis blokowy macierzy, który wyjaśniamy na przykładzie: Ą# ń# a11 a12 a13 a14 a15 a16 a17 ó# a21 a22 a23 a24 a25 a26 a27 Ą# ó# Ą# A11 A12 A13 ó# A = a31 a32 a33 a34 a35 a36 a37 Ą# = ó# Ą# A21 A22 A23 Ł# a41 a42 a43 a44 a45 a46 a47 Ś# a51 a52 a53 a54 a55 a56 a57 gdzie
a11 a12 a13 a14 a15 A11 = , A12 = , itp. a21 a22 a23 a24 a25 Przy powyższych oznaczeniach dodawanie macierzy można zapisać w postaci:
A11 A12 B11 B12 A11 + B11 A12 + B12 + = A21 A22 B21 B22 A21 + B21 A22 + B22 i w podobny sposób iloczyn macierzy:
A11 A12 B11 B12 A11B11 + A12B21 A11B12 + A12B22 = . A21 A22 B21 B22 A21B11 + A22B21 A21B12 + A22B22 Niezależnie od wprowadzanych oznaczeń, podziału macierzy A na podmacierze Aij, Bij operacje macierzowe wykonanywane są zgodne z poprzednio wprowa- dzonymi. 18 5.L6-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Na poprzednim wykładzie zajmowaliśmy się rozwiązywaniem układów dolnie trójkątnych. Podobne wzory, można wyprowadzić dla układów górnie trójkąt- nych. Obecnie zajmiemy się przypadkiem ogólnym: Ax = y. (5.10) Definicja Macierz A " Rnn jest dodatnio określona, gdy " A jest symetryczna, tzn. A = AT , " x = 0 =! xT Ax > 0, dla x " Rn.
Definicja Dla macierzy A " Rnn Ą# ń# a11 a12 a13 . . . a1n ó# a21 a22 a23 . . . a2n Ą# ó# Ą# ó# a31 a32 a33 . . . a3n Ą# A = ó# Ą# ó# Ą# . . . . . . . . . . Ł# Ś# . . . . . an1 an2 an3 . . . ann macierz Ak " Rkk, k =1, . . . , n postaci Ą# ń# a11 a12 a13 . . . a1k ó# a21 a22 a23 . . . a2k Ą# ó# Ą# ó# a31 a32 a33 . . . a3k Ą# Ak = ó# Ą# ó# Ą# . . . . . . . . . . Ł# Ś# . . . . . ak1 ak2 ak3 . . . akk nazywamy k tym minorem podstawowym macierzy A. Załóżmy, że dla macierzy A istnieją macierze trójkątne L, U Ą# ń# Ą# ń# l11 0 0 . . . 0 u11 u12 u13 . . . u1n ó# Ą# ó# l21 l22 0 . . . 0 0 u22 u23 . . . u2n Ą# ó# Ą# ó# Ą# ó# Ą# ó# l31 l32 l33 . . . 0 0 0 u33 . . . u3n Ą# L = ó# Ą# , U = ó# Ą# ó# Ą# ó# Ą# . . . . . . . . . . . . . . . . . . . . Ł# Ś# Ł# Ś# . . . . . . . . . . ln1 ln2 ln3 . . . lnn 0 0 0 . . . unn (5.11) takie, że A = LU. (5.12) Mając rozkład (5.12) rozwiązanie układu (5.10) sprowadza się do rozwiązania dwóch układów trójkątnych: Lz = b, Ux = z. Przy założeniu, iż wszystkie operacje są dobrze określone z (5.12), wykorzystu- jąc postać (5.11) (znaczna liczba elementów zerowych!) wyprowadzimy wzory określające L i U. 5.L6-met 19 Krzysztof Tabisz 25 marca 2003, 8:39am Metoda rozkładu macierzy A na czynniki trójkątne (LU rozkład): Z (5.12) mamy min(i,j) n
akj = lksusj + lkkukj, k +1d" j d" n (5.15) s=1 k-1
aik = lisusk + likukk, k +1d" i d" n (5.16) s=1 Równości (5.14), (5.15) i (5.16) są podstawą sukcesywnego wyliczenia elemen- tów macierzy L i U. Idea algorytmu polega na wyliczeniu w Kroku ą kolejnego względem k =1, . . . , n elementu leżącego na przekątnej, który w Kroku służy do wyliczenia pozostałych elementów odpowiedniego wiersza lub kolumny: Krok ą : k =1 a11 = l11u11 Krok : k =1 a1j = l11u1j, j =2, . . . , n ai1 = li1u11, i =2, . . . , n Krok ą : k =2 a22 = l21u12 +l22u22
znane Krok : k =2 a2j = l21u1j +l22u2j, j =3, . . . , n
znane ai2 = li1u12 +li2u22, i =3, . . . , n
znane Krok ą : k =3 2
a33 = l3sus3 +l33u33 s=1
znane Krok : k =3 2
a3j = l3susj +l33u3j, j =4, . . . , n s=1
znane 2
ai3 = lisus3 +li3u33, i =4, . . . , n s=1
znane . . . 20 5.L6-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Określenie lii i uii w Kroku ą nie jest jednoznaczne, dlatego też przyjęto nastę- pujące założenia dodatkowe: " lii =1, i =1, . . . , n algorytm Doolittl a; " uii =1, i =1, . . . , n algorytm Crout sa; 2 " U = LT =! lii = aii, i =1, . . . , n algorytm Cholesky ego. W ten sposób otrzymujemy wzory: k-1
lkkukk = lksusk s=1 k-1
akj - lksusj s=1 ukj = , j = k +1, . . . , n lkk k-1
aik - lisusk s=1 lik = , i = k +1, . . . , n ukk Twierdzenie Jeżeli wszystkie podstawowe minory macierzy A " Rnn są nie- osobliwe, to wtedy dla A istnieje LU rozkład. Dowód: Dowód zostanie przeprowadzony względem wymiaru ma- indukcyjnie cierzy n. Dla n =1 mamy a11 =0 jedyny minor macierzy A =[a11] " R11 i kładąc
l11 =1, u11 = a11 otrzymujemy: LU =[1] [a11] =[a11] =A. Załóżmy teraz, że dla n = k - 1 i A " R(k-1)(k-1) istnieją macierze trójkątne L, U spełniające: A = LU. Niech teraz A " Rkk. Z założenia minor podstawowy Ak-1 " R(k-1)(k-1) jest macierzą nieosobliwą. Z założenia indukcyjnego, istnieją macierze trójkątne Lk-1, Uk-1 " R(k-1)(k-1) o własności: Ak-1 = Lk-1Uk-1, det(Ak-1) =0, det(Uk-1) =0.
(5.19) det Uk-1 =0 i lk1, lk2, lk3, . . . , lkk-1 niewiadome. Oba te układy okre-
ślają jednoznacznie rozwiązanie (są Cramer owskie). Pozostaje więc wyznaczyć lkk, ukk (5.17) z równania: k-1
lksusk + lkkukk = akk s=1 gdzie można przyjąc, że lkk =1.
Twierdzenie (rozkład Cholesky ego) Jeżeli macierz A " Rnn jest dodatnio określona, to wtedy istnieje jednoznacznie określona macierz L dolnie trójkątna z lii > 0, i =1, . . . , n o własności: A = LLT . Dowód Z xT Ax > 0 dla x = (x1, x2, . . . , xk, 0, . . . , 0)T = 0 (macierz A jest
dodatnio określona) wynika, że każdy podstawowy minor Ak, k =1, 2, . . . , n jest macierzą nieosobliwą. Z poprzedniego twierdzenia wynika, że istnieją macierze trójkątne L, U o własności: A = LU. Z symetri macierzy A wynika, iż: LU = A = AT =(LU)T = UT LT . Skąd otrzymujemy: U(LT )-1 = L-1UT . 22 5.L6-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 W ostatniej równości macierz U(LT )-1 jest górnie trójkątna a L-1UT dolnie trójkątna. Istnieje więc macierz diagonalna D = U(LT )-1. Stąd otrzymujemy U = DLT i A = LDLT . Oznacza to, iż D jest dodatnio określona. Oznaczmy D1/2 macierz diagonalną o własności D = D1/2D1/2, aby otrzymać:
A = LLT , gdzie L = LD1/2.
Zadania na Ćwiczenia 1. Udowodnić, że: (a) Jeżeli L nieosobliwa jest dolnie trójkątna, to L-1 też jest dolnie trój- kątna. (b) Jeżeli L jest dolnie trójkątna i lii =1, dla i =1, . . . , n, toL-1 istnieje i na jej przekątnej są 1. (c) Iloczyn dwóch macierzy dolnie trójkątnych jest macierzą dolnie trój- kątną. 2. Udowodnić, że jeżeli A = LU jest nieosobliwa oraz lii =1, i =1, . . . n w L, to L i U są wyznaczone jednoznacznie. 3. Wykazać, że macierz górnie/dolnie trójkątna jest nieosobliwa, wtedy i tylko wtedy, gdy elementy leżące na głównej przekątnej są różne od zera. 4. Wyakazać, że dla macierzy
0 1 A = 1 1 nie istnieją macierze trójkątne L, U o własności: A = LU. 5. Z równości UU-1 = I wyprowadzić algorytm na wyznaczenie U-1, gdzie U jest górnie trójkątna. 6. Pokazać, że macierz
0 a A = 0 b ma LU rozkład. Zadania na Laboratoria 1. Opracować algorytm znajdywania L-1, gdzie L macierz dolnie trójkątna. 2. Korzystając z algorytmu Cholesky ego rozwiązać układ równań: ż# 0.05x1 + 0.07x2 + 0.06x3 + 0.05x4 = 0.23 # # # 0.07x1 + 0.10x2 + 0.08x3 + 0.07x4 = 0.32 0.06x1 + 0.08x2 + 0.10x3 + 0.09x4 = 0.33 # # # 0.05x1 + 0.07x2 + 0.09x3 + 0.10x4 = 0.31 3. Napisaćproceduręznajdującą LU rozkład macierzy A. 5.L6-met 23 Krzysztof Tabisz 25 marca 2003, 8:39am 6 Wykład 7 Wykład 7 4 Twierdzenie rozkład Choleskiego5 Jeżeli A jest dodatnio określoną macierzą, to istnieje L macierz dolnie trójkątna o własności: A = LLT . Dowód Niech Ą# ń# Ą# ń# a11 a12 . . . a1n l11 0 . . . 0 ó# ó# Ą# a21 a22 . . . a2n Ą# l21 l22 . . . 0 ó# Ą# ó# Ą# A = , L = . ó# Ą# ó# Ą# . . . . . . . . . . . . . . . . Ł# Ś# Ł# Ś# . . . . . . . . an1 an2 . . . ann ln1 ln2 . . . lnn Podzielimy macierze na bloki (patrz: wykład 6) Ą# ń# Ą# ń# a21 a22 . . . a2n ó# Ą# ó# Ą# . . . . . . . . ą = a11, a = , A" = . Ł# Ś# Ł# Ś# . . . . an1 an2 . . . ann oraz odpowiednio Ą# ń# Ą# ń# l21 l22 . . . 0 ó# Ą# ó# Ą# . . . . . . . . = l11, r = , L" = . Ł# Ś# Ł# Ś# . . . . ln1 ln2 . . . lnn Zgodnie z tymi oznaczeniami mamy, iż
ozn ą aT ozn 0 rT A = , L = oraz LT = . a A" r L" 0 LT " Z dodatniej określoności macierzy A wynika, iż ą>0 i że macierz A" (wymiar tej macierzy jest 1 mniejszy od wymiaru macierzy A, czyli znowu dowodzimy indukcyjnie !) jest też dodatnio określona. Przy takich oznaczeniach, rozkład
A = LLT zapisze się w postaci
ą aT 0 rT 2 rT = = . a A" r L" 0 LT r rrT + L"LT " " Z równości macierzy wynikają równości odpowiednich bloków ą = 2, aT = rT , A" = L"LT + rrT , " 4 Dowód z książki G. W. Stewarta 5 Inaczej liczył Tadeusz Banachiewicz (1882 1954), prof. UJ, który wprowadził rachunek krakowianowy 24 6.L7-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 i dalej przekształcając, otrzymujemy " = ą, (6.20) rT = -1aT , (6.21) L"LT = A" - rrT . (6.22) " Pierwsze dwa równania (6.20), (6.21) są podstawą do wyliczenia pierwszej ko- lumny macierzy L. Ponieważ ą > 0, to = 0 i element l11 = jest dobrze
określony. Z =0 i (6.21) wynika, iż macierz rT =[l21, . . . , ln1] jest jednoznacz-
nie określona. Pozostaje udowodnić, iż macierz z prawej strony równości (6.22) jest dodatnio określona. Niech 1 Ć A" ozn A" - rrT = A" - aaT . = ą Ć Symetryczność macierzy AT wynika z równości (oczywiście A" = AT ) " " Ć Ć AT =(A" - rrT )T = AT - (rT )T rT = A" - rrT = A". " " Niech y " Rn1, wtedy yT a = aT y i dla y =0 udowodnimy
1 1 Ć yT A"y = yT A"y - (yT a)(aT y) =yT A"y - (aT y)2 > 0. (6.23) ą ą Macierz A jest dodatnio określona, więc z formy jej zapisu dla dowolnego " R wynika
ą aT 0 < [ yT ] = ą2 +2aT y + yT A"y. (6.24) a A" y 1 Do nierówności (6.24) podstawiamy = - aT y ą 1 0 <ą2 +2aT y + yT A"y = yT A"y - (aT y)2 ą a to jest dokładnie warunek (6.23). W ten sposób udowodniliśmy krok induk- cyjny względem wymiaru macierzy.
W przypadku algorytmu Cholesky ego wzory na LU rozkład macierzy A (patrz wykład 6) redukują się do postaci:
k-1
akk - lks 2 lkk = s=1 k-1
aik - lislks s=1 lik = , i = k +1, . . . , n lkk 6.L7-met 25 Krzysztof Tabisz 25 marca 2003, 8:39am Metoda eleminacji Gaussa Tak jak w przypadku LU rozkładu w pierwszej kolejności, przy założeniu do- puszczalności wszystkich operacji, będziemy upraszczać układ równań ż# a11 x1 + a12 x2 + . . . + a1n xn = b1 # # # a21 x1 + a22 x2 + . . . + a2n xn = b2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . # # # an1 x1 + an2 x2 + . . . + ann xn = bn Postępując naiwnie otrzymujemy odpowiednio dzieląc i mnożąc, równoważny układ równań ż# a11x1 + a12x2 + . . . + a1nxn = b1 # # #
# # # a21a12 a21a1n a21b1 # 0 + a22 - x2 + . . . + a2n - xn = b2 - a11 a11 a11 # . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . # # # # # # an1a12 an1a1n an1b1 0 + an2 - x2 + . . . + ann - xn = bn - a11 a11 a11 Zapowiada się, iz jezeli będziemy odpowiednio wytrwali (wszystko przy zało- żeniu braku problemów z dzieleniem) to otrzymamy równoważny układ górnie trójkątny. W przypadku ogólnym, aby z macierzy A(k): Ą# ń# a(k) . . . a(k) a(k) . . . a(k) . . . a(k) 1 1,k-1 1k 1j 1n ó# Ą# ó# Ą# . . . . ó# Ą# . . . . . . . . ó# Ą# ó# Ą# ó# Ą# ó# Ą# a(k) a(k) . . . a(k) . . . a(k) k-1,k-1 k-1,k k-1,j k-1,n ó# Ą# ó# Ą# ó# 0 . . . 0 a(k) . . . a(k) . . . a(k) Ą# ó# Ą# kk kj kn ó# Ą# ó# Ą# ó# Ą# 0 . . . 0 a(k) . . . a(k) . . . a(k) (6.25) A(k) = k+1,k k+1,j k+1,n ó# Ą# ó# Ą# ó# . . . . . Ą# . . . . . ó# Ą# . . . . . ó# Ą# ó# Ą# ó# 0 . . . 0 a(k) . . . a(k) . . . a(k) Ą# ó# Ą# ik ij in ó# Ą# ó# Ą# . . . . . ó# Ą# . . . . . ó# . . . . . Ą# Ł# Ś# 0 . . . 0 a(k) . . . a(k) . . . a(k) nn nk nj otrzymać macierz A(k+1) należy zastosować wzory: ż# # a(k), dla i d" k, # ij # # # a(k) a(k+1) = ik (6.26) ij a(k) - a(k), dla i e" k +1, j e" k +1, ij # a(k) kj # kk # # # 0, dla i e" k +1, j d" k, oraz dla kolumny wyrazów wolnych ż# # b(k), dla i d" k, # i b(k+1) = . (6.27) i # a(k) ik # b(k) - b(k), dla i e" k +1 i a(k) k kk 26 6.L7-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Przedyskutujemy teraz dopuszczalność przedstawionej procedury, która jest warun- kowana przez a(k) =0. Z postaci macierzy (6.25) wynika, iż
kk
a(k) . . . a(k)
kk kn
. .
. . det(A(k)) =a(k) . . . a(k) . 11 k-1,k-1 . .
a(k) . . . a(k) nn nk Z własności wyznaczników oraz postaci wzorów (6.26) wynika, że det(A) =det(A(1)) =det(A(2)) = =det(A(n)). W ten sposób otrzymaliśmy, że jeżeli det(A) =0 to dla każdego kroku k istnieje
i " {k, . . . , n} takie, że a(k) = 0. Więc w przypadku a(k) = 0 przy założeniu
ik kk det(A) =0, należy zamienić wiersz k ty z wierszem i tym. Zamiana taka pro-
wadzi oczywiście do równoważnego układu równań. Udowodniliśmy w ten sposób: Twierdzenie Gauss Jeżeli det(A) = 0, to wzory (6.26) i (6.27) z uwzględ-
nieniem zamiany wierszy w przypadku, gdy w k-tym kroku a(k) =0 sprowadza kk układ Ax = b do równoważnego Ux = b(n) gdzie macierz U jest górnie trójkątna. Przykład:6 W przykładzie bedziemy tylko stosować wzory (6.26) i (6.27), bez zamiany wierszy w przypadku, gdy w k-tym kroku a(k) =0 . kk Więc przy takich założeniach nie poradzimy sobie z układem
0 1 x1 1 = 1 1 x2 2 Poprawimy ten układ. Więc dla 0 < do układu:
1 x1 1 = 1 1 x2 2 możemy już procedurę zastosować. Otrzymujemy wtedy
1 x1 1 = . 0 1 - -1 x2 2 - -1 Rozwiązując otrzymany układ równań ( górnie trójkątny ) otrzymujemy, gdzie znakiem H" zaznaczyliśmy wpływ przybliżeń komputerowych na obliczenia:
2--1 x2 = H" 1 1--1 1-x2 x1 = H" 0 -1 6 za D. Kincaid i W. Cheney 6.L7-met 27 Krzysztof Tabisz 25 marca 2003, 8:39am Sprowadzenie układu do postaci
1 -1 x1 -1 = 1 1 x2 2 nie poprawia sytuacji
1 -1 x1 -1 = 0 1 - -1 x2 2 - -1 otrzymujemy tak jak poprzednio
1 1 x1 2 = 0 1 - x2 1 - 2 który, nie jest już tak czuły na przybliżenia komputerowe.
1-2 x2 = H" 1 1- x1 = 2 - x2 H" 1 Aby wyrównać szeregi przed zapowiedzianym kolokwium z rozwiązywania rów- nań nieliniowych na tej liście nie ma Zadań na Ćwiczenia i Laboratoria. 28 6.L7-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 7 Wykład 8 i 9 Rozwiązania zadań z ćwiczeń A. Metoda Steffensen a z listy do wykładu z dnia ??.XI.2000r Niech f(xn) f(x + f(x)) - f(x) xn+1 = xn - , gdzie g(x) = . (7.28) g(xn) f(x) Przy odpowiednich założeniach, wykazać kwadratową zbieżność metody. Rozwiązanie Ze wzoru Lagrange a mamy: f(x + f(x)) - f(x) = f (x + f(x)), dla pewnego 0 < <1. f(x) Wprowadzamy oznaczenie: xn ozn. xn + nf(xn). Wtedy wzór (7.28) można Ż = zapisać w następujący sposób: f(xn) xn+1 = xn - . f (Ż xn) W tej postaci widać podobieństwo (7.28) ze wzorem iteracyjnym z metody New- ton a. Dalej będziemy korzystali z oznaczeń i rozumowań przeprowadzonych przy dowodzie zbieżności metody Newton a. Dla pierwiastka równania f(r) =0 oznaczamy błąd przybliżenia en = xn - r. Mamy wtedy n ozn. xn - r = (7.29) = Ż = xn - r + n(f(xn) - f(r)) = en(1 + nf (śn)), śn " (r, xn). Tak jak w dowodzie metody Newton a wyliczamy: f(xn) enf (Ż - f(xn) xn) en+1 = xn+1 - r = xn - - r = . (7.30) f (Ż f (Ż xn) xn) Z dowodu metody Newton a mamy: 1 enf (xn) - f(xn) = f (n)e2 , gdzie n " (r, xn). n 2 Ze wzoru Lagrange a zastosowanego do pochodnej funkcji y = f (x) otrzymu- jemy (z (7.29) i f(r) =0 mamy f(xn) =f(xn) - f(r) =f (śn)en): f (Ż =f (xn + nf(xn)) = xn) = f (xn) +f (n)nf(xn) =f (xn) +f (n)nf (śn)en, czyli 1 enf (Ż - f(xn) = f (n)e2 + f (n)nf (śn)e2 . xn) n n 2 7.L8-met 29 Krzysztof Tabisz 25 marca 2003, 8:39am Wstawiając ostatnie wyrażenie do (7.30) otrzymujemy: 1 f (n) +f (n)nf (śn) 2 en+1 = e2 . (7.31) n f (Ż ) x n ograniczone w otoczeniu r Występujące punkty pośrednie są wybrane w następujący sposób: n, śn " (r, xn), xn = xn + nf(xn), n " (xn, xn + nf(xn)), 0 <n < 1 Ż gdzie występujące odcinki mogą mieć zmienioną orientację. Określamy teraz tak jak w dowodzie metody Newton a funkcję
1 max |f (x)| + max |f (x)| 2 |x-r|d" |x-r|d" def c() = . min |f (x)| |x-r|d" Teraz, jeżeli y = f (x) jest funkcją ciągłą i f (r) =0, to
lim c() =0. 0 ozn. Istnieje więc taka >0, że = c() < 1. Proces iteracji zaczniemy z punktu x0 spełniającego:
|x0 - r| < oraz |x0 - r| 1+| max |f (x)| < 1. |x-r|d" Z ostatniej nierówności oraz z (7.29) i oszacowania |xn - r| d"n|x0 - r| wynika, że jeżeli lim xn = r to lim xn = r co pozwala kontrolować zachowanie miano- Ż n " n " wnika w (7.31). W ten sposób znalezliśmy się w warunkach dowodu zbieżności ciągu przybliżeń en z metody Newton a.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. Przybliżenie y = f (x) z listy do wykładu z dnia ??.XI.2000r Korzystając ze wzoru Taylor a uzasadnić k2f(x + h) - h2f(x + k) +(h2 - k2)f(x) f (x) H" (k - h)kh Rozwiązanie Będziemy korzystali z wzoru Taylor a z resztą w postaci całki7 x n
rn(x) 7 G.M. Fichtenholz, Rachunek różniczkowy i całkowy, T 2, PWN, Warszawa 1976, ss. 126. 30 7.L8-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Skorzystamy też z twierdzenia o wartości średniej dla rachunku całkowego: Jeżeli f, g : [a, b] - R są ciągłe i funkcja y = g(x) ma na [a, b] stały znak, to istnieje " [a, b] o własności b b f(x)g(x) dx = f() g(x) dx. a a Korzystamy ze wzoru Taylor a dla n =1 aby otrzymać: k2f(x + h) - h2f(x + k) +(h2 - k2)f(x) = (k - h)hk k2 [f(x) +f (x)h + r1(x + h)] - h2 [f(x) +f (x)k + r1(x + k)] = + (k - h)hk (h2 - k2)f(x) + = (k - h)hk f (x)(hk2 - h2k) +k2r1(x + h) - h2r1(x + k) = = (k - h)hk k2r1(x + h) - h2r1(x + k) = f (x) + . (k - h)hk
(h,k) Pokażemy, że dla = (h, k) k2r1(x + h) - h2r1(x + k) ozn. (h, k) = (7.32) (k - h)hk zachodzi lim (h, k) =0. (h,k)(0,0) Przyjmijmy, że 0<|k|<|h| (założenie to nie ogranicza ogólności rozważań). Wte- dy k2r1(x + h) - h2r1(x + k) = (7.33) x+h x+k
k2 h2 = f (t)(x + h - t) dt - f (t)(x + h - t) dt = 2 2 x x x+h x+k
k2 1 = f (t)(x + h - t) dt + f (t) k2(x + h - t) - h2(x + k - t) dt. 2 2 x+k x
I2 I1 Każdą z zaznaczonych całek I1, I2 będziemy szacowali oddzielnie. Pierwszą z nich szacujemy następująco x+h x+h
k2 f (1)k2 I1 = f (t)(x + h - t) dt = (x + h - t) dt = 2 2 x+k x+k
f (1)k2 (x + h - t)2 x+h f (1)
= - = k2(h - k)2 2 2 x+k 4 7.L8-met 31 Krzysztof Tabisz 25 marca 2003, 8:39am gdzie 1 " [x + k, x + h]. Drugą szacujemy podobnie x+k
1 I2 = f (t) k2(x + h - t) - h2(x + k - t) dt = 2 x x+k x+k
(k2 - h2) (k - h)hk = f(t)(x - t) dt + f(t) dt = 2 2 x x
(k2 - h2)f (2) (x - t)2 x+k (k - h)hk2
= - + f (3) = 2 2 2 x f (2) f (3) = -(k2 - h2)k2 +(k - h)hk2 gdzie 2, 3 " [x, x + k]. 4 2 Jeżeli założymy, że druga pochodna y = f (x) jest ograniczona, to otrzymamy następujące oszacowania dla całek I1, I2 |I1| d"k2(h - k)2C i |I2| d"(|(k2 - h2)k2| + |(k - h)hk2|)C. W ten sposób wykazaliśmy (z (7.32), (7.33) i założenia 0<|k|<|h|), iż: |I1| |I2| |(h, k)| d" + d" |(k - h)hk| |(k - h)hk|
" d" (|h - k| + |k + h| + |k|)C d" 3(|h| + |k|)C d" 3 2 h2 + k2C. Z ostatniej nierówności wynika lim (h, k) =0. (h,k)(0,0)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Wykład 8 Wrócimy jeszcze do rozwiązywania układu równań metodą eleminacji Gaussa przedstawionej na wykładzie w dniu 29.XI.2000r oraz problemu obliczeniowego wynikającego z przykładu. W trakcie eleminacji niewiadomych w metodzie Gaussa powstają macierze postaci Ą# ń# X . . . X X . . . X ó# . . . Ą# . . . . . ó# Ą# . . . . ó# Ą# ó# Ą# X X . . . X ó# Ą# ó# 0 . . . 0 akk . . . akn Ą# ó# Ą# ó# Ą# . . . . . . . . Ł# Ś# . . . . 0 . . . 0 ank . . . ann gdzie8 X em są zaznaczone elementy macierzy niekoniecznie równe zeru. W ko- lejnym k tym kroku wymagane jest aby akk = 0, gdyby ten warunek nie był
8 za J. Wilkinson em 32 7.L8-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 spełniony to wystarczy wśród akk, . . . , ank znalezć element niezerowy (istnienie wynika z nieosobliwości macierzy) i zamienić odpowiednie wiersze. Z uwagi na dokładność obliczeń wybór l-tego wiersza spełniającego |alk| = max |aik| (czę- kd"id"n ściowy wybór elementu podstawowego, partial pivoting ) poprawia dokładność obliczeń. Jeszcze większą dokładność można otrzymać uwzględniając w wybo- rze pozostałe kolumny macierzy |alm| = max max |aij| (pełny wybór elementu kd"id"n kd"jd"n podstawowego, complete pivoting ). Tak wybrany element przez zamianę wier- szy i kolumn umieszczamy na miejscu (k, k) itd. W przypadky rozwiązywania dużych układów równań wybór ten może być kłopotliwy. Rozważa się oddzielnie macierze z dużą ilością zer (wzory rozwiązujące pozo- stają te same, można jedynie oszczędniej zapisać macierz w komputerze). Takim przykładem jest macierz Hessenberga gdzie dla i >j +1 =! aij =0 i pasmowa (trójdiagonalna) gdzie z |i - j| > 1 =! aij =0: Ą# ń# Ą# ń# X X X . . . X X X d1 c1 ó# Ą# ó# Ą# X X X . . . X X X a1 d2 c2 ó# Ą# ó# Ą# ó# Ą# ó# Ą# 0 X X . . . X X X a2 d3 c3 ó# Ą# ó# Ą# , ó# Ą# ó# Ą# . 0 0 X . . . X X X . . . . . . ó# Ą# ó# Ą# . . . ó# Ą# ó# Ą# . . . . . . . . . . . . . . Ł# Ś# Ł# . . . . . . . a-2 dn-1 cn-1 Ś# 0 0 0 . . . 0 X X an-1 dn Analiza błędów9 10 11 12 13 7.1 Definicja Dla macierzy A " Rnn Ą# ń# a11 a12 . . . a1n ó# a21 a22 . . . a2n Ą# ó# Ą# A = ó# Ą# . . . . . . . . Ł# Ś# . . . . an1 an2 . . . ann funkcję : Rnn - R+ nazywamy normą macierzową, gdy ma następujące własności 1. A e"0, oraz A =0 !! A =0, 2. cA = |c| A , c " R, 3. A + B d" A + B , 4. AB d" A B . Typowe normy
n n
A 2 = a2 (norma euklidesowa, 2 norma) ij i=1 j =1 A S = max i(AAT ) (norma spektralna) i 9 D. Kincaid & W. Cheney, Numerical Analysis, ss. 181 221; 10 G.W. Stewart, Afternotes on Numerical Analysis, ss. 103 131; 11 A. Ralston, Wstęp do analizy numerycznej, ss. 365 392; 12 R.L. Burden, J.D. Faires & A.C. Reynolds, Numerical Analysis, ss. 289 317; 13 H. Rutishauser, Lectures on Numerical Mathematics, ss. 23 76; 7.L8-met 33 Krzysztof Tabisz 25 marca 2003, 8:39am Gdzie w definicji normy14 spektralnej A S symbolem i(AAT ) oznaczono i tą wartość własną macierzy AAT . Uwaga Z macierzą A " Rnn jest związane wzorem f(x) =Ax przekształcenie liniowe f : Rn - Rn. Przy tak przyjętych oznaczeniach mamy f(x) def A = f = sup f(x) =sup . x =0 x =1 x Inaczej, norma macierzy A jest normą operatorową przekształcenia y = f(x) wyznaczonego przez tą macierz. Ogólnie, dla przekształcenia f : B1 - B2, gdzie B1, B2 przestrzenie unor- mowane, rozważa się klasę przekształceń ograniczonych. Dla takich przekształceń f(x) def f = sup x x =0 jest normą. Normy w przestrzeniach Rn: n
x 1 = |xi| 1 norma, i=1
n
x 2 = x2 2 norma, norma euklidesowa, i i=1 x " = max |xi| " norma, 1d"id"n gdzie xT =(x1, x2, . . . , xn) " Rn. 7.2 Lemat Jeżeli przestrzeń Rn jest z normą x = max |xi|, x " Rn to dla A " Rnn 1d"id" zachodzi n
A " = max |aij|. 1d"id"n j=1 Uzasadnienie
A " = sup Au " = sup max |(Au)i| = 1d"id"n u "=1 u "=1
= max sup |(Au)i| = max sup |(Au)i| = 1d"id"n 1d"id"n u "=1 u "=1 n
= max |aij| 1d"id"n j=1
n n
gdzie kres górny sup aijuj = |aij|, i = 1, . . . , n jest realizowany
j=1 u "=1 j=1 przez u " =1 określone w następujący sposób uj = sgn aij, i =1, . . . , n. 14 I.M. Gelfand, Wykłady z algebry liniowej, PWN, Warszawa 1971. 34 7.L8-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003
Przykład I Niech x " Rn będzie rozwiązaniem dokładnym równania Ax = b a x " Rn równania zaburzonego
Bx = b.
Pytanie: Jaki jest błąd względny otrzymanego wyniku ? Rozwiązanie: Mamy x = A-1b rozwiązanie dokładne, x = B-1b rozwiązanie przybliżone.
Stąd wynika x - x = x - Bx = x - BA-1x = (I - BA-1)x d" I - BA-1 x
Co daje następujące oszacowanie błędu względnego x - x
d" I - BA-1 . x
Przykład II Macierze układów są dokładne, natomiast prawe strony są zaburzone, tzn. Ax = b równanie dokładne ,
A = b równanie z zaburzoną prawą stroną. x Pytanie: (tak jak poprzednio) Jaki jest błąd względny otrzymanego wyniku ? Mamy x = A-1b rozwiązanie dokładne, x = A-1 rozwiązanie przybliżone. b Teraz jest
x - x d" A-1b - A-1 = A-1(b - b) d" A-1 b - b . b Skąd wynika (wiemy, że Ax = b)
Bx b - b
b - b d" A-1 b - b d" A-1 A x b b Otrzymaliśmy w ten sposób następujące oszacowanie błędu względnego
x - x b - b
d" A-1 A x b
(A) 7.L8-met 35 Krzysztof Tabisz 25 marca 2003, 8:39am Liczba (A) ( condition number ) związana z macierzą nieosobliwą A " Rnn def (A) = A-1 A odgrywa zasadniczą rolę w oszacowaniu błędu względnego rozwiązań układu rów- nań przez błąd względny prawych stron
x - x b - b
d" (A) . x b
Przykład III (liczbowy) Dla >0 niech
1 1 + A = . 1 - 1 Mamy det(A) =2, więc macierz A jest nieosobliwa. Wtedy macierz odwrotna A-1 jest określona wzorem
1 1 -1 - A-1 = . 2 -1+ 1 Dla normy macierzowej " mamy 2+ A " =2 + i A-1 " = 2 skąd otrzymujemy, iż 2 2+ 4 (A) = > . 2 2 Jeżeli przyjmiemy, że =0.01 to otrzymamy (A) > 40 000. Daje to pogląd o niebezpieczeństwach związanych z rozwiązywaniem układów równań.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 7.L8-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Wykład 9 Jeżeli rozwiązujemy układ równań Ax = b numerycznie nie otrzymujemy dokładnego rozwiązania dokładnego x a jedynie przybliżone x. Oszacowanie błędów względnych jakie powstają jest dane przez
7.3 Twierdzenie Jeżeli x jest rozwiązaniem dokładnym układu Ax = b a x przybliżonym
Ać H" b, to wtedy błąd względny przybliżonego rozwiązania x spełnia
1 b - A x - x b - A x x d" d" (A) , (7.34) (A) b x b gdzie def (A) = A-1 A . Dowód Dowodzimy prawej nierówności w (7.34). Mamy x - x b = A-1(b - A Ax d" A-1 b - A A x x) x skąd otrzymujemy
x - x b - A x d" A-1 A . x b
(A) Teraz lewa strona (7.34). Zachodzi b - A x = A(x - x) A-1b d" A x - x A-1 b x co daje
b - A x - x x d" A-1 A . b x
(A)
W oszacowaniu błedu względnego numerycznego rozwiązania x poza dokładno-
ścią wykonywanych obliczeń, istotną rolę pełni macierz A układu. Liczba (A) określa czułość układu na błędy zaokrągleń. Macierze dla których liczba (A) jest duża nazywamy zle uwarunkowane. W przypadku analizy dokładności rozwiązania normę macierzową dobiera się w sposób umożliwiający przeprowadzenie obliczeń (lub oszacowań) na (A). Przeważnie jest to norma spektralna. 7.L8-met 37 Krzysztof Tabisz 25 marca 2003, 8:39am Zadania na Ćwiczenia 1. Pokazać, że funkcje x ", x 2, x 1 są normami w Rn. 2. Wykazać nierówności x " d" x 2 d" x 1 dla x " Rn. 3. Wykazać, że jeżeli f, g : Rn - R są ograniczone, to wtedy zachodzi sup[f(x) +g(x)] d" sup f(x) +supg(x). n
4. Jeżeli w Rn jest określona norma wektorowa x = |xi|, to odpowia- i=1 dająca norma macierzowa jest określona wzorem n
A " = max |aij|. 1d"jd"n i=1 5. Dla normy macierzowej " obliczyć (A) dla macierzy
7 8 A = . 9 10 6. Niech n =3 i dla macierzy Ą# ń# 4 -3 2 Ł# Ś# A = -1 0 5 2 6 -2 obliczyć sup Ax ". x "d"1 7. Obliczyć (A) dla norm 1, 2, " i macierzy
a +1 a 0 1 ą 1 a) , b) , c) . a a - 1 -2 0 1 1 1 n p 8. Dla p e" 1 i x " Rn określamy x p def ( |xi|) . Wykazać, że " = i=1 jest normą dla której zachodzi limp " x p = x ". Zadania na Laboratoria 1. Napisać program rozwiązujący układ równań metodą eleminacji Gaussa z pełnym wyborem elementu podstawowego. 2. Opracować procedurę wyliczającą dla danej macierzy A współczynnik (A). 3. Podać przykłady dobrze i zle uwarunkowanych macierzy wraz z otrzy- manymi wyznaczonymi przez nie numerycznymi rozwiązaniami układów równań. 38 7.L8-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 8 Wykład 10 Uzupełnienie wykładu 9 Oszacowanie liczby wykonywanych operacji długich (mnożenia i dzielenia) dla poszczególnych algorytmów. Uwaga: Dawniej obliczenia były wykonywane ręcznie lub na mechanicznych arytmometrach, gdzie wykonanie mnozenia czy też dzielenia stanowiło czyn- ność dość skomplikowaną. Stąd ilość wykonywanych operacji długich było mia- rą czasu jaki był wymagany przez przyjęty algorytm. Obecnie na podstawie eksperymentu stwierdziłem, iz mnożenie może być wykonywane szybciej niz do- dawanie!? Zalezy to od procesora w jaki jest wyposażony komputer. W praktyce czas realizacji obliczeń o wiele bardziej jest zależny od umiejętności wykorzy- stania przez programistę (w przypadku komputerów typu PC XT/AT) pamięci operacyjnej i rozszerzonej oraz na ograniczeniu operacji dyskowych. W przy- padku wykonywania obliczeń w systemach wielozadaniowych (Windows NT, Novell, UNIX, Linux) czas ten zależy dodatkowo od ich aktualnego obciążenia. Liczba operacji długich wynosi15 " Lx = b układów dolnie trójkątnych n(n +1) n2 H" , 2 2 " LLT = b dla algorytmu Choleskiego (wyznaczenie macierzy L) n3 , 6 " Ax = b w metodzie eleminacji Gauss a n2 . 3 Przykład (wpływ skalowania macierzy na dokładność) W przypadku teoretycznym rozwiązywania układu Ax = b, pomnożenie wiersza przez liczbę różną od zera nie wpływa na jego rozwiązanie. Dla numerycznego wyznaczenia rozwiązań jest inaczej. Zajmiemy się współczynnikiem (A) = A " A-1 " występującym w oszacowaniu błędów wzgędnych. Niech
1 1 2 -1 A = wtedy A-1 = oraz (A) =3 3 =9. 1 2 -1 1 Mnożymy teraz pierwszy wiersz układu Ax = b przez liczbę 10-4 aby otrzymać:
10-4 10-4 2 104 -1 = wtedy A-1 = 1 2 -104 1 oraz (A) =3 (2 104 +1) H" 6 104. 15 za: G.W. Stewart, Afternotes on Numerical Analysis, SIAM 1996, ss. 94 102. 8.L9-met 39 Krzysztof Tabisz 25 marca 2003, 8:39am Zgodnie z tym co zostało udowodnione na wykładzie w dniu 13.XII.2000r współ- czynnik (A) dobrego uwarunkowania macierzy, uczestniczył w oszacowaniu błędu względnego przybliżonego rozwiązania x
1 b - A x - x b - A x x d" d" (A) . (A) b x b Z przykładu widać jak dużo zależy od wyboru normy macierzowej, jak i też od skalownia macierzy. Wykład 10 Inetrpolacja wielomianowa Zadanie: Dla danych n +1 par (xi, yi) liczb x x0 x1 x2 . . . xn y y0 y1 y2 . . . yn wyznaczyć wielomian pn(x) możliwie najniższego stopnia o własności pn(xi) =yi dla 0 d" i d" n. 8.1 Twierdzenie Jeżeli x0, x1, . . . , xn są liczbami spełniającymi warunek xi = xj dla i = j oraz
y0, y1, . . . , yn są dowolnymi liczbami, to wtedy istnieje jedyny wielomian pn co najwyżej stopnia n o własności: pn(xi) =yi dla 0 d" i d" n. Dowód Jedyność. Niech pn, i qn będą dwoma takimi wielomianami. Wtedy wilomian pn - qn jest stopnia co najwyżej n oraz spełnione jest pn(xi) - qn(xi) =0 dla 0 d" i d" n. Ponieważ wielomian stopnia nie większego od n ma n +1 pierwiastków, więc zachodzi pn(x) - qn(x) a" 0. Istnienie (dowód indukcyjny). Wielomiany pk dla k =0, . . . , n są wyznaczane ze wzoru ż# y0, dla k =0, # pk(x) = # pk-1(x) +ck(x - x0)(x - x1) . . . (x - xk-1), dla k =1, . . . , n. (8.35) gdzie yk - pk-1(xk) ck = . (8.36) (xk - x0)(xk - x1) . . . (xk - xk-1)
Wielomiany interpolacyjne pk(x), k =1, . . . , nw (8.35) nazywają sie wielomianami interpolacyjnymi Newtona. Samo obliczanie warości wielomianów może odby- wać sie według schematu Hornera pk =( (((ck)dk-1 + ck-1)dk-2 + ck-2)dk-3 + . . . + c1)d0 + c0, 40 8.L9-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 gdzie dk = x - xk. Wielomiany interpolacyjne Lagrange a Metoda wynika z następującego przedstawienia n
pn(x) =y0(x) +y1 1(x) +. . . + yn n(x) = yk k(x). k=0 Gdzie wielomiany 0(x), 1(x), . . . , n(x) zależa od danych x0, x1, . . . , xn a nie od y0, y1, . . . , yn. Dla yi = ij, i, j =1, . . . n zachodzi więc n n
ij = pn(xj) = yk k(xk) = kj k(xj) = i(xj). k=0 k=0 Jeżeli dodatkowo założymy, że i(x) są wielomianami stopnia n o własności i(xj) =0 dla i = j,
oraz i(xi) =1. Z poczynionych założeń wynika, że wielomiany i(x), i =0, . . . n są postaci
i(x) =c(x - x0) . . . (x - xi) . . . (x - xn),
gdzie (x - xi) oznacza, iż ten czynnik został opuszczony. Ostatecznie otrzymu- jemy, iż n
x - xj i(x) = . xi - xj j=0 j =i Metoda bezpośrednia Z warunku pn(xi) =yi, dla i =0, . . . , n gdzie pn(x) =anxn + an-1xn-1 + . . . + a1x + a0 otrzymujemy układ równań Ą# ń# Ą# ń# Ą# ń# 1 x0 x2 . . . xn a0 y0 0 0 ó# Ą# ó# 1 x1 x2 . . . xn a1 Ą# ó# y1 Ą# 1 1 ó# Ą# ó# Ą# ó# Ą# = ó# Ą# ó# Ą# ó# Ą# . . . . . . . . . . . . . . . Ł# Ś# Ł# Ś# Ł# Ś# . . . . . . . 1 xn x2 . . . xn an yn n n Otrzymaliśmy w ten sposób układ równań typu Va = y na współczynniki aT =[a0, a1, . . . , an] gdzie macierz V jest macierzą Vandermonde a. Ten sposób wyznaczania wielomianów iterpolacyjnych nie jest zalecany. 8.L9-met 41 Krzysztof Tabisz 25 marca 2003, 8:39am Błąd w wyznaczaniu wielomianów iterpolacyjnych 8.2 Twierdzenie Niech funkcja f będzie klasy Cn+1[a, b] i niech p będzie wielomianem stopnia nie większego od n o własności p(xi) =f(xi) gdzie xi " [a, b] dla i =0, . . . , n. Wtedy dla każdego x " [a, b] istnieje x " (a, b) o własności n
1 f(x) - p(x) = f(n+1)(x) (x - xi). (8.37) (n +1)! i=0 Dowód Ustalamy x " [a, b], x = xi i =0, . . . n. Niech
n
w(t) = (t - xi), (t) =f(t) - p(t) - w(t), dla t " [a, b], (8.38) i=0 gdzie stała jest wyznaczona z warunku (x) =0. Skąd otrzymujemy, iż f(x) - p(x) = . w(x) Mamy teraz " Cn+1[a, b] i (t) = 0 dla t " {x, x0, x1, . . . , xn}. Stąd funk- cja = (t) ma na odcinku [a, b] (n +2) różnych zer. Z twierdzenia o war- tości średniej Rolle a funkcja = (t), " Cn[a, b] ma na odcinku (a, b) (n +1) różnych miejsc zerowych. Kontynuując ten proces otrzymujemy, iż funk- cja (n+1) = (n+1)(t) ma na odcinku (a, b) jedno zero x " (a, b). Wyliczamy teraz (n+1)(t) =f(n+1)(t) - p(n+1)(t) - w(n+1)(t) =f(n+1)(t) - (n +1)! . Skąd otrzymujemy f(x) - p(x) 0 =(n+1)(x) =f(n+1)(x) - (n +1)! = f(n+1)(x) - (n +1)! . w(x) Co po uwzględnieniu określenia w = w(t) z (8.38) daje (8.37).
42 8.L9-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 9 Wykład 11 Wykład 11 Wielomiany iterpolacyjne Czebyszewa16 Wielomiany Czebyszewa pierwszego radzaju określamy rekurencyjnie w nastę- pujący sposób: ż# 1, dla n =0; # Tn(x) = x, dla n =1; # 2xTn(x) - Tn-1(x), dla 1 Postać kilku dalszych wielomianów Czebyszewa: T2 = 2x2 - 1, T3 = 4x3 - 3x, T4 = 8x4 - 8x2 +1, T5 = 16x5 - 20x3 +5x, T6 = 32x6 - 48x4 +18x2 - 1, . . . 9.1 Twierdzenie Jeżeli x " [-1, 1], to wielomiany Czebyszewa są określone wzorem Tn(x) =cos(n arccos x). (9.39) Dowód Wiadomo, iż cos(ą + ) =cos ą cos - sin ą sin . Skąd otrzymujemy cos(n +1) = cos cos n - sin sin n, cos(n - 1) = cos cos n +sin sin n, oraz cos(n +1) = 2 cos cos n - cos(n - 1). (9.40) Wstawiając do (9.40) za = arccos x (czyli x =cos ) dostajemy, iż funkcja fn(x) =cos(n arccos x) spełnia ż# f0(x) = 1, # # # f1(x) = x, # # # fn+1(x) = 2xfn(x) - fn-1(x) dla 1 d" n. 16 1821 1894 9.L10-met 43 Krzysztof Tabisz 25 marca 2003, 8:39am
Ze wzoru (9.39) wynikają następujące własności wielomianów Czebyszewa: |Tn(x)|d" 1 dla x " [-1, 1]
k Tn cos Ą = (-1)k, k =0, 1, . . . , n (9.41) n
2k-1 Tn cos Ą = 0,k =1, 2, . . . , n. 2n 9.2 Twierdzenie Dla wielomianów Czebyszewa jest spełniony następujący warunek ortogonalno- ści: ż# 0, dla n = m,
# # # 1 dx " Tn(x) Tm(x) = Ą, n = m =0, 1-x2 # -1 # # Ą , n = m =0.
2 Dowód W dowodzie ortogonalności skorzystamy z postaci (9.39) wielomianów oraz podstawienia: x() = cos , dx() = - sin . d Zachodzi więc 1 dx Tn(x) Tm(x)" = 1 - x2 -1 1 dx " = cos(n arccos x) cos(m arccos x) = 1-x2 -1 Ą 0 - sin d " = cos(n)cos(m) = cos n cos m d = 1-cos2 Ą 0 ż# 0, dla n = m,
# # Ą #
cos(n+m)+cos(n-m) = = Ą, n = m =0, 2 0 # # # Ą , n = m =0.
2
Z twierdzenia wynika, iż fukcję y = f(x) można przedstawić w postaci szeregu Czebyszewa: " c0 f(x) <" + cnTn(x) 2 n=1 gdzie współczynniki (cn) " są określone wzorem: n=0 1 dx cn = f(x) Tn(x) " . 1 - x2 -1 Podobnie jak w przypadku szeregu Fouriera należy wykazać zbieżność otrzyma- nego szeregu do zadanej funkcji y = f(x). 44 9.L10-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 W następnym twierdzeniu skorzystamy z postaci n tej funkcji Czebyszewa dla której współczynnik przy najwyższej potędze wynosi 2n-1: Tn(x) =2n-1xn + an-1xn-1 + . . . + a1x + a0 dla 1 d" n gdzie nie ma znaczenia zerowanie się współczynnika an-1. 9.3 Twierdzenie Jeżeli pn(x) =1 xn + an-1xn-1 + . . . + a1x + a0, to wtedy pn " def max |pn(x)| e"21-n. = -1d"xd"1 Dowód (nie wprost) Załóżmy przeciwnie, iż |pn(x)| < 21-n, dla x " [-1, 1]. Niech kĄ qn(x) =21-nTn(x) oraz xk =cos dla k =0, . . . , n. n Z przyjętych założeń oraz z (9.41) wynika, iż (-1)k pn(xk) =|pn(xk)| < 21-n =(-1)kqn(xk), k =0, 1, . . . , n. Wtedy (-1)k[qn(xk) - pn(xk)] > 0 dla k =0, 1, . . . , n. Z tej nierówności wynika, iż wartości wielomianu w(x) =qn(x) - pn(x) wkolej- nych punktach xk, k = 0, 1, . . . , n są różnych znaków. Z twierdzenia Rolle a wynika, że wielomian w = w(x) ma na odcinku (-1, 1) n pierwiastków. Ponie- waż wielomiany pn = pn(x) i qn = qn(x) są n tego stopnia z współczynnikiem przy najwyższej potędze xn równym 1 toteż wielomian w = w(x) jest stopnia n - 1. Doszliśmy w ten sposób do sprzeczności, wielomian w = w(x) jest n - 1 stopnia i ma n pierwiastków.
Niech x0, x1, . . . , xn " [-1, 1] spełniają x0 wie twierdzenia z poprzedniego wykładu mamy dla wielomianu p = p(x) stopnia co najwyżej n następujące oszacowanie
n
1
max |f(x) - p(x)| d" max |f(n)(x)| max (x - xk) (9.42)
-1d"xd"1 (n +1)! -1d"xd"1 -1d"xd"1 k=0 Z twierdzenia 9.3 mamy następujące oszacowanie (dla wielomianu n +1 stop- nia)
n
1
max (x - xk) e" . (9.43)
-1d"xd"1 k=0 2n Ponieważ wielomian o współczynniku przy najwyższej potędze równym 1 jest jednoznacznie wyznaczony przez swoje pierwiastki, to z (9.41) dla
2k +1 xk =cos Ą , k =0, . . . , n 2n +2 9.L10-met 45 Krzysztof Tabisz 25 marca 2003, 8:39am zachodzi
n n
1 1
(x - xk) = Tn+1(x) i dalej też z (9.41) max (x - xk) = .
2n -1d"xd"1 k=0 2n k=0 W ten sposób zrealizowaliśmy dolne oszacowanie w (9.43). Wstawiając ostatnią równość do (9.42) otrzymaliśmy: 9.4 Twierdzenie Jeżeli f " Cn+1[-1, 1] oraz dla k = 0, 1, . . . n punkty odcinka [-1, 1] xk " [-1, 1] są pierwiastkami Tn+1(xk) =0 n +1 wielomianu Czebyszewa i wielo- mian p = p(x) spełniający f(xk) =p(xk) jest co najwyżej stopnia n, to wtedy zachodzi
1 f(n+1)(t) |f(x) - p(x)| d" max .
2n(n +1)! -1d"td"1
Poniższe dwa przykłady zwracają uwagę, że dla y = f(x) ciągłej, ciąg wielo- mianów iterpolacyjnych nie jest zbieżny jednostajnie f - pn " = max f(x) - pn(x) ad"xd"b 1 Przykład 117 Niech f(z) = dla z " Z zespolonego oraz dla n " N z rozważmy n te pierwiastki z jedności 1, 2, . . . , n. Wtedy dla wielomianu pn-1(z) =zn-1 zachodzi 1 1 n-1 n pn-1(k) =k = k = = f(k), dla k =1, 2, . . . , n k k oraz
1
f - pn-1 " = max |f(z) - pn-1(z)| = max - zn-1 =
|z|=1 |z|=1 z 1 = max |1 - zn| = 2, |z|=1 |z| gdzie max jest realizowane dla zn = -1. 1 Przykład 218 Dla f(x) = , x " [-5, 5] funkcji Rungego, rozważa się x2+1 10k punkty xk = -5+ , k =0, 1, . . . , n. Wtedy19 dla ciągu wielomianów inter- n polacyjnych ciąg norm f - pn " nie jest ograniczony. 9.5 Twierdzenie Jeżeli f : [a, b] - R jest funkcją ciągłą, to istnieje taki ciąg podziałów odcinka [a, b] a d" x(n) 0 1 n iż odpowiadający (pn) ciąg wielomianów iterpolacyjnych Czebyszewa jest zbieżny jednostajnie do funkcji y = f(x), tzn. lim f - pn =0. n " 17 1884, Meray 18 1901, Runge 19 J.F. Epperson, On the Runge example, AMM 4, 1987, 329 341. 46 9.L10-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Dowód Można znalezć w [D.Kincaid & W.Cheney] ss. 342.
9.6 Twierdzenie (aproksymacyjne Weierstrassa) Jeżeli f[a, b] - R jest funkcją ciągła i >0, to istnieje wielomian p = p(x) o własności |f(x) - p(x)| d", dla x " [a, b]. Dowód Wystarczy udowodnić dla x " [0, 1] i pózniej podstawić x = a+t(b-a). Aproksymacja jest dana przez wielomiany Bernstein a20
n
k n (Bnf)(x) = f xk(1 - x)n-k, n k k=0 gdzie ż# n!
# , dla 0 d" k d" n k! (n-k)! n = k # 0, poza tym. Tak określone wielomiany Bn = Bn(x) mają następujące własności: " Bn : C[0, 1] - C[0, 1] dla n " N, " Bn(a f + b g) =a Bnf + b Bng dla a, b " R i f, g " C[0, 1], " 0 d" f =! 0 d" Bnf dla f " C[0, 1]. Dla funkcji h0(x) a" 1 zachodzi
n
n (Bnh0)(x) = 1 xk(1 - x)n-k =(x +(1- x))n =1. k k=0 Dla h1(x) =x mamy (Bnh1)(x) =
n n
k n k n! = xk(1 - x)n-k = xk(1 - x)n-k = n k n k!(n - k)! k=0 k=0
n - 1 = x xk(1 - x)(n-1)-k = x. k k=0 Dla h2(x) =x2 otrzymujemy 2 n
k n (Bnh2)(x) = xk(1 - x)n-k = n k k=0 20 1912, Serge Berstein 9.L10-met 47 Krzysztof Tabisz 25 marca 2003, 8:39am
n n
k2 n! k n - 1 = xk(1 - x)n-k = xk(1 - x)n-k = n2 k!(n - k)! n k - 1 k=0 k=1
n
n - 1 k - 1 1 n - 1 = + xk(1 - x)n-k = n n - 1 n k - 1 k=1
n
n " n - 1 n - 2 x n - 1 x = x2 xk-2(1 - x)n-k + = x2 + ! x2. n k - 2 n n n k=2 Dowód twierdzenia wynika z następującego ogólniejszego twierdzenia.
9.7 Twierdzenie (Bohman Korovkin) Niech ciąg Ln, n =1, 2, 3, . . . spełnia " Ln : C[a, b] - C[a, b] dla n " N, " Ln(af + bg) =aLnf + b Lng dla a, b " R i f, g " C[a, b], " 0 d" f =! 0 d" Lnf dla f " C[a, b], oraz lim Lnf - f " =0 dla f "{1, x, x2}. n " Wtedy lim Lnf - f " =0 dla f " C[a, b]. n " Dowód Operator Ln : C[a, b] - C[a, b] ma własności g d" f =! 0 d" f - g =! 0 d" Ln(f - g) =! 0 d" Lnf - Lng =! Lng d" Lnf. Z własności wartości bezwzględnej | |: |a| d"b !! -b d" a d" b dla a, b " R, wynika |Lnf| d"Ln(|f|) dla f " C[a, b]. Niech h0(x) =1, h1(x) =x i h2(x) =x2. oraz ąn = Lnh0 - h0, n = Lnh1 - h1 i łn = Lnh2 - h2. Z założeń wynika, iż przy n "zachodzi ąn " - 0, n " - 0, łn " - 0. Niech f " C[a, b] i 0 < . Pokażemy, że istnieje takie m " N, iż dla każdego m Lnf - f " < 3. 48 9.L10-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Z ciągłości y = f(x) dla x " [a, b] wynika, że funkcja ta jest jednostajnie ciągła na [a, b], tzn. istnieje taka 0 <, iż dla każdych x, y " [a, b] zachodzi |x - y| < =!|f(x) - f(y)| <. ozn. 2 f " Dla c = zachodzi 2 (x - y)2 d"|x - y| =!|f(x) - f(y)| d"2 f " d" 2 f " = c(x - y)2. 2 W ten sposób dla każdego x, y " [a, b] zachodzi |f(x) - f(y)| d" + c(x - y)2. Nierówność tę można zapisać w postaci |f - f(y)h0| d"h0 + c[h2 - 2yh1 + y2h0]. Z monotoniczności operatora Ln wynika: |Lnf - f(y)Lnh0| d"Lnh0 + c[Lnh2 - 2yLnh1 + y2Lnh0]. Mamy teraz |(Lnf)(y) - f(y)(Lnh0)(y)| d"
= [1 + ąn(y)] + c y2 + łn(y) - 2y(y + n(y)) + y2(1 + ąn(y)) = = + ąn + cłn(y) - 2cyn(y) +cy2ąn(y) d" d" + ąn " + c łn " +2c h1 " n " + c h2 " ąn ". Skąd wynika, że istnieje takie m " N, iż dla każdego m d" n zachodzi Lnf - f Lnh0 " d" . Na koniec oszacujemy Lnf - f " d" Lnf - f Lnh0 " + f Lnh0 - f h0 " d" d" 2 + f " ąn ".
Informacyjnie Historyczne oznaczenia ( [2] s. 51, [3] s. 144 oraz [D.Kincaid & W.Cheney] ss. 350 367) 9.8 Definicja Dla danej funkcji y = f(x) oraz liczby h " R określamy indukcyjnie 9.L10-met 49 Krzysztof Tabisz 25 marca 2003, 8:39am " "kf(x), k tą różnicę progresywną funkcji f(x) ż# # f(x), dla k =0 "kf(x) = # "k-1f(x + h) - "k-1f(x), k =1, 2, . . . " "kf(x), k tą różnicę wsteczną funkcji f(x) ż# # f(x), dla k =0 "kf(x) = # "k-1f(x) -"k-1f(x - h), k =1, 2, . . . " kf(x), k tą różnicę centralną funkcji f(x) ż# f(x), dla k =0 # kf(x) = # 1 1 k-1f(x + h) - k-1f(x - h), k =1, 2, . . . 2 2 " dla x0 ż# f(x0), dla n =0 # # # # # f(x1)-f(x0) ,n =1 f[x0, x1, . . . , xn] = x1-x0 # # # # # f[x1,x2,...,xn]-f[x0,x1,...,xn-1] , n =2, 3, . . . xn-x1 Interpolacja Hermite a ( [2] s. 65), [D.Kincaid & W.Cheney] s. 363) Niech y = f(x) będzie wystarczająco regularna. Dla x0 < x1 < . . . < xn szukamy wielomianu p = p(x) o możliwie najmniejszym stopniu, spełniającego p(j)(xi) =f(j)(xi), j =0, 1, . . . , ki - 1, i =0, 1, . . . , n gdzie m +1=k0 + k1 + . . . + kn jest liczbą równań jakie mają być spełnione. Interpolacja splajnami ( [D.Kincaid & W.Cheney] ss. 374 414) Dla x0 węzły w punktach xi, i =0, . . . , n, gdy spełnia " S(xk) =f(xk) dla k =0, 1, . . . , n
" dla każdego k =1, . . . , n obcięcie S jest wielomianem co najwyżej
[xk-1,xk) stopnia k " S " Ck-1[x0, xn] Zadania na Ćwiczenia 1. Wyznaczyć wielomian interpolacyjny najmniejszego stopnia realizującu: 50 9.L10-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 x 3 7 (a) y 5 -1 x 7 1 2 (b) y 146 2 1 x 3 7 1 2 (c) y 10 146 2 1 x 3 7 1 2 (d) y 12 146 2 1 x 1.5 2.7 3.1 -2.1 -6.6 11.0 (e) y 0.0 0.0 0.0 1.0 0.0 0.0 2. Udowodnić, że dla 23 dowolnie wybranych punktów z odcinka [-1, 1] i wie- lomianu interpolacyjnego p = p(x) stopnia 22 oraz funkcji f(x) =cosh x zachodzi |p(x) - f(x)| sup d" 5 10-16. |f(x)| x"[-1,1] 3. Jeżeli p = p(x) jest wielomianem iterpolacyjnym funkcji f(x) = sinh x stopnia nie większego od n - 1 oraz jeden z punktów interpolacyjnych jest w 0, to 2n |p(x) - f(x)| d" |f(x)| dla x " [-1, 1]. n! 4. Dla wielomianów Czebyszewa wykazać równość Tn(x) =cosh(n cosh-1 x). 5. Współczynnik przy xn w wielomianie Tn(x) wynosi 2n-1. Wyliczyć współ- czynniki przy xn-1 oraz xn-2. 6. Wyznaczyć wielomian interpolacyjny Lagrange a i Newton a dla następu- jących danych x 2 0 3 x -2 0 1 a) b) f(x) 11 7 28 f(x) 0 1 -1 7. Udowodnić Ą# ń# 1 x0 x2 . . . xn 0 0 ó# Ą# 1 x1 x2 . . . xn 1 1 ó# Ą#
ó# Ą# 1 x2 x2 . . . xn 2 2 det ó# Ą# = (xk - xj). ó# Ą# . . . . . . . . . . 0d"jŁ# Ś# . . . . . 1 xn x2 . . . xn n n 8. Cena znaczka pocztowego w USA była następująca (w centach) rok 1885 1917 1919 1932 1958 1963 1968 1971 cena 2 3 2 3 4 5 6 8 9.L10-met 51 Krzysztof Tabisz 25 marca 2003, 8:39am rok 1974 1978 III.1981 X.1981 1985 1988 1991 1995 cena 10 15 18 20 22 25 29 32 Wyznaczyć wielomian interpolacyjny Newton a dla tych danych. Kiedy cena znaczka pocztowego wyniesie $1 i $10? 9. Jeżeli y = f(x) jest wielomianem stopnia k, to dla n>k zachodzi f[x0, x1, . . . , xn] =0. 10. Ze wzorów Cramer a wyprowadzić
1 xn . . . xn-1 f(xn) 1 xn . . . xn n n 11. Dla f(x) =xm, m " N, wykazać
1, dla n = m f[x0, x1, . . . , xn] = 0, dla n>m 12. Wykazać n
(fg)[x0, x1, . . . , xn] = f[x0, x1, . . . , xk]g[xk, xk+1, . . . , xn]. k=0 Zadania na Laboratoria 1 1. Dla n = 5, 10 i 15 wyznaczyć dla funkcji f(x) = wielomian iter- 1+x2 polacyjny Newton a pn = pn(x) na przedziale [-5, 5]. Przyjąć następujące 10k punkty podziału xk = -5+ , k =0, 1, . . . , n. Celem oceny zbieżności n wielomianów pn do funkcji f obliczyć 10i f(xi) - pn(xi) gdzie xi = -5+ dla i =0, 1, . . . , 30. 30 Literatura [D.Kincaid & W.Cheney] D. Kincaid & W. Cheney, Numerical Analysis, Brooks 1996. [A.Ralston] A. Ralston, Wtęp do analizy numerycznej, PWN, Warszawa 1983. [G.W.Stewart] G.W. Stewart, Afternotes on Numerical Analysis, SIAM, Phila- delphia 1996. 52 9.L10-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 10 Wykład 12 Wykład 12 Numeryczne całkowanie Przykład Zadanie polega na wyznaczeniu wartości funkcji: x 2 2 " erf(x) = e-t dt. Ą 0 " xn " 2 Widomo, iż ex = co po podstawieniu -x2 daje e-x = (-1)n x2n . n=0 n! n=0 n! Otrzymany szereg potęgowy jest zbieżny dla x " R oraz dla każdego x " R na odcinku [0, x] zbieżny jest jednostajnie. Stąd mamy erf(x) = x x x " "
2 x2n+1 = " (-1)n . Ą (2n +1) n! n=0 Dla małych x (np. |x| d"1) otrzymany szereg jest dobrze zbieżny i może służyć jako dobre przybliżenie całki. Przypomnienie (Całka Riemann a) Dla danego odcinka [a, b] " R i funkcji f : [a, b] - R określamy ": a = x1 (") = max (xi - xi-1) średnica podziały " odcinka [a, b] 1d"id"n Określenie całki Riemann a b n
def f(x) dx = lim f(i)(xi - xi-1), gdzie i " (xi-1, xi) (")0 i=1 a oraz granica jest brana po dowolnych podziałach " odcinka [a, b]. Jeżeli f : [a, b] - R jest funkcją ograniczoną, to warunkiem koniecznym i wystarczającym istnienia całki Riemann a jest aby zbiór punktów nieciągłości funkcji y = f(x) miał miarę (zewnętrzną Lebegue a) 0, tzn. ({x0 " [a, b] : lim f(x) = f(x0)}) =0.
xx0 Jeżeli f : [a, b] - R jest funkcją ciągłą, to b n
1 f(x) dz = lim f(i), n " n i=1 a b-a b-a gdzie i " [a + (i - 1), a+ i] dla i =1, 2, . . . , n. n n 10.L11-met 53 Krzysztof Tabisz 25 marca 2003, 8:39am Przykłady n 1 1 2 1 1. e-x dx = lim e- i2 , n n " 0 i=1 1 n 1
1 dx 2. lim = =ln(1 +x) =ln 2.
i+n 1+x n " 0 i=1 0 (Zasadnicze twierdzenie rachunku całkowego) Jeżeli f : [a, b] - R jest funkcją ciągłą, to b b ozn. f(x) dx = F (b) - F (a) = F (x) ,
a a
gdzie funkcja y = F (x) spełnia F (x) =f(x).
Dla danej funkcji y = f(x) funkcje y = F (x) o własności F (x) =f(x) nazy- wamy całką nieoznaczoną (lub całką Newton a) funkcji y = f(x) i oznaczamy w następujący sposób
F (x) = f(t) dt. Całka nieoznaczona jest wyznaczona z dokładnością do stałej. Własności całki Jeżeli f, g : [a, b] - R są funkcjami całkowalnymi i ą, " R to zachodzi b b b ąf(x) +g(x) dx = ą f(x) dx + g(x) dx a a a oraz b b | f(x) dx| d" |f(x)| dx. a a Jeżeli ciąg funkcyjny fn : [a, b] - R funkcji całkowalnych jest zbieżny jedno- stajnie do funkcji f : [a, b] - R całkowalnej, to b b lim fn(x) dx = lim fn(x) dx. n " n " a a Jeżeli f : [a, b] - R jest całkowalna oraz ": a = x0 dowolnym podziałem odcinka [a, b], to i b x n
f(t) dt = f(t) dt. i=1xi-1 a Jeżeli funkcja f : [a, b] - R jest całkowalna oraz w punkcie x0 " (a, b) jest ciągła, to funkcja x F (x) = f(t) dt a 54 10.L11-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003
jest w punkcie x0 różniczkowalna i zachodzi związek F (x0) =f(x0). Metoda trapezów Jeżeli f : [a, b] - R, to b f(x) dx H" a b b f(b) - f(a) f(b) - f(a) (x - a)2 H" (x - a) +f(a) dx = + f(a)x = b - a b - a 2 a a f(b) - f(a) (b - a)2 b - a = + f(a)(b - a) = [f(a) +f(b)]. b - a 2 2 Jeżeli ": a = x0 i b x n n
xi - xi-1 f(t) dt = f(t) dt H" [f(xi-1) +f(xi)]. 2 i=1xi-1 i=1 a Jeżeli funkcja f : [a, b] - R jest ciągła, to funkcja y = f(x) jest jednostajnie ciągła na odcinku [a, b] i ciąg łamanych f(xi) - f(xi-1) fn(x) = (x-xi-1)+f(xi), dla x " [xi-1, xi] oraz i =1, 2, . . . , n xi - xi-1 jest jednostajnie zbieżny do funkcji y = f(x) o ile (") 0. Przykładowo, dla funkcji y = f(x) ciągłej, podział " można określić w nastę- pujący sposób: b - a ": xi = a + i, gdzie i =0, 1, 2, . . . , n n i wtedy otrzymamy b n
1 1 x0 f(xn) f(t) dt H" [f(xi-1)+f(xi)] = + f(x1) +. . . + f(xn-1) + . n n 2 2 i=1 a Dalej będziemy rozważali tylko podziały " odcinka [a, b] na n równych podod- cinków a o funkcji y = f(x) będziemy zakładali, że jej druga pochodna jest ciągła. Oszacowanie błedy w metodzie trapezów Dla danego n " N niech f(xi) def - f(xi-1) b - a i(x) = (x-xi-1)+f(xi-1), xi = i+a dla i =0, 1, . . . , n xi - xi-1 n oraz A = max |f (t)|. ad"td"b 10.L11-met 55 Krzysztof Tabisz 25 marca 2003, 8:39am Wtedy dla danego i = 1, 2, . . . , n oraz x " [xi-1, xi] z twierdzenia o wartości średniej Lagrange a wynika |f(x) - i(x)| = f(xi) - f(xi-1) = |f(x) - (x - xi-1) - f(xi-1)| = xi - xi-1 = |f (2)(x - xi-1) - f (1)(x - xi-1)| = |f ()(2 - 1)(x - xi-1)| d" A d" (x - xi-1) n gdzie 2 " (xi-1, x), 1 " (xi-1, xi) oraz " (1, 2). Skąd wynika, iż i b x n
| f(t) dt - i(t) dt| d" i=1xi-1 a # ś# i i i x x x n n
# d" | f(t) dt - i(t) dt # | = | f(t) - i(t) dt| d" i=1 i=1xi-1 xi-1 xi-1 i i x x n n n
i A A (t - xi-1)2 x
d" |f(t) - i(t)| dt d" t - xi-1 dt = = n n 2 xi-1 i=1xi-1 i=1 i=1 xi-1 n n
A (xi - xi-1)2 A 1 A = = = . n 2 n 2n2 2n2 i=1 i=1 56 10.L11-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 11 Wykład 13 Wykład 13 Numeryczne całkowanie cd Jeszcze raz oszacowanie błędu w metodzie trapezów Niech f : [a, b] - R posiada drugą pochodną ciągłą oraz dla n " N dzielimy odcinek [a, b] na n równych części, tzn. b - a ": xi = a + i, gdzie i =0, 1, 2, . . . , n. n W metodzie trapezów funkcję y = f(x) na każdym odcinku [xi-1, xi], i = 1, 2, . . . , n interpolujemy prostą, tzn. f(xi) def - f(xi-1) i(x) = (x - xi-1) +f(xi-1) dla i =1, 2, . . . , n, xi - xi-1 oraz n
1, x " A gdzie A(x) = funkcja charakterystyczna zbioru A. Wtedy otrzy- 0, x " A
mujemy metodę na przybliżenie całki b b 1 f(x0) f(xn) f(x) dx H" fn(x) dx = + f(x1) +. . . + f(xn-1) + . n 2 2 a a 11.1 Twierdzenie Jeżeli f : [a, b] - R ma drugą pochodną y = f (x) ciągłą, to istnieje " [a, b] o własności b b (b - a)3f () f(x) dx - fn(x) dx = - . 12n2 a a Dowód W pierwszej kolejności oszacujemy błąd jaki popełniamy dla pojedynczego od- cinka [xi-1, xi] z podziału" interpolując funkcję y = f(x) prostą (Tw 2, wykład 12). Dla x " (xi-1, xi) określamy (t) =f(t) - i(t) - (t - xi-1)(t - xi) gdzie stała jest wyznaczona jest przez równość (x) =0. Otrzymujemy stąd, iż f(x) - i(x) = . (x - xi-1)(x - xi) Dla tak określonego funkcja y = (t) ma na odcinku [a, b] trzy pierwiastki (w punktach xi-1, x, xi). Oczywiście y = (t) jest tak regularna jak y = f(t) tzn. jej druga pochodna jest ciągła. 11.L12-met 57 Krzysztof Tabisz 25 marca 2003, 8:39am Teraz z twierdzenia Rolle a otrzymujemy, że istnieją 1 " (xi-1, x), 2 " (x, xi) o własności (1) = (2) = 0. Stosując jeszcze raz twierdzenie Rolle a do pochodnej y = (t) rozważanej na odcinku [1, 2] stwierdzamy, iż istnieje i " (1, 2) dla której (i) =0. Wyliczamy teraz drugą pochodną y = (t): f(x) - i(x) (t) =f (t) - 2 . (x - xi-1)(x - xi) Teraz z warunku (i) =0 otrzymujemy 1 f(x) - i(x) = f (i)(x - xi-1)(x - xi). 2 Funkcja x (x - xi-1)(x - xi) jest stałego znaku na odcinku [xi-1, xi]. Punkt
f(x)- i(x) i zależy od x i funkcja x f (i) =2(x-xi-1)(x-xi) jest ciągła na [xi-1, xi].
Z twierdzenia o wartości średniej dla rachunku całkowego otrzymujemy i x f(x) - i(x) dx = xi-1 i i x x 1 1 = f (i)(x - xi-1)(x - xi) dx = f (i) (x - xi-1)(x - xi) dx = 2 2 xi-1 xi-1 h h 1 1 t3 t2 (xi - xi-1)3 = f (i) t(t - h) dt = f (i) - h = - f (i) = 2 2 3 2 12 0 0 (b - a)3 = - f (i) 12n3 b-a gdzie skorzystaliśmy z podstawienia t := x - a, h := b - a oraz xi - xi-1 = . n Na koniec mamy b b f(x) dx - fn(x) dx = a a i x n n
(b - a)3 1 (b - a)3 = f(x) - i(x) dx = - f (i) = - f () 12n2 n 12n2 i=1xi-1 i=1 n
1 gdzie min f (x) d" f (i) d" max f (x) i z twierdzenia Darboux wynika, n x"[a,b] x"[a,b] i=1 n
1 iż istnieje " [a, b] dla której f() = f (i). n i=1
Uogólnienie metody trapezów, metoda Newton Cotes. W metodzie trapezów dla danego n " N odcinek [a, b] został podzielony na n 58 11.L12-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 b-a odcinków o długości równej . Na każdym odcinku funkcję y = f(x) interpo- n lowaliśmy prostą. W metodzie Newton Cotes na każdym odcinku interpolujemy funkcję wielomianem Lagrange a stopnia nie większego niż n. Wielomian interpolacyjny Lagrange a (wykład 11), dla ustalonych punktów x(i), x(i), . . . , x(i) " m 0 1 [xi-1, xi] jest postaci m
fi,m(x) = f(x(i)) i,k(x) k k=0 gdzie m
x - x(i) j i,k(x) = . x(i) j=0 - x(i) k j j =k Stąd otrzymujemy, iż i i i x x x m
b b - a a + b f(x) dx H" f(a) +4f + f(b) . 6 2 a Oszacujemy teraz błąd. 11.L12-met 59 Krzysztof Tabisz 25 marca 2003, 8:39am x
ozn b-a Niech h = i F (x) = f(t) dt. Wtedy 2 a a+2h
h F (a +2h) = f(x) dx H" [f(a) +4f(a + h) +f(a +2h)] = (11.44) 3 a 4 2 100 = 2hf(a) +2h2f (a) + h3f (a) + h4f (a) + h5f(4)(a) +. . . 3 3 3 5! Skorzystaliśmy tu ze wzoru Taylor a f (a) f (a) f (a) f(4)(a) f(a + h) = f(a) + h + h2 + h3 + h4 + . . . 1! 2! 3! 4!
Z podstawowego wzoru rachunku całkowego F (x) =f(x) i wzoru Taylor a za- stosowanego do F (a +2h) (oczywiście F (a) =0 !) otrzymujemy 4 2 32 F (a +2h) = 2hf(a) +2h2f (a) + h3f (a) + h4f (a) + h5f(4)(a) +. . . 3 3 5! Porównując otrzymane wyrażenie z prawą stroną przybliżenia (11.44) otrzymu- jemy a+2h
h 1 f(x) dx - [f(a) +4f(a + h) +f(a +2h)] = - h5f(4)(a) +. . . (11.45) 3 90 a co oktreśla błąd metody dla pojedynczego przedziału podziału. Teraz podobnie jak poprzednio dzielimy odcinek [a, b] na n " N (n parzyste) liczbę podocinków do każdego z nich stosujemy poprzednie rozumowanie. Ozna- czmy b - a xi = a + ih h = (0 d" i d" n). n Wtedy mamy x2i b n/2
h Ł#f(x0) +2 f(x2i-2) +4 f(x2i-1) +f(xn)Ś# . = 3 i=2 i=1 Jeżeli do (11.46) wstawimy (11.45) to otrzymamy następujący wzór na błąd metody (b - a)5f(4)() - gdzie " (a, b). 180n4 60 11.L12-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 12 Wykład 14 Wykład 14 Numeryczne całkowanie zakończenie Aby określić całkę Riemann a b f(x) dx a o funkcji f : [a, b] - R zakładamy, iż jest ograniczona. Całkę osobliwą z funkcji nieograniczonej na przedziale całkowania określa się w następujący sposób: 1 1 1 " " dx dx def " = lim " = lim x = lim (1 - ) = 1.
0+ 0+ 0+ 2 x 2 x 0 W przypadku obliczeń numerycznych całek z funkcji typu c <" f(x) " po przedziale (0, 1] = x postępujemy w następujący sposób: Wprowadzamy nową funkcję " def g(x) = xf(x) która na przedziale (0, 1] zachowuje się dobrze tzn. jest ciągła na całym [0, 1]. Stsujemy teraz jedną z poznanych metod całkowania numerycznego np. metodę 1 3 nieokreślonych współczynników zastosowanej do punktów x0 = , x1 = . Dla 4 4 g(x) = 1 otrzymujemy następujące równania na współczynniki A0, A1: 1 1 " 1 dx = 2 = A0 + A1 x 0 biorąc za g(x) = x dostajemy drugie równanie 1 1 2 1 3 " x dx = = A0 + A1. x 3 4 4 0 Wyliczając A0, A1 otrzymujemy przybliżoną formułę całkowania numerycznego 1 1 5 1 2 3 <" g(x) " dx g( ) + g( ). = x 3 4 3 4 0 Zmieniając odcinek całkowania z [0, 1] na [0, h] (przekształcenie: x = hy) otrzymujemy wzór h 1 " 1 1 g(x) " dx = x g(hx) " dx. x x 0 0 12.L13-met 61 Krzysztof Tabisz 25 marca 2003, 8:39am Aby porównać efektywność opisanej metody porównamy ją z metodą Newton Cotes zastosowaną do tych samych punktów dla wzór przybliżony jest następu- jący h h h 3h <" f(x) dx = f( ) +f( ) . 2 4 4 0 Przykład21 Niech h = 0.1 oraz funkcji " cos x " f(x) = - x sin x 2 x której całka nieoznaczona wynosi
" " cos x " - x sin xdx = x cos x. 2 x 21 za: G.W. Stewart, Afternotes on Numerical Analysis, ss. 168, SIAM, 1966. 62 12.L13-met 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 13 Metoda gradientów 13.L14-met 63 Krzysztof Tabisz 25 marca 2003, 8:39am A Uzupełnienia A.1 Analiza A.1 Definicja (zbiór ograniczony) Niepusty zbiór A " R nazywamy ograniczonym, gdy "M "a"A |a| M. A.2 Twierdzenie Uzupełnienia z analizy A.2 Algebra Uzupełnienia z algebry Notatki do wykładu z dnia 11.X.02r A.2.1 Struktura form dwuliniowych Opracowano na podstawie [1, Roz. XIV, ss. 382 409]. Niech E będzie przestrze- nią liniową nad ciałem liczb zespolonych C z iloczynem skalarnym ć% : E E - C. Pojęcia wstępne A.3 Definicja (forma symetryczna) Niech E będzie przestrzenią liniową na ciałem rzeczywistym R. Wtedy dwuli- niowe odwzorowanie g : E E - R nazywamy formą symetrycznym, gdy zachodzi g(x, y) =g(y, x) dla x, y " E. A.4 Definicja (forma hermitowska) Niech E będzie przestrzenią liniową nad ciałem C liczb zespolonych. Wtedy odwzorowanie g : E E - C jest formą hermitowską, gdy jest liniowe względem swego pierwszego argumentu, a antyliniowe względem drugiego, oraz g(x, y) =g(y, x). ozn Będziemy pisali g(x, y) = x, y , jeżeli wiadomo o jaką formę g chodzi. Tak więc dla formy hermitowskiej zachodzi x1 + x2, y = x1, y + x2, y dla , " C i x1, x2, y " E, x, y1 + y2 = x, y1 + x, y2 dla , " C i x, y1, y2 " E, x, y = x, y dla x, y " E. 64 A.AppMet1 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 A.5 Definicja (ortogonalność wektora do podprzestrzeni) Wektor x " E jest ortogonalny do podprzestrzeni F " E względem formy hermitowskiej , , co zapisujemy x Ą" F , gdy def x Ą" F !! "y"F x, y =0. A.6 Definicja (podprzestrzenie ortogonalne)
Podprzestrzenie F, F " E są ortogonale względem formy hermitowskiej , ,
co zapisujemy F Ą" F , gdy def
F Ą" F !! "y"F "y "F y, y =0. A.7 Definicja (niezdegenerowana forma hermitowska) Forma hermitowska , jest niezdegenerowana, gdy zachodzi ("y"E x, y =0) =! x =0 dla x " E. A.8 Uwaga Z Definicji A.5 ortogonalności wektora do podprzestrzeni względem formy hermito- wskiej , wynika, że Definicja A.7 niezdegenerowania formy hermitowskiej jest równoważna warunkowi: x Ą" E =! x =0 dla x " E. A.9 Definicja (suma prosta podprzestrzeni ortogonalnych)
Przestrzeń E jest sumą prostą podprzestrzeni ortogonalnych F, F , co zapisu- jemy
E = F Ą" F gdy spełnione są warunki
E = F " F ,
F Ą" F . A.10 Definicja (jądro formy) Jeżeli , jest formą symetryczną lub hermitowską, to zbiór E0 = {x " E : x Ą" E} nazywamy jądrem formy. A.11 Definicja (baza ortogonalna) Baza przestrzeni E {v1, v2, v3, . . . , vn} jest bazą ortogonalną jeżeli vi, vj =0 dla i = j.
A.AppMet1 65 Krzysztof Tabisz 25 marca 2003, 8:39am A.12 Stwierdzenie Niech E będzie przestrzenią wektorową nad ciałem R lub C, a g niech będzie formą symetryczną lub hermitowską. Załóżmy, że przestrzeń E ma rozkład or- togonalny E = E1 Ą" . . . Ą" Em. Wówczas forma g nie jest zdegenerowana na przestrzeni E wtedy i tylko wtedy, 0 gdy nie jest zdegenerowana na żadnej z przestrzeni Ei. Jeżeli Ei jest jądrem ograniczenia formy g do Ei, to jądrem formy g na przestrzeni E jest suma ortogonalna 0 0 E0 = E1 Ą" . . . Ą" Em. A.13 Stwierdzenie Niech E będzie przestrzenią skończenie wymiarową nad ciałem R lub C, a g formą symetryczną lub hermitowską na przestrzeni E. Załóżmy, że g jest formą niezdegenerowaną. Niech F będzie podprzestrzenią przestrzeni E. Forma g jest Ą" niezdegenerowana na F wtedy i tylko wtedy, gdy F + F = E, a także wtedy i Ą" tylko wtedy, gdy jest niezdegenerowana na F . Formy symetryczne A.14 Twierdzenie Niech E będzie przestrzenią wektorową nad ciałem R, a g niech będzie formą symetryczną określoną na przestrzeni E. Jeżeli dim E 1, to w przestrzeni E istnieje baza ortogonalna. A.15 Wniosek Niech {v1, . . . , vn} będzie bazą ortogonalną przestrzeni E. Załóżmy, że g(vi, vi) = 2 2 vi dla i r oraz vi = 0 dla i > r. Wówczas jądro przestrzeni E równa się (vr+1, . . . , vn). A.16 Twierdzenie (Sylvester) Niech E będzie przestrzenią wektorową nad ciałem R, z niezdegenerowaną formą symetryczną g. Istnieje taka liczba całkowita nieujemna r, że dla dowolnej bazy 2 2 ortogonalnej {v1, . . . , vn} przestrzeni E wśród elementów v1, . . . , vn dokładnie r jest dodatnich i dokładnie n - r ujemnych. A.17 Wniosek 2 Istnieje taka baza ortogonalna {v1, . . . , vn} przestrzeni E, że vi =1 dla i r 2 oraz vi = -1 dla i >r, przy czym liczba r jest wyznaczona jednoznacznie. A.18 Definicja (dodatnia określoność) Forma symetryczna g jest dodatnio określona, jeżeli x2 = g(x, x) > 0 dla każdego x " E, x =0.
A.19 Wniosek Przestrzeń wektorowa E ma taki rozkład ortogonalny E = E+ +E-, że forma g jest dodatnio określona na E+ i ujemnie określona na E-. Wymiar przestrzeni E+ (lub E-) jest taki sam we wszystkich takich rozkładach. 66 A.AppMet1 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Ortogonalizacja Grama Schmidta Dla formy dodatnio określonej, istnieje kanoniczna droga otrzymania bazy or- tonormalnej poprzez proces indukcyjny rozpoczynający się z dowolnej bazy {v1, . . . , vn}. Niech 1
Wówczas {v1, . . . , vn} jest bazą ortonormalną. A.20 Lemat (tożsamość polaryzacyjna) Niech A : E - E będzie przekształceniem liniowym. Wtedy, jeżeli Ax, y =0 dla każdego x " E, to A =0. Dalej zakładamy, że dim E < ". Niech A : E - E będzie odwzorowaniem liniowym. Dla ustalonego y " E odwzorowanie x - Ax, y jest funkcjonałem
liniowyn, a więc istnieje jednoznacznie określony element y" " E, taki że Ax, y = x, y" dla wszystkich x " E. A.21 Definicja (odwzorowanie sprzężone) Określamy odwzorowanie A" : E - E, przyjmując A"y = y". A.22 Definicja Odwzorowanie liniowe A" nazywamy odwzorowaniem samosprzężonym (albo odwzorowaniem hermitowskim), jeżeli A = A". A.23 Stwierdzenie Odwzorowanie A jest wtedy i tylko wtedy hermitowskie, gdy forma Ax, x jest rzeczywista dla każdego x " E. A.24 Stwierdzenie Niech odwzorowanie A będzie hermitowskie. Wówczas wszystkie wartości wła- sne odwzorowania A są rzeczywiste. Jeżeli , są wektorami własnymi =0, o
wartościach własnych odpowiednio , , to dla = jest Ą" .
A.AppMet1 67 Krzysztof Tabisz 25 marca 2003, 8:39am A.25 Lemat Nich A : E - E będzie odwzorowaniem liniowym oraz niech dim E 1. Wówczas odwzorowanie A ma co najmniej jeden wektor własny, różny od zera. A.26 Twierdzenie (spektralne przypadek hermitowski) Niech E będzie niezerową przestrzenią wektorową nad ciałem liczb zespolonych z dodatnio określoną formą hermitowską. Niech A : E - E będzie odwzoro- waniem liniowym hermitowskim. Wówczas przestrzeń E ma bazę ortogonalną, złożoną z wektorów własnych odwzorowania A. A.27 Wniosek Przy założeniach twierdzenia istnieje baza ortonormalna, złożona z wektorów własnych odwzorowania A. A.28 Wniosek Niech E będzie niezerową przestrzenią wektorów nad ciałem liczb zespolonych z dodatnio określoną formą hermitowską f. Niech g będzie drugą formą hermi- towską na przestrzeni E. Wówczas istnieje baza przestrzeni E, ortogonalna dla formy f oraz dla formy g. A.29 Definicja (odwzorowanie sprzężone) Teraz niech E będzie przestrzenią wektorową nad ciałem liczb rzeczywistych, a g symetryczną dodatnio określoną formą na przestrzeni E. Jeżeli A : E - E jest odwzorowaniem liniowym, to odwzorowanie AT , sprzężone do A względem formy g jest określone wzorem Ax, y = x, AT y dla x, y " E. A.30 Twierdzenie (spektralne przypadek symetryczny) Nich E będzie niezerową przestrzenią wektorową nad ciałem liczb rzeczywistych, z formą symetryczną określoną dodatnio. Niech A : E - E będzie symetrycz- nym odwzorowaniem liniowym. Wówczas przestrzeń E ma bazę ortogonalną składającą się z wektorów własnych odwzorowania A. A.3 Analiza funkcjonalna Uzupełnienia z analizy funkcjonalnej 68 A.AppMet1 25 marca 2003, 8:39am Metody numeryczne I, 2003/2003 Literatura [1] S. Lang, Algebra, PWN, Warszawa, 1973. [2] A. Ralston, Wstęp do analizy numerycznej, PWN, Warszawa, 1983, Tłumaczył R. Zuber, wydanie polskie opracował St. Paszkowski. Do wydania polskiego jest dedykacja autora. [3] G.W Stewart, Afternotes on numerical analysis, SIAM, 1996, Polecona przez A. Szustalewicza. LitMet 69 Krzysztof Tabisz 25 marca 2003, 8:39am Wrocław, dnia 25 marca 2003 ( ) Krzysztof Tabisz tabisz@math.uni.wroc.pl 70 Indeks