(wykład do wyboru)
Krzysztof Tabisz
Instytut Matematyczny
Uniwersytet Wrocławski
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 . . . 64A.2.1 Struktura form dwuliniowych . . . 64
A.3 Analiza funkcjonalna . . . 68
Literatura 69
Wstęp
1
Wykład 1
Wykład 1 Postawowe fakty z ANALIZY: • A ⊂ R – jest ograniczony, gdy ∃M∀a∈A|a| ≤ M ;
• Niech A ⊂ R ograniczony i A 6=. Wtedy α = sup A, gdy – ∀a∈Aa ≤ α,
– (∀a∈Aa ≤ β) =⇒ α ≤ β,
• 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 x→x0 f (x) = g gdy – (Cauchy) ∀ε>0∃δ(ε)>0∀x0 < |x − x0| < δ(ε) =⇒ |f (x) − f (x0)| < ε, – (Heine) ∀(xn), xn6=x0 lim n→+∞xn= x0 =⇒lim n→∞f (xn) = g . • f : R −→ R jest ciągła w x0, gdy lim
x→x0
f (x) = f (x0);
• 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)].
• określenie pochodnej y = f (x) w x0: f0(x0) = lim x→x0
f (x)−f (x0)
x−x0 ;
• warunek konieczny ekstremum y = f (x) przyjętego w x0: f0(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)f0(ξ) = 0,
– (Lagrange) ∃ξ∈(a,b)
f (x)−f (x0)
x−x0 = f
0(ξ),
– (wzór Taylor’a) jeżeli dodatkowo f : (a, b) −→ R jest (n + 1)–razy różniczkowalna, to wtedy ∀x,x0∈(a,b)∃ξ∈(x0,x) f (x) = n X k=0 f(k)(x0) k! (x − x0) k+ E n(x), gdzie En(x) = f(n+1)(ξ) (n + 1)! (x − x0) n+1
Przykład strzelanie z armaty na odległość d:
• model matematyczny:
y00(t) = −g, y(0) = 0, y0(0) = V0sin θ
,
• rozwiązanie teoretyczne: T = 2V0sin θ
g ,
• czas lotu pocisku: T = 2V0sin θ
g ,
• przebyta droga przez pocisk: 2V02cos θ sin θ
g = d,
• określenie nieliniowego równania: f (θ) ≡ 2V
2
0 cos θ sin θ
g − d = 0,
• 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. Pokazć, że |x2− 4| < ε dla 0 < |x − 2| < ε(5 + ε)−1 oraz korzystając tylko
z definicji udowodnić: lim
x→2x 2= 4.
2. Pokazć, że lim
x→1(4x + 2) = 6 wprost z ε − −δ definicji.
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.
5. Pokazać, że funkcja f (x) = x sin1x z f (0) = 0 jest ciągła w 0 ale nie jest 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ędwzględ-nęj jaką można uzyskać przy wykonywaniu obliczeń przy zastosowaniu typów rzeczywistych Real, 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;
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
xn+1= xn− f (xn) f0(x n) , n ≥ 0 to lim n→∞xn = r oraz |xn+1− r| ≤ C(xn− r)2, n ≥ 0. Dowód:
r jest jednokrotnym zerem funkcji y = f (x), gdy f (r) = 0 oraz f0(r) 6= 0. Niech
en= xn− r, oznacza błąd przybliżenia. Wtedy mamy en+1= xn+1− r = xn− f (xn) f0(x n) − r (2.1) = en− f (xn) f0(x n) = enf 0(x n) − f (xn) f0(x n)
Ze wzoru Taylor’a mamy
0 = f (r) = f (xn− en) = f (xn) − enf0(xn) + 1 2e 2 nf00(ξn), gdzie ξn ∈ (r, xn). Skąd otrzymujemy enf0(xn) − f (xn) = 1 2f 00(ξ n)e2n. Co wstawiając do (1) otrzymujemy en+1= 1 2 f00(ξn) f0(x n) e2n≈ 1 2 f00(r) f0(r)e 2 n= Ce2n. (2.2) Dla C ≈ 1 i en ≈ 10−4 z (2) mamy: en+1 ≈ 10−8, en+2 ≈ 10−16. . . . W
ten sposób błąd en+1 jest w przybliżeniu równy e2n co oznacza, iż mamy do
czynienia z tzw. kwadratową zbieżnością. Przejdźmy teraz do formalnego dowodu zbieżności metody. Niech
c(δ)def= 1 2 max |x−r|≤δ|f 00(x)| min |x−r|≤δ|f 0(x)|, δ > 0. (2.3)
Ponieważ f0(r) 6= 0 oraz y = f00(x) jest funkcją ciągłą, to lim
δ→0δc(δ) = 0. Więc
x0, spełniającego |x0− r| ≤ δ, to wtedy |e0| ≤ δ i |ξ0− r| ≤ δ i z określenia c(δ) otrzymujemy 1 2 f00(ξ0) f0(x 0) ≤ c(δ). Tym samym z (3) mamy
|x1− r| = |e1| ≤ e20c(δ) = |e0||e0|c(δ) ≤ |e0|δc(δ) = |e0|% < |e0| ≤ δ.
To pokazuje, że kolejny punkt x1jest w δ otoczeniu r i iteracja może być
kon-tynuowana |e1| ≤ %|e0| |e2| ≤ %|e1| ≤ %2|e0| |e3| ≤ %|e2| ≤ %3|e0| .. . Ogólnie mamy |en| ≤ %n|e0|.
Poniewż 0 ≤ % < 1, to mamy lim
n→∞%
n= 0 i lim
n→∞en = 0.
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 xn+1= xn− f (xn) f0(x n) , n = 0, 1, 2, . . . przybliżamy f0(xn) ≈ f (xn) − f (xn−1) xn− xn−1 , n = 1, 2, . . . otrzymuje wtedy (metoda siecznych)
xn+1= xn− f (xn) h x n−xn−1 f (xn)−f (xn−1) i =xn−1f (xn)−xnf (xn−1) f (xn)−f (xn−1) , n = 1, 2, 3, . . . Komentarz: Jeśli określimy ϕ(u, v) = u − f (u)h u−v f (u)−f (v) i , dla u 6= v u −ff (u)0(u), dla u = v
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,
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
en+1= xn+1− r =
f (xn)xn−1− f (xn−1)xn
f (xn) − f (xn−1)
− r =f (xn)en−1− f (xn−1)en f (xn) − f (xn−1)
Dzieląc przez enen−1i wstawiajć
xn−xn−1 xn−xn−1 = 1 otrzymujemy en+1= xn− xn−1 f (xn) − f (xn−1) × f (xn) en − f (xn−1) en−1 xn− xn−1 × enen−1 (2.4)
Z wzoru Taylor’a mamy
f (xn) = f (r + en) = f (r) + enf0(r) + 1 2e 2 nf 00(r) + O(e3 n). Ponieważ f (r) = 0, to f (xn) en = f0(r) +1 2enf 00(r) + O(e2 n)
oraz dla wskaźnika n − 1 f (xn−1) en−1 = f0(r) +1 2en−1f 00(r) + O(e2 n−1). Ostatecznie f (xn) en −f (xn−1) en−1 = 1 2(en− en−1)f 00(r) + O(e2 n−1).
czywiście xn− xn−1= en− en−1. Możemy teraz f (xn) en − f (xn−1) en−1 xn− xn−1 ≈ 1 2f 00(r).
Pierwszy czynnik w wyrażeniu (4) xn− xn−1
f (xn) − f (xn−1)
≈ 1
f0(r).
Pokazaliśmy w ten sposób, iż en+1≈
1 2
f00(r)
f0(r)enen−1= Cenen−1. (2.5)
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
Tzn. α–rząd zbieżności metody oznacza, iż lim n→∞ |en+1| A|en|α = 1. Wstawiając |en| ∼ A|en−1|α i |en−1| ∼ (A−1|en|) 1 α do (5), otrzymujemy A1+α1|C|−1∼ |en|1−α+1α.
Ponieważ prawa strona zależy od en, lim
n→∞en = 1, to 1 − α + 1
α = 0, więc
α = 1+√5
2 ≈ 1.62 . W ten sposób otrzymaliśmy, iż szybkość zbieżności metody
siecznych jest superliniowa tzn. lepsza niż liniowa. Z równania 1 + α1 = α otrzymujemy A = |C| 1 1+ 1α = |C|1α = |C|α−1= |C|0.62= f00(r) 2f0(r) 0.62 . Z tak wyliczonym A otrzymujemy
|en+1| ≈ A|en|
1+√5 2 .
Zadania na Ćwiczenia
1. Dla równania tan x1 = 0 znaleźć najmniejszy dodatni punkt startowy me-tody Newton’a, dla którego metoda ta jest jeszcze zbieżna.
2. Metoda Steffensen’a. Niech xn+1= xn−
f (xn)
g(xn)
, gdzie g(x) = f (x + f (x)) − f (x)
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: xn+1= xn− fnfn0 (f0 n)2− fnfn00 2 , gdzie fn= f (xn). Przez podstawienie f −→ √f
f0 wyprowadzić powyższy wzór z metody
New-ton’a.
Zadania na Laboratoria
1. Metodą Newton’a znaleźć rozwiązania równania x = tan x najbliźsze punk-tów 4.5 i 7.7.
2. Znaleźć najmniejsze dodatnie x0 dla którego funkcja f (x) = tan xx−2 osiąga
swoje minium. Zastosować metodę Newton’a do f0.
3. Równanie 2x4+ 24x3+ 61x2− 16x + 1 = 0 ma dwa pierwiastki w otoczeniu
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:
ϕ(u) = u − f (u) f0(u), a dla siecznych ϕ(u, v) = vf (u)−uf (v) f (u)−f (v) , dla u 6= v
u −ff (u)0(u), dla u = v
Pierwiastek r równania (3.6) jest granicą ciągu r = lim
n→∞xn,
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
n→∞xn = r.
Mówimy, że w punkcie r metoda jest rzędu p, jeśli istnieje liczba rzeczywista p ≥ 1 taka, że lim n→∞ |xn+1− r| |xn− r| p = lim n→∞ |en+1| |en| p = C 6= 0.
Stałą C 6= 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-maliśmy p = 1+
√ 5
2 = 1.618 . . . .
Arytmetyka na liczbach zmiennopozycyjnych. 1. Przedstawienie liczb rzeczywistych w komputerze:
• typy rzeczywiste w TP ver. 6.0
typ zakres dokładność
w cyfrach zajętość w bajtach Real 2.9 × 10−39 . . . 1.7 × 1038 11 − 12 6 Single 1.5 × 10−45 . . . 3.4 × 1038 7 − 8 4 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− 1 19 − 20 8
• typy rzeczywiste w TC ver. 3.1
typ zakres dokładność
w cyfrach zajętość w bajtach float 3.4 × 10−38 . . . 3.4 × 1038 7 4 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 ≤ |m| < 1 dla x 6= 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: −∞, +∞,
• wynik nieokreślony (np. dla 0
0): NaN.
• obcięcie („chopping”, „truncation”), • zaokrąglenie („rounding”),
• błąd zaokrąglenia dla a 6= 0 liczby o podstawie 2 („unit roundoff error” lub „machine epsilon”):
|b − a| |a| =
2−n, zaokrąglenie
2−n+1, obcięcie gdzie n– liczba cyfr w mantysie,
fl(a)def= a(1 + ε), gdzie ε = b − a
a , a 6= 0, • 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+ ε2ε1) ≈ x(y + z)(1 + ε1+ ε2) = x(y + z)(1 + ε3)
5. Analiza błędów zaokrągleń (przykład)
Twierdzenie Niech x0, x1, . . . , xnbędą dodatnimi liczbami rzeczywistymi,
które są przedstawione w komputerze z „roundoff error” ε. Wtedy błąd względny przedstawienia
n
P
i=0
xi wynosi (1 + ε)n− 1 ≈ nε.
Zadania na Ćwiczenia do metody „siecznych”
1. W metodzie siecznych wykazać, że jeżeli xn→ q gdy n → ∞ i f0(q) 6= 0,
to wtedy f (q) = 0.
2. Korzystając ze wzoru Taylor’a uzasadnić następujące przybliżenie na f0(x):
f0(x) ≈ k
2f (x + h) − h2f (x + k) + (h2− k2)f (x)
(k − h)kh .
3. Metodę siecznych zastosować do funkcji f (x) = x2− 2 z x0= 0 i x1= 1.
Wyliczyć x2.
4. Jakie jest x2z metody siecznych, gdy x0= 1, x1= 2, f (x0) = 2 oraz f (x1) =
5. Równość asymptotyczną dwóch ciągów określamy w następujący sposób: xn ∼ yn ⇐⇒ lim n→∞ xn yn = 1. Udowodnij, że jeżeli xn ∼ yn, un∼ vn i c 6= 0, to wtedy
a) cxn ∼ cyn, b) xcn ∼ ync, c) xnun ∼ ynvn,
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:
a) sinx
2− 1; b) e
x− tan x; c) x3− 12x2+ 3x + 1.
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)10na [0, 1]; d) x3+ 10−4 na [−0.75, 0.5]; e) x19+ 10−4 na [−0.75, 0.5]; f ) x5na [−1, 10]; g) x9na [−1, 10]; h) xe−x2 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.
1. Przedstawić 35 w postaci 12.a1a2. . . a24. Jaki jest „roundoff error” tego
przedstawienia? 2. Liczbę 2
3 1 − 2
−24 zapisać w postaci „single precision”.
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≥ 2−24Si, to
wtedy |S∗ n− Sn| Sn ≤n − 1 224 . 5. Uzasadnić nierówność x − x∗ x ≤ 1 1 + 224, x 6= 0
4
Wykład 5
Wykład 5 Podstawowe fakty z ALGEBRY LINIOWEJ2
Macierzą A ∈ Rn×m o n wierszach i m kolumnach nazywamy układ nm liczb:
A = a11 a12 . . . a1m a21 a22 . . . a2m .. . ... . .. ... an1 an2 . . . anm
oraz przyjmujemy, iż aij def
= (A)ij. Wektory x ∈ Rm, b ∈ Rn zapisujemy w
postaci x = x1 x2 .. . xm , b = b1 b2 .. . 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 Eioznacza 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
• 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 ∈ Rn×p i B ∈ Rp×m: (AB)ij= p X k=1 aik· bkj.
Uwaga: Istnieją macierze A, B takie, że AB 6= 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 = BTAT.
Definicja: Macierz identycznościowa (jednostkowa) I ∈ Rn×n
Idef= 1 0 0 . . . 0 0 1 0 . . . 0 0 0 1 . . . 0 .. . ... ... . .. ... 0 0 0 . . . 1
Własność macierzy identycznościowej
AI = IA = A, dla A ∈ Rn×n.
Definicja: Lewostronna A ∈ Rn×m i prawostronna B ∈ Rm×n macierz odwrotna jest określona wzorem
AB = I.
Twierdzenie Jeżeli A, B ∈ Rn×n i AB = I, to BA = I.
Określenie Jeżeli dla A ∈ Rn×n (macierzy kwadratowej) istnieje B ∈ Rn×n: AB = I, to macierz 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:
Związek przekształceń liniowych z macierzami
Jeżeli f : Rn −→ Rmjest przekształceniem liniowym oraz (e
1, . . . , en) jest bazą Rn a (e01, . . . , e0m) – bazą R m, to macierz (A)ij = (f (ei))j, tzn. f (ei) = ai1e01+ . . . + aime0m, ma własność: Ax = f (x) dla x ∈ Rn. I odwrotnie, jeżeli A ∈ Rn×m – jest macierzą, to
f (x)def= Ax
jest przekształceniem liniowym f : Rn −→ Rmwyznaczonym przez A.
Definicja: Jeżeli A ∈ Rn×n oraz
Ax = λ · x, gdzie λ ∈ R, x ∈ Rn i x 6= 0
to λ nazywamy wartością własną macierzy A a x – wektorem własnym tej ma-cierzy.
Definicja: Wyznacznika macierzy A ∈ Rn×n:
det(A)def= X
σ∈Sn
sgn(σ) · a1σ(1)· a2σ(2)· . . . · anσ(n)
gdzie Sn – zbiór wszystkich permutacji zbioru {1, 2, . . . , n} a sgn(σ)
jest znakiem permutacji σ.
Twierdzenie: (Cauchy’ego) Dla A, B ∈ Rn×n zachodzi wzór:
det(AB) = det(A) det(B)
Twierdzenie Rozwinięcie Laplace’a wyznacznika macierzy A (może teź być defini-cją indukcyjną wyznacznika): det(A) = n X j=1 aijAij,
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 Abpowstała z A przez
dodanie kolumny wyrazów wolnych b.
Rozwiązanie układu (5.10) dla n = m i det(A) 6= 0 dane jest wzorem (wzory Cramera):
xj =
det(A0j)
gdzie macierz A0jjest utworzona z A przez zastąpienie j–tej kolumny kolumną wyrazów wolnych b.
Definicja: Macierz A ∈ Rn×n 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 ∈ Rn×n następujące warunki są równoważne: • macierz odwrotna A−1 istnieje, tzn. macierz A jest nieosobliwa;
• det(A) 6= 0;
• wiersze macierzy A tworzą bazę Rn;
• kolumny macierzy A tworzą bazę Rn;
• jeżeli Ax = 0, to x = 0;
• odwzorowanie f (x)def= Ax jest przekształceniem „różnowartościowym” f : Rn
−→ Rn;
• odwzorowanie f (x)def= Ax jest przekształceniem „na” f : Rn
−→ Rn;
• ∀b∈Rn∃!x∈RnAx = b;
• A jest iloczynem macierzy elementarnych; • 0 nie jest wartością własną macierzy A.
Łatwe do rozwiązania układy równań. Systemy równań „dolnie trójkąt-nych”: a11 0 0 . . . 0 a21 a22 0 . . . 0 a31 a32 a33 . . . 0 .. . ... ... . .. ... an1 an2 an3 . . . ann x1 x2 x3 .. . xn = b1 b2 b3 .. . bn . (4.9)
Rozwiązanie (5.11) jest określone rekurenyjnie wzorami (oczywiście zakładamy, iż det(A) = a11a22· . . . · ann6= 0): x1 = b1 a11 , x2 = b2− a21x1 a22 , .. .
xi = bi− i−1 P j=1 aijxj aii , .. . xn = bn− n−1 P j=1 anjxj ann .
Co daje następujący algorytm3rozwiązywania (5.11) „trójkątnych dolnie” ukła-dów równań liniowych: for(i=0;i<n;i++) { for(j=0;j<i;j++) b[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 X i=1 i X j=1 1 = n X i=1 i = n(n + 1) 2 ∼ = n 2 2
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. Znaleźć największą liczbę n ∈ N, dla której 1.0 + 2−n 6= 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. Niech x:=1.0/3;. Dla różnych deklaracji x wydrukować 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;
• pokazać, że [1
2, 1] zawiera tę samą ilość liczb co [1, 2] z odstępem ε 2.
5
Wykład 6
Wykład 6
Dla x, y ∈ Rn×1 w przyjętym zapisie kolumnowym wektorów
x = x1 x2 x3 .. . xn , y = y1 y2 y3 .. . yn mamy x ∈ Rn×1, yT ∈ R1×n i xyT ∈ Rn×n xyT = x1 x2 x3 .. . xn y1 y2 y3 . . . yn = x1y1 x1y2 x1y3 . . . x1yn x2y1 x2y2 x2y3 . . . x2yn x3y1 x3y2 x3y3 . . . x3yn .. . ... ... . .. ... xny1 xny2 xny3 . . . xnyn oraz xT y ∈ R1×1 xTy = x1 x2 x3 . . . xn y1 y2 y3 .. . yn = n X i=1 xiyi
Aby uprościć zapis wzorów macierzowych wprowadzono zapis blokowy macierzy, który wyjaśniamy na przykładzie:
A = a11 a12 a13 a14 a15 a16 a17 a21 a22 a23 a24 a25 a26 a27 a31 a32 a33 a34 a35 a36 a37 a41 a42 a43 a44 a45 a46 a47 a51 a52 a53 a54 a55 a56 a57 = A11 A12 A13 A21 A22 A23 gdzie A11= a11 a12 a21 a22 , A12= a13 a14 a15 a23 a24 a25 , itp.
Przy powyższych oznaczeniach dodawanie macierzy można zapisać w postaci: A11 A12 A21 A22 + B11 B12 B21 B22 = A11+ B11 A12+ B12 A21+ B21 A22+ B22
i w podobny sposób iloczyn macierzy: A11 A12 A21 A22 B11 B12 B21 B22 = A11B11+ A12B21 A11B12+ A12B22 A21B11+ A22B21 A21B12+ A22B22 . Niezależnie od wprowadzanych oznaczeń, podziału macierzy A na podmacierze Aij, Bij– operacje macierzowe wykonanywane są zgodne z poprzednio
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 ∈ Rn×n jest dodatnio określona, gdy • A – jest symetryczna, tzn. A = AT,
• x 6= 0 =⇒ xT
Ax > 0, dla x ∈ Rn.
Definicja Dla macierzy A ∈ Rn×n
A = a11 a12 a13 . . . a1n a21 a22 a23 . . . a2n a31 a32 a33 . . . a3n .. . ... ... . .. ... an1 an2 an3 . . . ann macierz Ak∈ Rk×k, k = 1, . . . , n postaci Ak= a11 a12 a13 . . . a1k a21 a22 a23 . . . a2k a31 a32 a33 . . . a3k .. . ... ... . .. ... ak1 ak2 ak3 . . . akk
nazywamy k–tym minorem podstawowym macierzy A. Załóżmy, że dla macierzy A istnieją macierze trójkątne L, U
L = l11 0 0 . . . 0 l21 l22 0 . . . 0 l31 l32 l33 . . . 0 .. . ... ... . .. ... ln1 ln2 ln3 . . . lnn , U = u11 u12 u13 . . . u1n 0 u22 u23 . . . u2n 0 0 u33 . . . u3n .. . ... ... . .. ... 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.
Metoda rozkładu macierzy A na czynniki trójkątne (LU–rozkład): Z (5.12) mamy aij = n X s=1 lisusj= min(i,j) X s=1 lisusj. (5.13) Otrzymujemy stąd równości akk = k−1 X s=1 lksusk+ lkkukk (5.14) akj = k−1 X s=1 lksusj+ lkkukj, k + 1 ≤ j ≤ n (5.15) aik = k−1 X s=1 lisusk+ likukk, k + 1 ≤ i ≤ n (5.16)
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 | {z } znane +l22u22 Krok β : k = 2 a2j = l21u1j | {z } znane +l22u2j, j = 3, . . . , n ai2 = li1u12 | {z } znane +li2u22, i = 3, . . . , n Krok α : k = 3 a33 = 2 X s=1 l3sus3 | {z } znane +l33u33 Krok β : k = 3 a3j = 2 X s=1 l3susj | {z } znane +l33u3j, j = 4, . . . , n ai3 = 2 X s=1 lisus3 | {z } znane +li3u33, i = 4, . . . , n .. .
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;
• U = LT =⇒ l2
ii= aii, i = 1, . . . , n– algorytm Cholesky’ego.
W ten sposób otrzymujemy wzory:
lkkukk = k−1 X s=1 lksusk ukj = akj− k−1 P s=1 lksusj lkk , j = k + 1, . . . , n lik = aik− k−1 P s=1 lisusk ukk , i = k + 1, . . . , n
Twierdzenie Jeżeli wszystkie podstawowe minory macierzy A ∈ Rn×n są nie-osobliwe, to wtedy dla A istnieje LU– rozkład.
Dowód: Dowód zostanie przeprowadzony:::::::::::indukcyjnie względem wymiaru ma-cierzy n.
Dla n = 1 mamy a11 6= 0– jedyny minor macierzy A = [a11] ∈ R1×1 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 ∈ Rk×k. 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) 6= 0, det(Uk−1) 6= 0. Z równości (5.13) mamy: Ak−1 a1k .. . ak1 . . . akk = Lk−1 0 .. . lk1 . . . lkk Uk−1 u1k .. . 0 . . . ukk
Ak−1 a1k .. . ak1 . . . akk = Lk−1Uk−1 k−1 P s=1 l1susk .. . k−1 P s=1 lksus1 . . . k−1 P s=1 lksusk+ lkkukk (5.17)
Skąd otrzymujemy układy równań: l11u1k + l12u2k + . . . + l1k−1uk−1k = a1k l21u1k + l22u2k + . . . + l2k−1uk−1k = a2k l31u1k + l32u2k + . . . + l3k−1uk−1k = a3k .. . ... ... ... ... lk−11u1k + lk−12u2k + . . . + lk−1k−1uk−1k = ak−1k (5.18) oraz lk1u11 + lk2u21 + . . . + lkk−1uk−11 = ak1 lk1u12 + lk2u22 + . . . + lkk−1uk−12 = ak2 lk1u13 + lk2u23 + . . . + lkk−1uk−13 = ak3 .. . ... ... ... ... lk1u1k−1 + lk2u2k−1 + . . . + lkk−1uk−1k−1 = akk−1 (5.19)
W układzie (5.18) det Lk−16= 0 i u1k, u2k, u3k, . . . , uk−1k są niewiadomymi a w
(5.19) det Uk−16= 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
X
s=1
lksusk+ lkkukk= akk
gdzie można przyjąc, że lkk= 1.
Twierdzenie (rozkład Cholesky’ego) Jeżeli macierz A ∈ Rn×n 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 xTAx > 0 dla x = (x1, x2, . . . , xk, 0, . . . , 0)T 6= 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 = UTLT.
Skąd otrzymujemy:
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 = eLeLT, gdzie L = LDe 1/2.
Zadania na Ćwiczenia
1. Udowodnić, że:
(a) Jeżeli L–nieosobliwa jest dolnie trójkątna, to L−1też jest dolnie trój-kątna.
(b) Jeżeli L jest dolnie trójkątna i lii= 1, dla i = 1, . . . , n, to L−1istnieje
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
A =
0 1 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
A = 0 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
6
Wykład 7
Wykład 7
Twierdzenie 4rozkład Choleskiego5Jeżeli A jest dodatnio określoną macierzą, to istnieje L macierz dolnie trójkątna o własności:
A = LLT. Dowód Niech A = a11 a12 . . . a1n a21 a22 . . . a2n .. . ... . .. ... an1 an2 . . . ann , L = l11 0 . . . 0 l21 l22 . . . 0 .. . ... . .. ... ln1 ln2 . . . lnn .
Podzielimy macierze na bloki (patrz: wykład 6)
α = a11, a = a21 .. . an1 , A∗= a22 . . . a2n .. . . .. ... an2 . . . ann . oraz odpowiednio ρ = l11, r = l21 .. . ln1 , L∗= l22 . . . 0 .. . . .. ... ln2 . . . lnn .
Zgodnie z tymi oznaczeniami mamy, iż Aozn= α aT a A∗ , Lozn= ρ 0 r L∗ oraz LT = ρ rT 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 a A∗ = ρ 0 r L∗ ρ rT 0 LT ∗ = ρ2 ρrT rρ rrT + L ∗LT∗ . Z równości macierzy wynikają równości odpowiednich bloków
α = ρ2, aT = ρrT, A∗ = L∗LT∗ + rr
T,
4Dowód z książki G. W. Stewarta
5Inaczej liczył Tadeusz Banachiewicz (1882–1954), prof. UJ, który wprowadził rachunek
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 ρ 6= 0 i element l11 = ρ jest dobrze
określony. Z ρ 6= 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 ˆ A∗ ozn = A∗− rrT = A∗− 1 αaa T. Symetryczność macierzy ˆAT
∗ wynika z równości (oczywiście A∗= AT∗)
ˆ
AT∗ = (A∗− rrT)T = AT∗ − (rT)TrT = A∗− rrT = ˆA∗.
Niech y ∈ Rn×1, wtedy yTa = aTy i dla y 6= 0 udowodnimy
yTAˆ∗y = yTA∗y − 1 α(y Ta)(aTy) = yTA ∗y − 1 α(a Ty)2> 0. (6.23)
Macierz A jest dodatnio określona, więc z formy jej zapisu dla dowolnego η ∈ R wynika 0 < [ η yT ] α aT a A∗ η y = αη2+ 2aTyη + yTA∗y. (6.24) Do nierówności (6.24) podstawiamy η = −1 αa Ty 0 < αη2+ 2ηaTy + yTA∗y = yTA∗y − 1 α(a Ty)2
a to jest dokładnie warunek (6.23). W ten sposób udowodniliśmy „krok induk-cyjny” względem wymiaru macierzy.
2 W przypadku algorytmu Cholesky’ego wzory na LU–rozkład macierzy A (patrz wykład 6) redukują się do postaci:
lkk = v u u takk− k−1 X s=1 l2 ks lik = aik− k−1 P s=1 lislks lkk , i = k + 1, . . . , n
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 0 + a22−a21aa12 11 x2 + . . . + a2n−a21aa1n 11 xn = b2−aa21b1 11 . . . . 0 + an2−an1a11a12 x2 + . . . + ann−an1a11a1n xn = bn−an1a11b1
Zapowiada się, iź jeźeli 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)1 . . . a(k)1,k−1 a1k(k) . . . a(k)1j . . . a(k)1n . .. ... ... ... a(k)k−1,k−1 a(k)k−1,k . . . a(k)k−1,j . . . a(k)k−1,n 0 . . . 0 a(k)kk . . . a(k)kj . . . a(k)kn 0 . . . 0 a(k)k+1,k . . . a(k)k+1,j . . . a(k)k+1,n .. . ... ... ... ... 0 . . . 0 a(k)ik . . . a(k)ij . . . a(k)in .. . ... ... ... ... 0 . . . 0 a(k)nk . . . a(k)nj . . . a(k)nn (6.25)
otrzymać macierz A(k+1) należy zastosować wzory:
a(k+1)ij = a(k)ij , dla i ≤ k, a(k)ij −a (k) ik a(k)kka (k) kj, dla i ≥ k + 1, j ≥ k + 1, 0, dla i ≥ k + 1, j ≤ k, (6.26)
oraz dla kolumny wyrazów wolnych
b(k+1)i = b(k)i , dla i ≤ k, b(k)i −a (k) ik a(k)kkb (k) k , dla i ≥ k + 1 . (6.27)
Przedyskutujemy teraz dopuszczalność przedstawionej procedury, która jest warun-kowana przez a(k)kk 6= 0. Z postaci macierzy (6.25) wynika, iż
det(A(k)) = a(k)11 · . . . · a(k)k−1,k−1· a(k)kk . . . a(k)kn .. . ... a(k)nk . . . a(k)nn .
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) 6= 0 to dla każdego kroku k istnieje i ∈ {k, . . . , n} takie, że a(k)ik 6= 0. Więc w przypadku a(k)kk = 0 przy założeniu det(A) 6= 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) 6= 0, to wzory (6.26) i (6.27) z uwzględ-nieniem zamiany wierszy w przypadku, gdy w k-tym kroku a(k)kk = 0 sprowadza układ
Ax = b do równoważnego
U x = b(n) gdzie macierz U jest górnie trójkątna.
Przykład:6W przykładzie bedziemy tylko stosować wzory (6.26) i (6.27), bez zamiany wierszy w przypadku, gdy „w k-tym kroku a(k)kk = 0”.
Więc przy takich założeniach nie poradzimy sobie z układem 0 1 1 1 x1 x2 = 1 2
Poprawimy ten układ. Więc dla 0 < ε do układu: ε 1 1 1 x1 x2 = 1 2
możemy już procedurę zastosować. Otrzymujemy wtedy ε 1 0 1 − ε−1 x1 x2 = 1 2 − ε−1 .
Rozwiązując otrzymany układ równań („górnie trójkątny”) otrzymujemy, gdzie znakiem ≈ zaznaczyliśmy wpływ przybliżeń komputerowych na obliczenia:
( x2 = 2−ε −1 1−ε−1 ≈ 1 x1 = 1−xε−12 ≈ 0 6za D. Kincaid i W. Cheney
Sprowadzenie układu do postaci 1 ε−1 1 1 x1 x2 = ε−1 2
nie poprawia sytuacji 1 ε−1 0 1 − ε−1 x1 x2 = ε−1 2 − ε−1
otrzymujemy tak jak poprzednio
x2 = 2−ε
−1
1−ε−1 ≈ 1
x1 = ε−1− ε−1x2≈ 0 Dopiero zamina wierszy
1 1 ε 1 x1 x2 = 2 1 prowadzi do układu 1 1 0 1 − ε x1 x2 = 2 1 − 2ε
który, nie jest już tak „czuły” na przybliżenia komputerowe. x
2 = 1−2ε1−ε ≈ 1 x1 = 2 − x2≈ 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.
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 xn+1= xn− f (xn) g(xn) , gdzie g(x) = f (x + f (x)) − f (x) f (x) . (7.28)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 0(x + θf (x)), dla pewnego 0 < θ < 1. Wprowadzamy oznaczenie: ¯xn ozn. = xn + θnf (xn). Wtedy wzór (7.28) można
zapisać w następujący sposób:
xn+1= xn−
f (xn)
f0(¯x n)
.
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 ¯ en ozn. = ¯xn− r = (7.29) = xn− r + θn(f (xn) − f (r)) = en(1 + θnf0(ζn)), ζn∈ (r, xn).
Tak jak w dowodzie metody Newton’a wyliczamy: en+1= xn+1− r = xn− f (xn) f0(¯x n) − r = enf 0(¯x n) − f (xn) f0(¯x n) . (7.30)
Z dowodu metody Newton’a mamy: enf0(xn) − f (xn) =
1 2f
00(ξ
n)e2n, gdzie ξn ∈ (r, xn).
Ze wzoru Lagrange’a zastosowanego do pochodnej funkcji y = f0(x) otrzymu-jemy (z (7.29) i f (r) = 0 mamy f (xn) = f (xn) − f (r) = f0(ζn)en): f0(¯xn) = f0(xn+ θnf (xn)) = = f0(xn) + f00(ηn)θnf (xn) = f0(xn) + f00(ηn)θnf0(ζn)en, czyli enf0(¯xn) − f (xn) = 1 2f 00(ξ n)e2n+ f00(ηn)θnf0(ζn)e2n.
Wstawiając ostatnie wyrażenie do (7.30) otrzymujemy: en+1= 1 2f 00(ξ n) + f00(ηn)θnf0(ζn) f0(¯x n) | {z } ograniczone w otoczeniu r e2n. (7.31)
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ę
c(δ)def= max |x−r|≤δ|f 00(x)| 1 2+ max |x−r|≤δ|f 0(x)| min |x−r|≤δ|f 0(x)| .
Teraz, jeżeli y = f00(x) jest funkcją ciągłą i f0(r) 6= 0, to
lim
δ→0δc(δ) = 0.
Istnieje więc taka δ > 0, że ρozn.= δc(δ) < 1. Proces iteracji zaczniemy z punktu x0spełniającego: |x0− r| < δ oraz |x0− r| 1 + | max |x−r|≤δ|f 0(x)|< 1.
Z ostatniej nierówności oraz z (7.29) i oszacowania |xn− r| ≤ ρn|x0− r| wynika,
że jeżeli lim
n→∞xn= r to limn→∞x¯n= r co pozwala kontrolować zachowanie
miano-wnika w (7.31). W ten sposób znaleźliśmy się w warunkach dowodu zbieżności ciągu przybliżeń en z metody Newton’a.
2 . . . .
B.
Przybliżenie y = f0(x) – z listy do wykładu z dnia ??.XI.2000r Korzystając ze wzoru Taylor’a uzasadnićf0(x) ≈ k
2f (x + h) − h2f (x + k) + (h2− k2)f (x)
(k − h)kh Rozwiązanie
Będziemy korzystali z wzoru Taylor’a z resztą w postaci całki7
f (x) = n X i=0 f(i)(x 0) i! (x − x0) i+ 1 n! x Z x0 f(n+1)(t)(x − t)ndt | {z } rn(x) .
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 Z a f (x)g(x) dx = f (ξ) b Z a g(x) dx.
Korzystamy ze wzoru Taylor’a dla n = 1 aby otrzymać: k2f (x + h) − h2f (x + k) + (h2− k2)f (x) (k − h)hk = = k 2[f (x) + f0(x)h + r 1(x + h)] − h2[f (x) + f0(x)k + r1(x + k)] (k − h)hk + +(h 2− k2)f (x) (k − h)hk = = f 0(x)(hk2− h2k) + k2r 1(x + h) − h2r1(x + k) (k − h)hk = = f0(x) +k 2r 1(x + h) − h2r1(x + k) (k − h)hk | {z } ρ(h,k) .
Pokażemy, że dla ρ = ρ(h, k) ρ(h, k)ozn.= k 2r 1(x + h) − h2r1(x + k) (k − h)hk (7.32) zachodzi lim (h,k)→(0,0) ρ(h, k) = 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) = k 2 2 x+h Z x f00(t)(x + h − t) dt −h 2 2 x+k Z x f00(t)(x + h − t) dt = = k 2 2 x+h Z x+k f00(t)(x + h − t) dt | {z } I1 +1 2 x+k Z x f00(t)k2(x + h − t) − h2(x + k − t) dt. | {z } I2
Każdą z zaznaczonych całek I1, I2 będziemy szacowali oddzielnie. Pierwszą
z nich szacujemy następująco
I1= k2 2 x+h Z x+k f00(t)(x + h − t) dt = f 00(ξ 1)k2 2 x+h Z x+k (x + h − t) dt = = f 00(ξ 1)k2 2 −(x + h − t) 2 2 x+h x+k = k2(h − k)2f 00(ξ 1) 4
gdzie ξ1∈ [x + k, x + h]. Drugą szacujemy podobnie I2= 1 2 x+k Z x f00(t)k2(x + h − t) − h2(x + k − t) dt = = (k 2− h2) 2 x+k Z x f (t)(x − t) dt +(k − h)hk 2 x+k Z x f (t) dt = = (k 2− h2)f00(ξ 2) 2 −(x − t) 2 2 x+k x +(k − h)hk 2 2 f 00(ξ 3) = = −(k2− h2)k2f00(ξ2) 4 + (k − h)hk 2f00(ξ3) 2 gdzie ξ2, ξ3∈ [x, x + k]. Jeżeli założymy, że druga pochodna y = f00(x) jest ograniczona, to otrzymamy następujące oszacowania dla całek I1, I2
|I1| ≤ k2(h − k)2C i |I2| ≤ (|(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ż: |ρ(h, k)| ≤ |I1| |(k − h)hk|+ |I2| |(k − h)hk| ≤ ≤ |k 2(h − k)2| |(k − h)hk| + |(k2− h2)k2| + |(k − h)hk2| |(k − h)hk| C ≤ ≤ (|h − k| + |k + h| + |k|)C ≤ 3(|h| + |k|)C ≤ 3√2ph2+ k2C.
Z ostatniej nierówności wynika lim (h,k)→(0,0) ρ(h, k) = 0. 2 . . . . 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
gdzie8X–em są zaznaczone elementy macierzy niekoniecznie równe zeru. W ko-lejnym k–tym kroku wymagane jest aby akk 6= 0, gdyby ten warunek nie był
spełniony to wystarczy wśród akk, . . . , ankznaleźć 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
k≤i≤n|aik|
(czę-ś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
k≤i≤nk≤j≤nmax |aij| (pełny wybór elementu
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 X X X . . . X X X 0 X X . . . X X X 0 0 X . . . X X X .. . ... ... . .. . .. ... ... 0 0 0 . . . 0 X X , d1 c1 a1 d2 c2 a2 d3 c3 . .. . .. . .. a−2 dn−1 cn−1 an−1 dn .
Analiza błędów
9 10 11 12 13 7.1 Definicja Dla macierzy A ∈ Rn×n A = a11 a12 . . . a1n a21 a22 . . . a2n .. . ... . .. ... an1 an2 . . . ann funkcję k · k : Rn×n−→ R+ nazywamy normą macierzową, gdy ma następujące
własności
1. kAk ≥ 0, oraz kAk = 0 ⇐⇒ A = 0, 2. kcAk = |c|kAk, c ∈ R, 3. kA + Bk ≤ kAk + kBk, 4. kABk ≤ kAk kBk. Typowe normy kAk2 = s n P i=1 n P j=1 a2
ij (norma euklidesowa, `2–norma)
kAkS = max
i pλi(AA
T) (norma spektralna)
9D. Kincaid & W. Cheney, Numerical Analysis, ss. 181–221; 10G.W. Stewart, Afternotes on Numerical Analysis, ss. 103–131; 11A. Ralston, Wstęp do analizy numerycznej, ss. 365–392;
12R.L. Burden, J.D. Faires & A.C. Reynolds, Numerical Analysis, ss. 289–317; 13H. Rutishauser, Lectures on Numerical Mathematics, ss. 23–76;
Gdzie w definicji normy14spektralnej kAk
S symbolem λi(AAT) oznaczono i–tą
wartość własną macierzy AAT.
Uwaga Z macierzą A ∈ Rn×n jest związane wzorem f (x) = Ax
przekształcenie liniowe f : Rn−→ Rn. Przy tak przyjętych oznaczeniach mamy
kAk = kf kdef= sup
kxk=1
kf (x)k = sup
x6=0
kf (x)k kxk .
Inaczej, norma macierzy kAk 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ń kf kdef= sup x6=0 kf (x)k kxk jest normą. Normy w przestrzeniach Rn: kxk1 = n P i=1 |xi| `1–norma, kxk2 = s n P i=1 x2
i `2–norma, norma euklidesowa,
kxk∞ = max
1≤i≤n|xi| `∞–norma,
gdzie xT = (x1, x2, . . . , xn) ∈ Rn.
7.2 Lemat
Jeżeli przestrzeń Rn jest z normą kxk = max
1≤i≤|xi|, x ∈ R n to dla A ∈ Rn×n zachodzi kAk∞= max 1≤i≤n n X j=1 |aij|. Uzasadnienie kAk∞ = sup kuk∞=1 kAuk∞= sup kuk∞=1 max 1≤i≤n|(Au)i| = = max 1≤i≤n ( sup kuk∞=1 |(Au)i| ) = max 1≤i≤n ( sup kuk∞=1 |(Au)i| ) = = max 1≤i≤n n P j=1 |aij|
gdzie kres górny sup
kuk∞=1 n P j=1 aijuj = n P j=1 |aij|, i = 1, . . . , n jest realizowany
przez kuk∞= 1 określone w następujący sposób uj = sgn aij, i = 1, . . . , n.
2 Przykład I
Niech x ∈ Rn będzie rozwiązaniem dokładnym równania
Ax = b a ˜x ∈ Rn równania zaburzonego
B ˜x = 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 kx − ˜xk = kx − Bxk = kx − BA−1xk = k(I − BA−1)xk ≤ kI − BA−1k kxk Co daje następujące oszacowanie błędu względnego
kx − ˜xk
kxk ≤ kI − BA
−1k.
2 Przykład II
Macierze układów są dokładne, natomiast prawe strony są zaburzone, tzn. Ax = b – równanie „dokładne”,
A˜x = ˜b – równanie z zaburzoną prawą stroną.
Pytanie: (tak jak poprzednio) Jaki jest błąd względny otrzymanego wyniku ? Mamy x = A−1b – rozwiązanie dokładne, ˜ x = A−1˜b – rozwiązanie przybliżone. Teraz jest kx − ˜xk ≤ kA−1b − A−1˜bk = kA−1(b − ˜b)k ≤ kA−1k kb − ˜bk. Skąd wynika (wiemy, że Ax = b)
kb − ˜bk ≤ kA−1kkBxk
kbk kb − ˜bk ≤ kA
−1k kAk kxk kb − ˜bk
kbk Otrzymaliśmy w ten sposób następujące oszacowanie błędu względnego
kx − ˜xk kxk ≤ kA −1k kAk | {z } κ(A) kb − ˜bk kbk
Liczba κ(A) („condition number”) związana z macierzą nieosobliwą A ∈ Rn×n
κ(A)def= kA−1k kAk
odgrywa zasadniczą rolę w oszacowaniu błędu względnego rozwiązań układu rów-nań przez błąd względny prawych stron
kx − ˜xk
kxk ≤ κ(A) kb − ˜bk
kbk .
2 Przykład III (liczbowy)
Dla ε > 0 niech A = 1 1 + ε 1 − ε 1 .
Mamy det(A) = ε2, więc macierz A jest nieosobliwa. Wtedy macierz odwrotna
A−1 jest określona wzorem
A−1= 1 ε2 1 −1 − ε −1 + ε 1 . Dla normy macierzowej k · k∞mamy
kAk∞= 2 + ε i kA−1k∞= 2 + ε ε2 skąd otrzymujemy, iż κ(A) = 2 + ε ε2 2 > 4 ε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ń.
2 . . . .
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˜a ≈ b,
to wtedy błąd względny przybliżonego rozwiązania ˜x spełnia
1 κ(A) kb − A˜xk kbk ≤ kx − ˜xk kxk ≤ κ(A) kb − A˜xk kbk , (7.34) gdzie
κ(A)def= kA−1k kAk. Dowód Dowodzimy prawej nierówności w (7.34). Mamy
kx − ˜xk kbk = kA−1(b − A˜x)k kAxk ≤ kA−1k kb − A˜xk kAk kxk skąd otrzymujemy kx − ˜xk kxk ≤ kA −1k kAk | {z } κ(A) kb − A˜xk kbk . Teraz lewa strona (7.34). Zachodzi
kb − A˜xk kxk = kA(x − ˜x)k kA−1bk ≤ kAk kx − ˜xk kA−1k kbk co daje kb − A˜xk kbk ≤ kA −1k kAk | {z } κ(A) kx − ˜xk kxk . 2 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 źle 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.
Zadania na Ćwiczenia
1. Pokazać, że funkcje kxk∞, kxk2, kxk1 są normami w Rn.
2. Wykazać nierówności
kxk∞≤ kxk2≤ kxk1
dla x ∈ Rn.
3. Wykazać, że jeżeli f, g : Rn −→ R są ograniczone, to wtedy zachodzi sup[f (x) + g(x)] ≤ sup f (x) + sup g(x).
4. Jeżeli w Rn jest określona norma wektorowa kxk = Pn
i=1
|xi|, to
odpowia-dająca norma macierzowa jest określona wzorem kAk∞= max 1≤j≤n n X i=1 |aij|.
5. Dla normy macierzowej k · k∞ obliczyć κ(A) dla macierzy
A = 7 8 9 10 .
6. Niech n = 3 i dla macierzy
A = 4 −3 2 −1 0 5 2 6 −2 obliczyć sup kxk∞≤1 kAxk∞.
7. Obliczyć κ(A) dla norm k · k1, k · k2, k · k∞ i macierzy
a) a + 1 a a a − 1 , b) 0 1 −2 0 , c) α 1 1 1 . 8. Dla p ≥ 1 i x ∈ Rn określamy kxk p def = (Pn i=1|xi|) 1 p. Wykazać, że k · k ∞
jest normą dla której zachodzi limp→∞kxkp= kxk∞.
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 źle uwarunkowanych macierzy wraz z otrzy-manymi wyznaczonymi przez nie numerycznymi rozwiązaniami układów równań.
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 mnoźenia 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, iź mnożenie może być wykonywane szybciej niź do-dawanie!? Zaleźy 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)
2 ≈
n2 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) = kAk∞kA−1k∞
występującym w oszacowaniu błędów wzgędnych. Niech A = 1 1 1 2 wtedy A−1= 2 −1 −1 1 oraz κ(A) = 3 · 3 = 9. Mnożymy teraz pierwszy wiersz układu Ax = b przez liczbę 10−4 aby otrzymać:
ˆ A = 10−4 10−4 1 2 wtedy A−1= 2 · 104 −1 −104 1 oraz κ(A) = 3 · (2 · 104+ 1) ≈ 6 · 104.
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 κ(A) kb − A˜xk kbk ≤ kx − ˜xk kxk ≤ κ(A) kb − A˜xk kbk .
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 ≤ i ≤ n.
8.1 Twierdzenie
Jeżeli x0, x1, . . . , xn są liczbami spełniającymi warunek xi 6= xj dla i 6= 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 ≤ i ≤ 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 ≤ i ≤ n.
Ponieważ wielomian stopnia nie większego od n ma n + 1 pierwiastków, więc zachodzi
pn(x) − qn(x) ≡ 0.
Istnienie (dowód indukcyjny). Wielomiany pk dla k = 0, . . . , n są wyznaczane
ze wzoru pk(x) = y0, dla k = 0, pk−1(x) + ck(x − x0)(x − x1) . . . (x − xk−1), dla k = 1, . . . , n. (8.35) gdzie ck= yk− pk−1(xk) (xk− x0)(xk− x1) . . . (xk− xk−1) . (8.36) 2 Wielomiany interpolacyjne pk(x), k = 1, . . . , n w (8.35) nazywają sie wielomianami
interpolacyjnymi Newtona. Samo obliczanie warości wielomianów może odby-wać sie według schematu Hornera
gdzie
dk = x − xk.
Wielomiany interpolacyjne Lagrange’a Metoda wynika z następującego przedstawienia
pn(x) = y0(x) + y1`1(x) + . . . + yn`n(x) = n
X
k=0
yk`k(x).
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
δij = pn(xj) = n X k=0 yk`k(xk) = n X k=0 δkj`k(xj) = `i(xj).
Jeżeli dodatkowo założymy, że `i(x) są wielomianami stopnia n o własności
`i(xj) = 0 dla i 6= 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 − x\i) oznacza, iż ten czynnik został opuszczony. Ostatecznie
otrzymu-jemy, iż `i(x) = n Y j=0 j6=i x − xj xi− xj . 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 x20 . . . xn0 1 x1 x21 . . . xn1 .. . ... ... . .. ... 1 xn x2n . . . xnn a0 a1 .. . an = y0 y1 .. . yn .
Otrzymaliśmy w ten sposób układ równań typu V a = y na współczynniki aT = [a
0, a1, . . . , an] gdzie macierz V jest macierzą Vandermonde’a. Ten sposób
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
f (x) − p(x) = 1 (n + 1)!f (n+1)(ξ x) n Y i=0 (x − xi). (8.37)
Dowód Ustalamy x ∈ [a, b], x 6= xii = 0, . . . n. Niech
w(t) =
n
Y
i=0
(t − xi), ϕ(t) = f (t) − p(t) − λw(t), dla t ∈ [a, b], (8.38)
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, x
0, 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 ϕ0 = ϕ0(t), ϕ0 ∈ 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 0 = ϕ(n+1)(ξx) = f(n+1)(ξx) − (n + 1)! λ = f(n+1)(ξx) − (n + 1)! f (x) − p(x) w(x) .
Co po uwzględnieniu określenia w = w(t) z (8.38) daje (8.37).
9
Wykład 11
Wykład 11 Wielomiany iterpolacyjne Czebyszewa16
Wielomiany Czebyszewa pierwszego radzaju określamy rekurencyjnie w nastę-pujący sposób: Tn(x) = 1, dla n = 0; x, dla n = 1; 2x Tn(x) − Tn−1(x), dla 1 < n.
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) = 2x fn(x) − fn−1(x) dla 1 ≤ n. 161821–1894
2 Ze wzoru (9.39) wynikają następujące własności wielomianów Czebyszewa:
|Tn(x)| ≤ 1 dla x ∈ [−1, 1] Tn cosnkπ = (−1)k, k = 0, 1, . . . , n Tn cos2k−12n π = 0, k = 1, 2, . . . , n. (9.41) 9.2 Twierdzenie
Dla wielomianów Czebyszewa jest spełniony następujący warunek ortogonalno-ści: 1 R −1 Tn(x) Tm(x)√1−xdx 2 = 0, dla n 6= m, π, n = m = 0, π 2, n = m 6= 0.
Dowód W dowodzie ortogonalności skorzystamy z postaci (9.39) wielomianów oraz podstawienia: x(ϕ) = cos ϕ, dx(ϕ) dϕ = − sin ϕ. Zachodzi więc 1 Z −1 Tn(x) Tm(x) dx √ 1 − x2 = = 1 R −1
cos(n arccos x) cos(m arccos x) √dx 1−x2 =
=
0
R
π
cos(nϕ) cos(mϕ)√− sin ϕ dϕ
1−cos2ϕ = π R 0 cos nϕ cos mϕ dϕ = = cos(n+m)ϕ+cos(n−m)ϕ2 π 0 = 0, dla n 6= m, π, n = m = 0, π 2, n = m 6= 0. 2 Z twierdzenia wynika, iż fukcję y = f (x) można przedstawić w postaci szeregu Czebyszewa: f (x) ∼c0 2 + ∞ X n=1 cnTn(x)
gdzie współczynniki (cn)∞n=0 są określone wzorem:
cn = 1 Z −1 f (x) Tn(x) dx √ 1 − x2.
Podobnie jak w przypadku szeregu Fouriera należy wykazać zbieżność otrzyma-nego szeregu do zadanej funkcji y = f (x).
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 ≤ 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 kpnk∞ def = max −1≤x≤1|pn(x)| ≥ 2 1−n.
Dowód (nie wprost) Załóżmy przeciwnie, iż
|pn(x)| < 21−n, dla x ∈ [−1, 1].
Niech
qn(x) = 21−nTn(x) oraz xk= cos
kπ
n dla k = 0, . . . , n. Z przyjętych założeń oraz z (9.41) wynika, iż
(−1)kpn(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) w
kolej-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.
2 Niech x0, x1, . . . , xn∈ [−1, 1] spełniają x0< x1< . . . < xn. Wtedy na
podsta-wie tpodsta-wierdzenia z poprzedniego wykładu mamy dla podsta-wielomianu p = p(x) stopnia co najwyżej n następujące oszacowanie
max −1≤x≤1|f (x) − p(x)| ≤ 1 (n + 1)!−1≤x≤1max |f (n)(x)| max −1≤x≤1 n Y k=0 (x − xk) (9.42)
Z twierdzenia 9.3 mamy następujące oszacowanie (dla wielomianu n + 1 stop-nia) max −1≤x≤1 n Y k=0 (x − xk) ≥ 1 2n. (9.43)
Ponieważ wielomian o współczynniku przy najwyższej potędze równym 1 jest jednoznacznie wyznaczony przez swoje pierwiastki, to z (9.41) dla
xk= cos
2k + 1 2n + 2π
