Numerické metody. Vratislava Mošová

Numerické metody Vratislava Mošová

1

Předmluva S rozvojem počítačů vzrostl význam numerických metod. Řada výpočetních postupů, které sem řadíme, vznikla sice už dávno předtím, ale teprve nástup výpočetní techniky vytvořil předpoklady pro vznik a rozvoj nových metod, které ke své realizaci potřebují velký objem výpočtů. Vedle klasických metod jako jsou Gaussova a Jakobiova metoda pro soustavy lineárních rovnic, Newtonova metoda pro nelineární rovnice, Lagrangeova interpolace nebo obdélníková metoda výpočtu integrálu, zde teď stojí multigridní metody pro řešení soustav lineárních rovnic, aproximace křivek Bézierovými polynomy nebo metoda konečných prvků pro řešení diferenciálních rovnic. Algoritmů k numerickému řešení nejrůznějších matematických problémů je mnoho. Ve skriptech, která se vám dostávají do rukou, se setkáte s těmi nejzákladnějšími. K pochopení textu je nutná znalost základních pojmů z lineární algebry a matematické analýzy. Samotný text je rozdělen do šesti kapitol. V úvodní kapitole se budeme zabývat problematikou vzniku a šíření chyb, ke kterým dochází při zpracování numerické úlohy na počítači. Také zde připomeneme některá fakta o maticích a uvedeme některé pojmy a tvrzení z funkcionální analýzy, které budou užitečné v následujícím výkladu. Další tři kapitoly budou věnovány numerickým metodám řešení úloh algebry, jako je řešení soustav lineárních rovnic, určování vlastních čísel matice a výpočet kořenů nelineárních rovnic. Poslední dvě kapitoly budou zasvěceny numerickým metodám řešení úloh matematické analýzy, zejména problematice aproximace funkcí a otázkám numerické kvadratury. Na konci skript je seznam doporučené literatury. K připomenutí si potřebných pojmů a vět z matematické analýzy může posloužit [4]. Tituly [1], [3], [7], [5], [6], [8] jsou základní literaturou, která se týká studované problematiky. Ucelený přehled současných numerických metod poskytují knihy [9] a [11]. Text [2] obsahuje příklady k procvičení probíraných metod a ve skrip3

tech [10] najde čtenář další rozšíření studované problematiky včetně odkazů na numerický software. Cílem skript zdaleka není podat vyčerpávající obraz o metodách numerického řešení úloh. Jde pouze o první seznámení s těmito metodami. Osvojení si efektivních postupů by však mělo jít ruku v ruce se snahou proniknout do podstaty studované problematiky. Jedině ten, kdo si dokáže uvědomit úskalí, se kterými se může v průběhu výpočtu setkat, dokáže také efektivně využívat softwarové produkty, které jsou v současné době k dispozici.

4

Obsah Předmluva

3

1 Úvodní pojmy 7 Numerické metody a matematické modelování . . . . . . . . . . . . 7 O nepřesnostech při řešení problému . . . . . . . . . . . . . . . . . 9 Některá fakta z lineární algebry a funkcionální analýzy . . . . . . . 14 2 Řešení soustav lineárních rovnic Přímé metody . . . . . . . . . . . . . . . . . Gaussova eliminační metoda . . . . . . Metoda LU-rozkladu . . . . . . . . . . Iterační metody . . . . . . . . . . . . . . . . Jacobiova iterační metoda . . . . . . . Gaussova-Seidelova iterační metoda . . Podmíněnost soustav lineárních rovnic

. . . . . . .

. . . . . . .

. . . . . . .

3 Vlastní čísla a vlastní vektory matice Částečný problém vlastních čísel . . . . . . . . . . Mocninná metoda . . . . . . . . . . . . . . . Metoda Rayleighova podílu . . . . . . . . . Úplný problém vlastních čísel . . . . . . . . . . . Přímý výpočet vlastních čísel . . . . . . . . Určení vlastních čísel metodou LU-rozkladu Určení vlastních čísel metodou ortogonálních

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . transformací

. . . . . . . . . . . . . .

. . . . . . .

27 28 28 33 41 42 46 55

. . . . . . .

59 61 61 63 65 65 67 69

4 Řešení nelineárních rovnic 73 Řešení nelineární rovnice f (x) = 0 . . . . . . . . . . . . . . . . . . . 73 Metoda bisekce . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5

Metoda prosté iterace . Metoda regula falsi . . Newtonova metoda . . Metoda sečen . . . . . Řešení rovnic Pn (x) = 0 . . Bernoulliova metoda . Graefova metoda . . . Laguerrova metoda . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

5 Aproximace funkcí Interpolace . . . . . . . . . . . . . . . . . . . Lagrangeův interpolační polynom . . . . Newtonův interpolační polynom . . . . . Extrapolace . . . . . . . . . . . . . . . . Splajny . . . . . . . . . . . . . . . . . . Aproximace trigonometrickými polynomy Hermitova interpolace . . . . . . . . . . Bézierovy křivky . . . . . . . . . . . . . . . . Metoda nejmenších čtverců . . . . . . . . . . 6 Numerická kvadratura a derivace Newtonovy-Cotesovy kvadraturní vzorce Složené kvadraturní vzorce . . . . . Gaussovy kvadraturní vzorce . . . . . . Numerická derivace . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

76 80 82 85 86 89 90 93

. . . . . . . . .

. . . . . . . . .

95 96 97 101 106 108 115 119 120 124

. . . .

129 . 130 . 133 . 139 . 142

Literatura

144

Rejstřík

145

6

Kapitola 1 Úvodní pojmy Numerické metody a matematické modelování V životě potřebuje člověk řešit velké množství praktických problém˚ u. Aby se s nimi vypořádal, vytvořil si řadu nástroj˚ u, které mu to umožní. Matematika je jedním z nich. Tak například úlohu týkající se bezpečné konstrukce mostu lze modelovat pomocí diferenciálních rovnic. Jinou úlohu m˚ uže představovat problém, jak optimálně řídit stavbu nějakého areálu, což se dá s určitými zjednodušeními zachytit prostřednictvím soustavy lineárních rovnic. Řešit takové úlohy přímo však není možné. To je jeden z d˚ uvod˚ u, proč je zapotřebí rozvíjet numerickou matematiku. Dalším důvodem je fakt, že algebra a matematická analýza sice umožňují rozhodnout o existenci a jednoznačnosti řešení nějakého problému (viz základní věta algebry), ale pouze v některých případech dávají konkrétní návod, jak tyto problémy řešit. (Např. v důkazu věty o existenci a jednoznačnosti řešení Cauchyova počátečního problému se konstruuje posloupnost postupných aproximací řešení.) Někdy zase je způsob řešení uveden, ale z hlediska realizace výpočtu nemusí být zrovna optimální. (Např. když budeme na počítači řešit soustavy lineárních rovnic Cramerovým pravidlem, bude výpočet hodnot determinantů vyšších řádů působit problémy.) Cílem numerických metod je proto vytvořit efektivní algoritmy pro řešení nejrůznějších matematických problémů. Formulace úloh a způsob jejich řešení je závislý na skutečnosti, že pracujeme s počítačem. Práce na počítači vyžaduje, abychom zadali na vstup konečný počet číselných údaj˚ u a postup, 7

prostřednictvím kterého po konečném počtu krok˚ u dojdeme k výsledku. Ten je dán opět konečným počtem výstupních číselných údaj˚ u. Jinými slovy, naším cílem je pomocí vhodného algoritmu vyřešit numerickou úlohu. Definice 1.1 Numerická úloha je jednoznačný funkční popis vztahu mezi konečným počtem vstupních a konečným počtem výstupních dat. Jednotlivé matematické úlohy však ještě nemusejí být úlohami numerickými. Příklad Řešit soustavu lineárních rovnic, najít kořeny polynomu, určit vlastní čísla matice jsou numerické úlohy. Avšak integrování a derivování funkce nebo řešení soustavy diferenciálních rovnic nejsou numerické úlohy, protože je nelze zadat konečným počtem vstupních dat. Definice 1.2 Algoritmus je jednoznačně zadaná konečná posloupnost operací, která m-tici přípustných vstupních dat přiřadí n-tici dat výstupních. Neřeší jen jedinou speciální úlohu, ale poskytuje řešení pro celou skupinu úloh téhož typu. √ Příklad Napište algoritmus pro výpočet odmocniny a, a > 0. √ Řešení: Když označíme x = a, pak x2 = a, neboli 2x2 − x2 = a. A tedy pro x 6= 0 máme 2x = x + xa a x = 12 (x + xa ). Odtud už obdržíme iterační formuli, tzv. Heron˚ uv vzorec, 1 a xn = (xn−1 + ), n = 1, 2, . . . 2 xn−1 Algoritmus (Heron˚ uv vzorec) Vstup: a, x0 , n pro k = 1, 2, . . . , n a xk = 21 (xk−1 + xk−1 ) Výstup: xn - přibližná hodnota odmocniny

8

√

a

O nepřesnostech při řešení problému Z posledního příkladu se vidí, že hodnota odmocniny, kterou získáme pomocí Heronova vzorce, nedává přesný výsledek. Tato vlastnost je charakteristická pro řadu numerických výpočtů. Některé numerické metody sice umožňují dojít po konečném počtu kroků k přesnému řešení (např. Gaussova eliminace), ale velká část numerických řešení má pouze aproximativní charakter (např. obdélníkové pravidlo, metoda bisekce). Z toho vyplývá nutnost zabývat se takovými otázkami, jako je odhad chyby, podmíněnost úlohy či stabilita algoritmu. Navíc v souvislosti s matematickým modelováním dochází při vytváření matematického modelu i při jeho řešení k celé řadě dalších nepřesností. Zdrojem chyb mohou být 1. chyby vzniklé při vytváření modelu, protože — naměřená data díky přístroji nemusejí být přesná, — matematický model je zjednodušením už zjednodušeného fyzikálního modelu a neodpovídá přesně realitě. 2. chyby vzniklé při zpracování úlohy na počítači, protože — čísla se zobrazují v počítači v semilogaritmickém tvaru a celá mantisa se nemusí vždy do počítače vejít. Navíc množina počítačových čísel není vzhledem k algebraickým operacím uzavřená, což m˚ uže zp˚ usobit další kumulování chyb, — úloha je špatně podmíněná (viz dále strana 12), — algoritmus je nestabilní (viz dále strana 13).

Druhému typu chyb se nyní budeme věnovat podrobněji. Nejprve si všimněme, jak je to s chybami, které vznikají při zobrazení čísla v počítači. Číslo x zapsané v semilogaritmickém tvaru x = sgn x (

x1 x2 + 2 + . . . ) q b , x1 6= 0, q q 9

se v q-adickém t-místném počítači M (q, t) uchovává jako x˜ = sgn x (

x1 x2 xt + 2 + · · · + t ) qb, q q q

kde q > 1 je základ, b je exponent a xi ∈ {1, 2, . . . , q − 1}, i = 1, . . . t, jsou číslice mantisy. Pokud se číslo nem˚ uže celé v počítači zobrazit, dochází k jeho řezání nebo zaokrouhlení, a tím i ke vzniku chyby zobrazení. Přitom pro absolutní chybu ∆x platí b

∆x = |x − x˜| = q |

∞ X xi i=1

qi

−

½

t X xi i=1

q

q b q −t 1 b −t qq 2

|≤ i

při řezání, při zaokrouhlování.

x Pro relativní chybu | x−˜ | na vstupu je x

qb| x − x˜ | |= x neboť |

P∞

xi i=1 q i |

P∞

≥|

Pt xi xi i=1 q i − i=1 q i | P ∞ q b | i=1 xqii |

P∞

1 i=1 q i |

½ ≤

q 1−t 1 1−t q 2

při řezání, při zaokrouhlování,

= 1q .

Díky zaokrouhlovacím chybám vznikají také chyby ve výsledcích aritmetických operací. Příklad. V M (10, 3) sečtěte 1,24 a 0,0221. Řešení: Sčítání probíhá ve sčítacím registru, který má navíc jedno rezervní místo. Nejprve se provede porovnání exponent˚ u s případnou denormalizací. Pak se sečtou mantisy a výsledek se vyjádří v normálním tvaru. . 0, 124 · 101 + 0, 221 · 10−1 = (0, 124|0 + 0, 002|21) · 101 = 0, 126 · 101 .

Je také třeba si uvědomit, že pro aritmetické operace prováděné počítačem obecně neplatí ani asociativní ani komutativní zákon. Příklad. V M (10, 3) pomocí částečných součtů odhadněte součet řady 10

P∞ n=0

0, 9n .

Řešení: Protože se jedná o součet geometrické řady s kvocientem q = 0, 9, 1 = 10. Částečné součty spočteme dvěma je přesná hodnota součtu s = 1−0,9 způsoby. Označme sk = ((1 + 0, 9) + 0, 92 ) + · · · + 0, 9k ) — sčítání odpředu, rk = ((0, 9k + 0, 9k−1 ) + · · · + 1) — sčítání odzadu. Výsledky, které získáme pro různá k, zapišme do tabulky k 50 100 150

sk rk s − sk 9.98 9.97 0.02 9.98 10.0 0.02 9.98 10.0 0.02

s − rk 0.03 0 0

V prvním případě, kdy sčítáme od největších členů, se nejdřív rychle blížíme k součtu geometrické řady, ale od jistého k počínaje je příspěvek dalších členů řady ve srovnání se stávajícím součtem příliš malý a nezapočítává se. Oproti tomu přičítání odzadu umožňuje do celkového součtu zahrnout příspěvek všech sčítanců. Definice 1.3 Řekneme, že v normalizovaném tvaru aproximace t X xi x˜ = sgn x ( ) · 10b i 10 i=1

čísla

∞ X xi x = sgn x ( ) · 10b i 10 i=1

je s číslic platných, když vztah |x − x˜| ≤ 0, 5 · 10b−j platí jedině pro j ≤ s. Příklad Určete počet platných číslic dekadického rozvoje Eulerova čísla, když x˜ = 2, 718. 11

Řešení: x = e = 2, 718218 · · · = 10 · 0, 2718218 . . . , |x − x˜| ≤ 10−3 · 0, 218 · · · ≤ 0, 5 · 101−4 , V dané aproximaci Eulerova čísla jsou 4 číslice platné. Přirozené důvody nás vedou k požadavku, aby úloha, kterou se zabýváme, byla jednoznačně řešitelná a její řešení aby spojitě záviselo na vstupních datech. Jinými slovy, aby úloha byla korektní a dobře podmíněná. Definice 1.4 Nechť X a Y jsou Banachovy prostory 1 , X je množina vstupních dat a Y množina výstupních dat. Řekneme, že úloha U je korektní, když ∀x ∈ X ∃! y ∈ Y : y = U x, xn → x, yn = U (xn ) ⇒ U (xn ) → U (x) = y.

Poznámka Můžeme říci, že numerická úloha je korektní, když pro každý vstup existuje právě jedno její řešení, které spojitě závisí na vstupních datech. Definice 1.5 Úloha je dobře podmíněná, když je korektní a číslo podmíněnosti relativní chyba na výstupu cp = relativní chyba na vstupu má hodnotu blízkou 1. Poznámka Pokud cp ≈ 1, znamená to, že malá změna na vstupu vyvolá pouze malou změnu na výstupu. Příklad Určete číslo podmíněnosti pro úlohu nalézt hodnotu polynomu p(x) = x2 + x − 1150 v bodě x = 33. Volte x ≈ x˜ = 1

100 . 3

Banachův prostor viz definice 1.16.

12

Řešení: Protože −28− 50

| −28 9 | 50 100 p(33) = −28, p( ) = − , cp = 33− = 50 3 9 | 33 9 |

22,4 28 1/3 33

. = 79,

je úloha špatně podmíněná. Úspěch numerického řešení problému také závisí na volbě algoritmu. Pokud malá nepřesnost ve vstupních datech začne generovat velké chyby během výpočtu, nepovede zvolený postup ke správnému výsledku. Zda určitý algoritmus je vhodný pro daný problém či ne, souvisí s otázkou stability algoritmu. Algoritmus je stabilní, když je málo citlivý na poruchy v datech a na vliv zaokrouhlovacích chyb během výpočtu. Stabilitu algoritmů lze studovat pomocí numerických experimentů. Příklad Je dána rovnice x2 −2bx+c = 0. Její kořeny lze určit prostřednictvím algoritmu √ √ A1: x1 = b + b2 − c, x2 = b − b2 − c nebo algoritmu √ √ √ 2 −c)(b+ b2 −c) √ A2: x1 = b + b2 − c, x2 = (b− bb+ = b+√cb2 −c . (viz [3]) b2 −c Ukážeme, že A1 je citlivější na vliv zaokrouhlovacích chyb během výpočtu, a tedy je méně stabilní než A2. K prověření přesnosti řešení využijeme vztahu, který platí mezi koeficienty a kořeny kvadratické rovnice: c = x1 x2 . Řešení: Pro rovnici x2 − 105 x + 1 = 0 je v M (10, 8) : x1 = 100000, 00. Když spočteme druhý kořen pomocí algoritmu A1, dostaneme x2 = 0 a x1 x2 = 0 6= c. Když využijeme algoritmus A2, je x2 = 0, 0000100000000 a x1 x2 = 1 = c. Vidíme, že tomto případě je algoritmus A2 stabilnější než algoritmus A1.

Na závěr si shrňme, jak budeme postupovat při řešení praktického problému. — Ke sledovanému problému musíme přiřadit odpovídající matematický model. To znamená zadat jednoznačný a srozumitelný funkční vztah mezi zadanými a hledanými matematickými objekty (čísly, maticemi, funkcemi, křivkami). 13

— Převedeme matematickou úlohu, pokud není numerická, na numerickou (viz definice 1.1). — Najdeme vhodný algoritmus řešení a analyzujeme ho. — Provedeme rozbor přesnosti získaných výsledk˚ u.

Některá fakta z lineární algebry a funkcionální analýzy V dalších kapitolách se nám vedle poznatků získaných v základních kurzech algebry a matematické analýzy budou hodit ještě některé další skutečnosti, které si zde ve stručnosti uvedeme. V řadě algoritmů pro řešení úloh lineární algebry lze těžit z toho, že matice, se kterou pracujeme, má určitou speciální strukturu. Zajímavé pro nás budou zejména následující typy matic: Definice 1.6 Čtvercová matice A = (aij )ni,j=1 , aij ∈ R, je ostře diagonálně dominantní, když pro všechna i = 1, . . . , n je |aii | >

n X

|aik |.

(1.1)

k=1, k6=i

Symetrická, když pro všechna i, j = 1, . . . , n platí aij = aji . Ortogonální, když pro všechna i, j = 1, . . . , n platí n X

2

aik akj = δij .

(1.2)

k=1

Pozitivně definitní, když n X n X

aij xi xj > 0 pro libovolný vektor ~x = (x1 , . . . , xn ).

i=1 j=1

14

(1.3)

Čtvercová matice je (p, q)-pásová, když existují čísla p, q < n tak, že pro prvky matice A platí aij = 0 pokud j > i + p nebo i > j + q.

(1.4)

Číslo p + q + 1 je pak šířka pásu. Příklad Matice



 2 −1 0 5 2  A =  −1 0 2 3

je ostře diagonálně dominatní, symetrická, pozitivně definitní, pásová. Významnou roli při rozhodování o konvergenci iteračních metod pro řešení soustav lineárních rovnic hrají vlastní čísla. Definice 1.7 Nechť A je čtvercová matice řádu n, I jednotková matice řádu n. Řekneme, že λ je vlastním číslem matice A, když homogenní soustava (A − λI)~v = ~0,

(1.5)

má netriviální řešení ~v = (v1 , v2 , . . . , vn )T . Poznámka Soustava (1.5) má nenulové řešení, právě když det(A − λI) = 0. Uvedený determinant představuje tzv. charakteristický polynom p(λ) = (−1)n λn + b1 λn−1 + · · · + bn−1 λ + bn .

(1.6)

Pokud A je matice řádu n, je p(λ) polynom n-tého stupně a má tedy právě n kořen˚ u λ1 , λ2 , . . . , λn , které nemusejí být navzájem různé. Pak se ale může stát, že počet lineárně nezávislých vlastních vektorů je menší než je řád matice (viz [3] strana 107). 2

½ δij =

0 1

pro i 6= j jinak

15



 1 1 0 Příklad Stanovte vlastní čísla a vlastní vektory matice A =  0 2 0  . 0 0 3 Řešení:

¯ ¯ ¯ ¯1 − λ 1 0 ¯ ¯ 2−λ 0 ¯¯ = (2 − λ)(λ − 3)(λ − 1) = 0. |A − λI| = ¯¯ 0 ¯ 0 0 3 − λ¯

Odtud dostaneme, že λ1 = 1, λ2 = 2 a λ3 = 3 jsou vlastní čísla matice A. Dále spočteme vlastní vektory. Pro λ1 = 1 budeme řešit soustavu Aλ1 = 0 tzn.   µ ¶ 0 1 0 |0  0 1 0 | 0  −→ 0 1 0 | 0 , 0 0 2 |0 0 0 2 |0 Vidíme, že obecným řešením soustavy je vektor ~v1 = (t, 0, 0)T , t ∈ R. Pro λ2 = 2 je



 µ ¶ −1 1 0 | 0  0 0 0 | 0  −→ −1 1 0 | 0 , 0 0 1 |0 0 0 1 |0

~v2 = (r, r, 0)T , r ∈ R. Pro λ3 = 3 je



 µ ¶ −2 1 0 | 0 1 0 |0  0 −1 0 | 0  −→ −2 , 0 −1 0 | 0 0 0 0 |0

~v3 = (0, 0, s)T , s ∈ R. Poznámka Dále připomeňme, že symetrická matice má pouze reálná vlastní čísla. Pokud tato vlastní čísla jsou ještě navíc navzájem různá, jsou k nim příslušné vlastní vektory ortogonální. Symetrická pozitivně definitní matice má všechna vlastní čísla reálná a kladná. Když ~v je vlastní vektor matice A příslušný k vlastnímu číslu λ, pak platí ~v T A~v λ= T . ~v ~v 16

Definice 1.8 Pokud ~x není vlastním vektorem matice A, pak podíl ~x T A~x ~xT ~x

(1.7)

nazýváme Rayleighovým podílem. Problematika numerického výpočtu vlastních čísel matice je spojena s pojmem podobné matice. Definice 1.9 Řekneme, že matice A je podobná matici B, když existuje regulární matice T tak, že T −1 AT = B. Píšeme pak A ∼ B. Při konstrukci numerických metod pro výpočet vlastních čísel matice se využívá následující věta. Věta 1.1 Když A a B jsou podobné matice, pak mají stejná vlastní čísla. Důkaz. det(B − λI) = det (T −1 (A − λI)T ) = det(A − λI). Poznámka Opak platit nemusí. Totiž matice, které mají stejná vlastní čísla, ještě nemusejí být podobné. Např. matice ¯    1 1| 0 1 0 0 ¯¯ A =  0 2 0 ¯¯ a B =  0 2 | 0  0 0 2 0 0 2 ¯ sice mají stejná vlastní čísla, ale nejsou podobné, protože jsou tvořeny navzájem různými Jordanovými bloky

Následující systém definic a vět představuje obecný aparát, který využijeme při odvozování numerických metod řešení daných problémů. 17

Definice 1.10 Množinu X nazveme lineárním prostorem, když pro každé dva prvky x, y ∈ X a libovolné konstanty α, β ∈ R platí αx + βy ∈ X. Příklad Lineární prostory tvoří Rn — množina n rozměrných sloupcových vektorů, které mají reálné složky Mn — množina čtvercových matic řádu n, Pn — množina polynomů n-tého stupně, C n (a, b) — množina funkcí n-krát spojitě diferencovatelných na (a, b), Rb Lp (a, b) — množina funkcí, pro které je (Lebesgueův) integrál a |f (x)|p dx konečný. Definice 1.11 Nechť X je lineární prostor. Zobrazení ||.|| : X → R, které pro všechna x, y ∈ X, α ∈ R splňuje ||x|| ≥ 0, (||x|| = 0 ⇔ x = 0), ||αx|| = |α| ||x||, ||x + y|| ≤ ||x|| + ||y||, se nazývá norma. Lineární prostor s normou budeme nazývat normovaným lineárním prostorem. Příklad Pro ~x ∈ Rn představuje k~xk∞

v u n uX x2i = max {|xi |} nebo k~xkE = t i=1,...,n

i=1

normu vektoru. V prostoru Mn můžeme normu matice definovat jedním z následujících způsobů: Řádková norma ||A||∞ = max

i=1,...,n

Sloupcová norma ||A||S = max

j=1,...,n

18

n X

|aij |.

(1.8)

|aij |.

(1.9)

j=1

n X i=1

Eukleidovská norma

v uX n u n X t ||A||E = ( a2ij ).

(1.10)

i=1 j=1

Pro f ∈ C 0 (a, b) lze vzít ||f (x)||∞ = max |f (x)|. x∈

(1.11)

Pro f ∈ Lp (a, b), p = 1, 2, . . . je Z

b

||f (x)||p = (

1

|f (x)|p dx) p .

(1.12)

a

Poznámka Říkáme, že dvě normy k.k1 a k.k2 v lineárním prostoru X jsou ekvivalentní, právě když existují kladné konstanty c1 , c2 ∈ R tak, že pro všechna x ∈ X platí c1 kxk1 ≤ kxk2 ≤ c2 kxk1 . Protože v konečnědimenzionálních prostorech jsou všechny normy ekvivalentní, jsou ekvivalentní také normy (1.8), (1.9), (1.10). Poznámka Když A, B ∈ Mn a ~x ∈ Rn , pak pro odpovídající si normy platí ||AB|| ≤ ||A|| ||B||.

(1.13)

||A~x|| ≤ ||A|| ||~x||.

(1.14)

Definice 1.12 Nechť X je lineární prostor. Zobrazení z X do R, které pro všechna x, y ∈ X, α, β ∈ R splňuje (x, x) ≥ 0, ((x, x) = 0 ⇔ x = 0), (x, y) = (y, x), (αx + βy, z) = (αx, z) + (βy, z), se nazývá skalární součin. Lineární prostor se skalárním součinem je tzv. Hilbertův prostor. 19

Příklad (Spojitý případ) Když funkce f, g ∈ L2 (a, b), je možné definovat skalární součin3 Z b (f, g) = f (x)g(x) dx. (1.15) a

(Diskrétní případ) Když funkce f a g jsou definovány v uzlových bodech x0 , x1 , . . . , xn ∈ ha, bi, je skalární součin (f, g) =

n X

f (xj )g(xj ) dx.

(1.16)

j=0

Poznámka V Hilbertově prostoru X lze definovat normu vztahem p ||x|| = (x, x), x ∈ X. √ Příklad Nechť i = −1 je imaginární jednotka. V L2 (−π, π) stanovte normu funkce f (x) = eix . √ Rπ Řešení: keix k2 = −π eix e−ix dx = 2π, keix k = 2π. Definice 1.13 Nechť X je Hilbertův prostor. Řekneme, že prvky xi , xj ∈ X jsou ortogonální, když ½ = 0 pro i 6= j, (xi , xj ) 6= 0 pro i = j. Speciálně když platí (xi , xj ) = δij , jsou prvky xi , xj ∈ X ortonormální. Poznámka Ve spojitém případě je ortogonalita závislá na volbě váhové funkce. V diskrétním případě ještě navíc na volbě tabulkových bodů. Příklad Jako příklad ortogonálních soustav mohou sloužit následující systémy algebraických a trigonometrických polynomů: a) Funkce {xj }, j = 0, 1, . . . , jsou ortogonální na intervalu (−∞, ∞). b) Funkce {1, cos x, sin x, cos 2x, sin 2x, . . . } jsou ortogonální na h−π, πi. 3

Funkce g(x), je komplexně sdružená k funkci g(x).

20

c) Funkce {eijx }, i2 = −1, j = 0, ±1, ±2. . . . , jsou ortogonální na intervalu h−π, πi. d) Legenderovy polynomy P0 (x) = 1, P1 (x) = x, (n + 1)Pn+1 (x) = (2n + 1)xPn (x) − nPn−1 (x), které můžeme také definovat vztahem 1 dn (x2 − 1)n Pn (x) = n , 2 n! dxn jsou ortogonální na intervalu h−1, 1i. Ukažme si nyní, že funkce ϕj = eijx (j = 0, ±1, . . . ) jsou opravdu na intervalu h−π, πi ortogonální. Spojitý případ: ½ Z π 0 pro j 6= k, ijx −ikx (ϕj (x), ϕk (x)) = e e dx = 2π pro j = k. −π Diskrétní případ: Když xs = h−π, πi, pak

2πs , n+1

s = 0, 1, . . . n, jsou uzlové body z intervalu

(ϕj (x), ϕk (x)) =

n X

2πs

ei(j−k) n+1

s=0 2π

je n-tým částečným součtem geometrické řady, která má kvocient q = ei(j−k) n+1 . j−k je celé, je q = 1 a (ϕj (x), ϕk (x)) = n + 1. V opačném případě je Pokud n+1 (ϕj (x), ϕk (x)) = 0. Definice 1.14 Nechť X je normovaný prostor. Řekneme, že posloupnost prvků xn ∈ X konverguje k x ∈ X, když ||xn − x|| → 0, xn → ∞. Definice 1.15 Nechť X je normovaný prostor. Posloupnost {xn } je cauchyovská v X, když ||xm − xn || → 0, n, m → ∞. Definice 1.16 Normovaný prostor, ve kterém je každá cauchyovská posloupnost konvergentní, je úplný normovaný neboli Banachův prostor. 21

Příklad Prostor C (0) (a, b) s normou ||.||∞ je úplný. Konvergence posloupnosti spojitých funkcí v normě ||.||∞ je totiž ekvivalentní se stejnoměrnou konvergencí těchto funkcí. Limita posloupnosti spojitých funkcí je tedy opět spojitá funkce. Prostor Lp (a, b) s normou ||.||p je úplný. Definice 1.17 Nechť X, Y jsou lineární prostory. Zobrazení A : X → Y je lineární operátor, když pro všechna x, y ∈ X, α, β ∈ R platí f (αx + βy) = αf (x) + βf (y).

Definice 1.18 Nechť X, Y jsou normované prostory. Řekneme, že lineární operátor A : X → Y je omezený, když existuje konstanta c > 0 taková, že ||Ax|| ≤ c||x||. Nejmenší z konstant, která splňuje uvedenou nerovnost, je tzv. norma operátoru A. Platí ||A|| = sup ||Ax||. ||x||=1

Definice 1.19 Lineární operátor A : X → Y je spojitý, když pro všechna xn ∈ X, xn → x je Axn → Ax. Věta 1.2 omezený.

V normovaném prostoru je lineární operátor spojitý, když je

Důkaz. Pokud je A omezený operátor, je ||Axn − Ax|| = ||A(xn − x)|| ≤ ||A|| ||xn − x|| → 0 pro ||xn − x|| → 0.

Definice 1.20 Nechť U je podprostor normovaného prostoru X. Operátor A : U → X nazveme kontrakcí, když existuje kontanta q ∈ h0, 1) tak, že ||Ax − Ay|| ≤ q||x − y|| 22

(1.17)

pro všechna x, y ∈ U. Definice 1.21 Nechť X je Banachův prostor. Řekneme, že x ∈ X je pevným bodem operátoru A : X → X, když Ax = x. Věta 1.3 (Banach) Nechť U je úplný podprostor normovaného prostoru X a operátor A : U → U je kontrakce, pak A má právě jeden pevný bod. Důkaz. Předpokládejme, že A je kontrakce (viz definice 1.20), x0 ∈ U je pevně dáno a pro n = 1, 2, . . . je xn+1 = Axn . Pak ||xn+1 − xn || = ||Axn − Axn−1 || ≤ q ||xn − xn−1 || = q ||Axn−1 − Axn−2 || ≤ ≤ q 2 ||xn−1 − xn−2 || ≤ · · · ≤ q n ||x1 − x0 ||. Pro m > n dostaneme ||xn − xm || ≤ ||xn − xn+1 || + ||xn+1 − xn+2 || + · · · + ||xm−1 − xm || ≤ ≤ (q n + q n+1 + · · · + q m−1 )||x1 − x0 ||. (1.18) Protože A je kontrakce, máme pro n → ∞ ||xn − xm || ≤

qn ||x1 − x0 || → 0. 1−q

To znamená, že posloupnost {xn } je cauchyovská v úplném prostoru U, a tedy i konvergentní. Ze spojitosti operátoru A máme x = lim xn+1 = lim Axn = Ax. n→∞

n→∞

Je tedy x = limn→∞ xn+1 pevným bodem operátoru A. Tento bod je právě jeden. Kdyby existovaly dva takové pevné body x a y, pak 0 6= ||x − y|| = ||Ax − Ay|| ≤ q||x − y||. To je možné, jen když q ≥ 1. To je ale spor s tím, že konstanta q < 1. 23

Poznámka Banachova věta nejenže udává podmínku pro existenci jediného pevného bodu operátoru A, ale její důkaz poskytuje i návod, jak pevný bod najít. Věta 1.4 Když A je kontrakce úplného normovaného prostoru U do sebe, pak aproximace x(k+1) = Ax(k) konvergují pro libovolné x0 k pevnému bodu x operátoru A. Jestliže q je kontrakční konstanta, pak pro všechna n ∈ N dostaneme apriorní odhad qn ||xn − x|| ≤ ||x1 − x0 || (1.19) 1−q a aposteriorní odhad ||xn − x|| ≤

q ||xn − xn−1 ||. 1−q

(1.20)

Důkaz. Existence pevného bodu plyne z Banachovy věty 1.3. Apriorní odhad obdržíme ze vztahu (1.18), když m → ∞. Aposteriorní odhad plyne z apriorního, když začneme od x0 = xn−1 . Poznámka Apriorní odhad slouží jako horní mez při určování počtu kroků potřebných k dosažení požadované přesnosti. Když je ε požadovaná přesnost, pak z (1.19) qn ||x1 − x0 || < ε, 1−q 1−q qn < ε, ||x1 − x0 || n≥

(1 − q) ln ε˜ , kde ε˜ = ε. ln q ||x1 − x0 ||

Vidíme, že čím je menší q, tím je k dosažení vyšší přesnosti zapotřebí méně kroků. Na druhé straně aposteriorní odhad poskytuje lepší odhad přesnosti než odhad apriorní. 24

Věta 1.5 Když X je Banachův prostor, I : X → X je identita a H : X → X je omezený operátor takový, že ||H|| < 1, pak rovnice (I − H)x = z má pro každé z ∈ X právě jedno řešení x ∈ X. Posloupnost aproximací xn+1 = Hxn + z pak konverguje pro libovolné x0 ∈ X k řešení x a platí apriorní odhad ||xn − x|| ≤

||H||n ||x1 − x0 || 1 − ||H||

a aposteriorní odhad ||xn − x|| ≤

||H|| ||xn − xn−1 ||. 1 − ||H||

Důkaz. Pro libovolné z ∈ X pevné definujme operátor A : X → X vztahem Ax = Hx + z. Pro všechna x, y ∈ X je ||Ax − Ay|| = ||Hx − Hy|| ≤ ||H|| ||x − y||. Protože q = ||H|| < 1, je A kontrakce. Když zvolíme x0 = z, pak pro postupné aproximace máme ||xn || ≤

n X i=0

i

||H z|| ≤

n X

||H i || ||z|| ≤

i=0

1 ||z||. 1 − ||H||

Protože pro n → ∞ konverguje xn → x = (I −H)−1 z, dostaneme pro všechna z∈X ||z|| . ||(I − H)−1 z|| ≤ 1 − ||H||

25

Definice 1.22 Nechť X je normovaný prostor, U ⊂ X a x ∈ X. Prvek v ∈ U nazveme nejlepší aproximací x vzhledem k U, když ||x − v|| = inf ||x − u||. u∈U

(1.21)

Poznámka Prvek v ∈ U je nejlepší aproximací x ∈ X, právě když pro všechna u ∈ U je (x − v, u) = 0. (1.22)

Poznámka Ke každému prvku normovaného prostoru existuje nejlepší aproximace vzhledem k nějakému jeho konečnědimenzionálnímu podprostoru. Věta 1.6 Nechť U je konečnědimenzionální podprostor Hilbertova prostoru X a u1 , . . . , un jeho báze. Prvek v=

n X

αi ui

(1.23)

i=1

je nejlepší aproximací prvku x vzhledem k U, právě když koeficienty α1 , . . . , αn vyhovují soustavě n X

αi (ui , uj ) = (x, uj ), j = 1. . . . , n.

(1.24)

i=1

Důkaz. Vztah (1.24) obdržíme dosazením (1.23) do (1.22) Věta 1.7 Nechť U je podprostor Hilbertova prostoru X a u1 , . . . , un jeho báze. Nejlepší aproximaci v prvku u pak lze psát n X (x, ui )ui . v= i=1

Důkaz. Tvrzení plyne z věty 1.6. 26

(1.25)

Kapitola 2 Řešení soustav lineárních rovnic Řešení soustav lineárních rovnic je velice frekventovaný problém. Na úlohu najít řešení soustavy lineárních rovnic vedou nejr˚ uznější optimalizační úlohy nebo k ní dospějeme při numerickém řešení parciálních diferenciálních rovnic. To podtrhuje význam tohoto odstavce. Přestože soustavy lineárních rovnic mohou mít obecnější tvar, omezíme se pouze na soustavy lineárních rovnic s regulární čtvercovou maticí A řádu n A~x = ~b, (2.1) kde ~x je n-členný sloupcový vektor řešení a ~b je n-dimenzionální sloupcový vektor absolutních členů. Připomeňme si, že matice A je regulární, právě když det A 6= 0. V takovém případě soustava A~x = ~b má jediné řešení. Metody, pomocí kterých toto řešení budeme hledat, rozdělíme do dvou skupin: • Metody přímé (Gaussova eliminační metoda, LU rozklad). Tyto metody poskytují přesné řešení v konečném počtu kroků, pokud během výpočtu nezaokrouhlujeme. • Metody nepřímé tzv. iterační (Gaussova-Seidelova metoda, Jacobiova metoda). Těmito metodami získáme pouze aproximativní řešení. Obecně nelze říci, který z uvedených typů pro který okruh řešených úloh je nejvýhodnější použít. Pro první orientaci lze zhruba uvést, že soustavy s plnou maticí se spíš řeší pomocí přímých metod a soustavy s řídkou maticí, v nichž převládají nulové prvky, zase pomocí iteračních metod. 27

Přímé metody

a) Gaussova eliminační metoda (GEM) Nejprve si metodu odvodíme. Rozepišme soustavu (2.1) na a011 x1 + a012 x2 + · · · + a01n xn = a01,n+1 ,

(2.2)

a021 x1 + a022 x2 + · · · + a02n xn = a02,n+1 , ........................... a0n1 x1 + a0n2 x2 + · · · + a0nn xn = a0n,n+1 . Nyní budeme postupně eliminovat prvky pod hlavní diagonálou matice soustavy. Nejprve 1. rovnici využijeme k eliminaci 1. neznámé ze zbývajících n − 1 a0 rovnic: 1. řádek vynásobíme čıslem m021 = − a21 a přičteme k 2. řádku, 0 1. řádek vynásobený číslem získáme soustavu

m031

=

a0 − a31 0 11

11

přičteme ke 3. řádku a tak dále, až

a011 x1 + a012 x2 + . . . + a01n xn = a01,n+1 , a122 x2 + . . . + a12n xn = a12,n+1 , ................................. a1n2 x2 + . . . + a1nn xn = a1n,n+1 . Obecně v k-tém kroku (k = 1, 2, . . . , n−1) vynulujeme prvky v k-tém sloupci ak−1

ik přičteme k itím, že k-tý řádek vynásobený multiplikátorem mk−1 = − ak−1 ik kk

tému řádku (i = k + 1, . . . n). Matici A tak pomocí ekvivalentních úprav převedeme na horní trojúhelníkovou matici. Obdržíme soustavu a011 x1 + a012 x2 + a122 x2 +

... ... a233 x3 + . . .

28

+a01n xn = a01,n+1 , +a12n xn = a12,n+1 , +a23n xn = a23,n+1 , .................. n−1 an−1 nn xn = an,n+1 ,

z níž zpětnou eliminací získáme hodnoty neznámých xn , xn−1 , . . . , x1 . Algoritmus (GEM) Vstup: n, A = (a0ij ), ~b = (a0i,n+1 ) pro k = 1, 2, . . . , n − 1 pro i = k + 1, k + 2, . . . , n ak−1 k−1 mk−1 = − aik 6= 0) k−1 (pokud akk ik kk pro j = k + 1, k + 2, . . . , n + 1 k−1 akij = ak−1 + mk−1 ij ik akj

pro i = n, n − 1, . . . , 1 P n i−1 1 xi = ai−1 (ai−1 i,n+1 − j=i+1 aij xj ) ii

Výstup: ~x = (x1 , x2 , . . . , xn ) Poznámka Protože platí n n(n + 1) X 2 n(n + 1)(2n + 1) k = k= , , 2 6 k=1 k=1

n X

dojdeme k závěru, že v přímém chodu GEM, kdy matici převádíme na trojúhelníkový tvar, potřebujeme k výpočtu multiplikátorů mk−1 ik n−1 X (n − 1)n n(n − 1) (n − k) = (n − 1)n − = 2 2 k=1

dělení a k výpočtu prvků akij potřebujeme n−1 n−1 X X n3 − n (n − k)(n − k + 1) = n(n + 1) − k(2n + 1) + k 2 = 3 k=1 k=1

násobení a stejně tolik sčítání. Při zpětném chodu GEM, tj. eliminaci proměnných, je třeba n2 + n 2 29

násobení a

n2 − n 2 3 sčítání. Celkem v našem případě tedy potřebujeme n3 + n2 − n3 násobení a 3 2 dělení a n3 + n2 − 5n sčítání. Z hlediska doby potřebné k realizaci algoritmu 6 jsou operace násobení a dělení v podstatě totožné. Protože časová náročnost násobení je několikanásobně vyšší než ta, kterou vyžaduje operace sčítání, budeme o efektivnosti výpočtu rozhodovat v závislosti na počtu násobení. Definice 2.1 Řekneme, že pro x → a je funkce f stejného řádu jako funkce g, když existuje konečná limita f (x) lim . x→a g(x) Píšeme pak f (x) = O(g(x)), x → a. 3

Poznámka K realizaci algoritmu GEM je potřeba n3 + O(n2 ) násobení. To znamená, že když v soustavě počet rovnic zdvojnásobíme, vzroste doba výpočtu osmkrát. Příklad Gaussovou eliminační x1 + x1 + 2x1 +

metodou řešte soustavu 2x2 − 2x3 = 1, x2 + x3 = 3, 2x2 + x3 = 5.

Řešení: Rozšířenou matici soustavy uvedeme na trojúhelníkový tvar. V našem případě řádky postupně násobíme multiplikátory m021 = −1, m031 = −2; m132 = −2. Dostaneme       1 2 −2 | 1 1 2 −2 | 1 1 2 −2 | 1  1 1 1 | 3  −→  0 −1 3 | 2  −→  0 −1 3 | 2  . 2 2 1 | 5 0 −2 5 | 3 0 0 −1 | −1 Zpětnou eliminací spočteme řešení: x3 = 1, x2 = −2 + 3x3 = 1, x3 = 1 − 2x2 + 2x3 = 1. 30

Poznámka Při řešení pomocí GEM jsme vyloučili případ, kdy v k-tém kroku k−1 je hlavní prvek ak−1 kk = 0. Pokud taková situace nastane a některé akk se v počítači zobrazí jako nula, můžeme zkusit použít Gaussovu eliminační metodu s výběrem hlavního prvku. Máme tři možnosti: i) GEM se sloupcovým výběrem — za hlavní prvek v k-tém kroku bereme největší prvek v k-tém sloupci. Vybíráme mezi řádky, z nichž jsme dosud nevzali vedoucí prvek. Hlavní prvek je v řádku p a platí pro něj |apk | = max |ak−1 ik |. i

ii) GEM s řádkovým výběrem — za hlavní prvek v k-tém kroku bereme největší prvek v k-tém řádku. Vybíráme mezi sloupci, z nichž jsme dosud nevzali vedoucí prvek. Hlavní prvek je ve sloupci q a platí pro něj |akq | = max |ak−1 kj |. j

iii) GEM s úplným výběrem — za hlavní prvek v k-tém kroku bereme největší prvek v tom řádku a sloupci, ve kterém jsme dosud nevybrali vedoucí prvek. Hlavní prvek se nachází v p-tém řádku a q-tém sloupci a platí pro něj |apq | = max |ak−1 ij |. ij

Příklad Gaussovou eliminační metodou řešte soustavu x1 + 2x2 − 2x3 = 1, x1 + 2x2 + x3 = 4, 2x1 + 2x2 + x3 = 5. Řešení:



   1 2 −2 | 1 1 2 −2 | 1  1 2 1 | 4  −→  0 0 3 | 3  . 2 2 1 | 5 0 −2 5 | 3

GEM nelze přímo aplikovat, protože ve druhém kroku je a122 = 0. Pokusme se řešit soustavu pomocí GEM se sloupcovým výběrem. V 1. sloupci p˚ uvodní matice má největší absolutní hodnotu prvek a031 . Proto přesuneme 31

3. řádek na začátek a dále už pokračujeme standardní Gaussovou eliminační metodou:     1 2 −2 | 1 2 2 1 | 5  1 2 1 | 4  →  1 2 −2 | 1  , 1 2 1 | 4 2 2 1 | 5       2 2 1 | 5 2 2 1 | 5 2 2 1 | 5  1 2 −2 | 1  −→  0 1 − 5 | − 3  −→  0 1 − 5 | − 3  , 2 2 2 2 0 1 12 | 32 0 0 62 | 26 1 2 1 | 4 x3 = 1, x2 = 1, x1 = 1. Jak bylo řečeno v předchozím odstavci, ne každé číslo se v počítači zobrazí přesně. Následující dva příklady jsou ukázkou toho, jak GEM s výběrem hlavního prvku může přispět ke snížení kumulování zaokrouhlovacích chyb v průběhu výpočtu. Příklad Gaussovou eliminační metodou v M (10, 2) řešte soustavu x1 + 200x2 = 100, x1 + x2 = 1. Řešení: Pokud během výpočtu nezaokrouhlujeme, obdržíme ¶ µ ¶ µ 1 200 | 100 1 200 | 100 −→ , 0 −199 | −99 1 1 | 1 99 = 0, 4974 . . . , 199 99 100 x1 = 100 − 200 = = 0, 5025 . . . 199 199 x2 =

Když během výpočtu zaokrouhlujeme, pak v M (10, 2) bude µ ¶ µ ¶ 1 200 | 100 1 200 | 100 −→ , 1 1 | 1 0 −200 | −99 99 . = 0, 50, 200 x1 = 100 − 200 · 0, 50 = 0. x2 =

32

Vidíme, že zaokrouhlovací chyby úplně znehodnotily výpočet. Příklad Soustavu z předchozího příkladu řešte v M (10, 2) Gaussovou eliminační metodou s řádkovým výběrem. Řešení: Protože největší prvek v 1. řádku se nachází ve 2. sloupci, zaměníme v původní matici první a druhý sloupec a pokračujeme eliminační metodou ¶ µ ¶ µ ¶ µ 200 1 | 100 200 1 | 100 1 200 | 100 → −→ . 1 1 | 1 0 0, 995 | 0, 455 1 1 | 1 Při zaokrouhlení v M (10, 2) 0, 455 . = 0, 46, 0, 995 100 − 0, 46 . x2 = = 0, 50. 200 x1 =

Při použití Gaussovy eliminační metody s řádkovým výběrem jsme tak docílili lepších výsledk˚ u než při řešení prostřednictvím Gaussovy eliminační metody bez výběru.

b) Metoda LU-rozkladu Gaussova eliminační metoda není nic jiného než převedení původní matice soustavy prostřednictvím lineární transformace na horní trojúhelníkovou matici, ze které pak zpětnou eliminací stanovíme řešení. Když A je p˚ uvodní matice soustavy, U je příslušná trojúhelníková matice a T = Tn−1 Tn−2 . . . T2 T1 , kde pro k = 1, 2, . . . , n − 1 je   1 0 ... 0 0 0 ... 0  0 1 ... 0 0 0 ... 0      ... ... . . . . . . . . . . . . . . . . . .   k−1 0 0 ... 0  0 mk+1,k 1 Tk =  ,  0 k−1  0 1 0 ... 0  0 mk+2,k 0    ... ... ... ... ... ... ... ...  0 0 0 ... 1 0 0 mk−1 n,k 33

pak proces převedení matice A na trojúhelníkovou matici lze zapsat jako T A = U. Matice T je regulární, a proto k ní existuje matice inverzní L = T −1 . Platí tedy A = T −1 U = LU. Tím jsme vyjádřili matici A jako  l11  l21  L=  l31  ... ln1

součin dolní trojúhelníkové matice  0 0 ... 0 l22 0 . . . 0   l32 l33 . . . 0   ... ... ... ...  ln2 ln3 . . . lnn

a horní trojúhelníkové matice  u11 u12  0 u22  0 U =  0  ... ... 0 0

u22 u23 u33 ... 0

... ... ... ... ...

u1n u2n u3n ... unn

   .  

Proto zde také hovoříme o tzv. LU-rozkladu. Uvedený rozklad je ekvivalentní s Gaussovou eliminační metodou, a proto je možné ho realizovat, právě když lze realizovat GEM. µ

0 −5 Příklad LU-rozklad nelze aplikovat přímo, když A = 3 0 0 matici pivot a11 = 0, a proto nelze stanovit multiplikátor m021 . µ

¶ . V této

¶ 2 5 Příklad Pro matici A = LU-rozklad existuje. Lehce ověříme, že 4 1 µ ¶µ ¶ µ ¶µ ¶ 1 0 2 5 2 0 1 52 A= nebo také A = . 2 1 0 −9 4 −1 0 9 34

Poznámka Z posledního příkladu je vidět, že rozklad A = LU není určen jednoznačně. Ze soustavy A = LU, která má n2 rovnic, je totiž třeba spočítat n(n+1) neznámých l11 , l21 , . . . , lnn , u11 , u12 , . . . , unn . To znamená, že řešení obsahuje n volných nezámých. Jednoznačnost rozkladu si zajistíme třeba tak, že zadáme prvky na hlavní diagonále matice L. Bývá zvykem volit lii = 1, i = 1, 2, . . . n. Definice 2.2 Nechť A je regulární matice řádu n. Řekneme, že dolní trojúhelníková matice L s jedničkami na hlavní diagonále a horní troúhelníková matice U tvoří LU-rozklad matice A, když platí A = LU. Příklad Najděte vztah mezi prvky LU-rozkladu pro matici třetího řádu. Řešení: 

    a11 a12 a13 1 0 0 u11 u12 u13  a21 a22 a23  =  l21 1 0   0 u22 u23  . a31 a32 a33 l31 l32 1 0 0 u33 u11 = a11 , l21 u11 = a21 =⇒ l21 = ua21 , 11 a31 l31 u11 = a31 =⇒ l31 = u11 , u12 = a12 , l21 u12 + u22 = a22 =⇒ u22 = a22 − l21 u12 , 31 u12 l31 u12 + l32 u22 = a32 =⇒ l32 = a32 −l , u22 u13 = a13 , l21 u13 + u23 = a23 =⇒ u23 = a32 − l21 u13 , l31 u13 + l32 u23 + u33 = a33 =⇒ u33 = a33 − l31 u13 − l32 u23 . Algoritmus (LU-rozklad) Vstup: n, A = (aij ) pro j = 1, 2, . . . , n pro i = 1, 2, . .P .,j uij = aij − i−1 r=1 lir urj pro i = j + 1, j +P 2, . . . , n 1 lij = ujj (aij − j−1 s=1 lis usj ) (pokud ujj 6= 0) Výstup: L = (lij ), U = (uij ) 35

K realizaci uvedeného algoritmu je třeba

n3 3

+ O(n2 ) násobení.

Věta 2.1 Nechť A = LU je LU-rozklad matice A. Vektor ~x je řešením soustavy (2.1), právě když je řešením soustavy L~y = ~b

(2.3)

U~x = ~y .

(2.4)

Důkaz. Nechť A = LU. Pak A~x = ~b ⇔ L(U~x) = ~b ⇔ L~y = ~b, U~x = ~y . Poznámka Ze soustav (2.3), (2.4) můžeme přímo eliminovat vektory řešení ~y a ~x, protože L a U jsou trojúhelníkové matice. Příklad Pomocí LU-rozkladu řešte soustavu x1 + 2x2 − 2x3 = 1, x1 + x2 + x3 = 3, 2x1 + 2x2 + x3 = 5. Řešení: A = LU, kde       1 2 −2 1 0 0 1 2 −2 A =  1 1 1  , L =  1 1 0  , U =  0 −1 3  . 2 2 1 2 2 1 0 0 −1 Eliminace ze soustavy L~y = ~b : 

 1 0 0 | 1  1 1 0 | 3 , 2 2 1 | 5

y1 = 1, y2 = 2, y3 = −1. Zpětná eliminace z U~x = ~y :   1 2 −2 | 1  0 −1 3 | 2  , 0 0 −1 | −1 36

x3 = 1, x2 = 1, x1 = 1. Poznámka Metodu LU-rozkladu je výhodné použít v případě, že máme řešit víc soustav se stejnou maticí soustavy A. Rozklad se pak počítá jen jednou. Při řešení celé řady úloh, mezi které patří například řešení parciálních diferenciálních rovnic metodou konečných prvk˚ u, se setkáváme s maticemi, které mají speciální tvar. Na závěr si proto uveďme dva příklady algoritm˚ u LUrozkladu pro takové matice. i) Pokud je matice A symetrická pozitivně definitní (viz definice 1.6), dá se ukázat, že existuje horní trojúhelníková matice U s kladnými prvky na hlavní diagonále tak, že A = U T U. Toho využijeme v následujícím algoritmu. Algoritmus (Choleského pro matici A, která je uložená ve sloupcích.) Vstup: n, A = (aij ) pro j = 1, n q2, . . . , P 2 ujj = ajj − j−1 r=1 ujr pro i = j + 1, j +P 2, . . . , n 1 uij = ujj (aij − j−1 s=1 uis ujs ) Výstup: U = (uij ) V Choleského algoritmu je zapotřebí

n3 6

+ O(n2 ) násobení.

Příklad Choleského algoritmus použijte k řešení soustavy x1 + x2 + x3 = 2, x1 + 5x2 + 5x3 = 6, x1 + 5x2 + 14x3 = 6. 

 1 1 1 Řešení: Označme A =  1 5 5  , ~b = (2, 6, 6)T . 1 5 14 37



   1 0 0 1 1 1 Platí A = U T U, když U T =  1 2 0  , U =  0 2 2  . 1 2 3 0 0 3 Řešení soustavy U T ~y = ~b :   1 0 0 | 2  1 2 0 | 6  , ~y = (2, 2, 0)T . 1 2 3 | 6 Ze soustavy U~x = ~y máme   1 1 1 | 2  0 2 2 | 2  , ~x = (1, 1, 0)T . 0 0 3 | 0

ii) Dalším typem speciálních matic, které se v praktických problémech často vyskytují, jsou pásové matice. (Pásová matice viz definice 1.6.) V takovém případě stačí pracovat pouze s nenulovými prvky matice. Sníží se pak počet operací potřebných k realizaci výpočtu. Přednosti tohoto postupu vyniknou, když p << n. Algoritmus (LU-rozklad pro (p, p)-pásovou matici) Vstup: n, p, A = (aij ) pro i = 1, 2, . . . , n γ = min(n, i + p) pro j = i, i + 1, . . . , γ α = max(1,P j − p) uij = aij − i−1 r=α lir ujr pro j = i + 1, i + 2, . . . , γ α = max(1, j −Pp) lji = u1ii (aji − i−1 r=α ljr uri ) Výstup: L = (lij ), U = (uij ) Počet násobení a dělení potřebných k realizaci algoritmu je np(p + 1). 38

Příklad Pomocí LU-rozkladu pro třídiagonální matici řešte soustavu −x1 + 2x1 +

x2 x2 + 2x3 3x2 − x3 + 3x4 4x3 + x4



−1  2 Řešení: Označme A =   0 0

= = = =

0, 5, 5, 5.

 1 0 0 1 2 0   , ~b = (0, 5, 5, 5)T . 3 −1 3  0 4 1



  1 0 0 0 −1  −2 1   0 0 0  aU = Platí, že A = LU, když L =   0 1  0 1 0  0 0 −4/3 1 0 ~ Řešení soustavy L~y = b :   1 0 0 0 | 0  −2 1 0 0 | 5    , ~y = (0, 5, 0, 5)T .  0 1 1 0 | 5  0 0 −4/3 1 | 5 Ze soustavy U~x = ~y obdržíme  −1 1 0  0 3 2   0 0 −3 0 0 0

0 0 3 5

| | | |

 1 0 0 3 2 0  . 0 −3 3  0 0 5

 0 5   , ~x = (1, 1, 1, 1)T . 0  5

Gaussovu eliminační metodu a LU-rozklad lze s výhodou použít i k jiným účel˚ um než jen k řešení soustav lineárních rovnic. Například m˚ užeme vypočítat hodnotu determinantu nebo určit tvar inverzní matice. i) Stanovení hodnoty determinantu: Pomocí GEM převedeme matici A na horní trojúhelníkovou matici. Pak n−1 det A = a011 a122 . . . ann .

39

Pokud jsme provedli LU-rozklad matice A, je det A = detL detU = 1 · u11 u22 . . . unn . 3

Počet násobení při Gaussově eliminaci je řádově n3 . To je při výpočtu determinantů vyšších řádů o hodně méně než n! násobení, která bychom museli realizovat při výpočtu determinantu rozvojem. 

 1 −1 1 Příklad Spočítejte determinant matice A =  2 −1 2  pomocí Gaus3 −1 1 sovy metody a metody LU-rozkladu. Řešení:



  1 −1 1 1    2 −1 2 0 Gaussova metoda: A = −→ 3 −1 1 0   1 0 0 1 −1    2 1 0 0 1 LU-rozklad: A = LU = 3 2 1 0 0

 −1 1 1 0  , det A = −2. 0 −2  1 0  , det A = −2. −2

ii) Obě metody je možné také využít při stanovení inverzní matice: Modifikací GEM je Gaussova-Jordanova metoda. Matici tvořenou p˚ uvodní maticí A a jednotkovou maticí I upravujeme pomocí ekvivalentních úprav tak dlouho, až (A | I) přejde na tvar (I | A−1 ). Pokud chceme určit tvar inverzní matice pomocí LU-rozkladu, využijeme skutečnosti, že když A = LU, pak A−1 = U −1 L−1 . 3

K realizaci Gaussovy-Jordanovy metody potřebujeme řádově n2 +O(n2 ) operací, což je více než při GEM. Proto se Gaussova-Jordanova metoda využívá při praktických výpočtech jen zřídka. Příklad Prostřednictvím Gaussovy-Jordanovy metody   a pomocí LU-rozkladu 1 −1 1 stanovte inverzní matici k matici A =  2 −1 2  . 3 −1 1 40

Řešení: Gaussova-Jordanova metoda:     1 −1 1 | 1 0 0 1 0 0 | −1/2 0 1/2 0  = (I|A−1 ). (A|I) =  2 −1 2 | 0 1 0  −→  0 1 0 | −2 1 3 −1 1 | 0 0 1 0 0 1 | −1/2 1 −1/2 LU-rozklad:     1 −1 1 1 0 0 1 −1 1 0 , A =  2 −1 2  =  2 1 0   0 1 3 2 1 0 0 −2 3 −1 1 



A−1 = U −1 L−1

    1 0 0 1 −1 1/2 −1/2 0 1/2 0  =  −2 1 0  . =  −2 1 0   0 1 1 −2 1 0 0 −1/2 −1/2 0 1/2

Iterační metody Předností iteračních metod je, že nekladou takové nároky na paměť počítače jako metody přímé a že jejich použití je univerzální. Můžeme pomocí nich nejen získat řešení soustavy lineárních rovnic, ale také je můžeme využít ke zpřesnění řešení, které jsme získali již dříve. Základní myšlenka konstrukce iterační metody je následující: Z každé rovnice soustavy A~x = ~b vyjádříme právě jednu neznámou v závislosti na zbývajících neznámých. Soustava A~x = ~b tak přejde na tvar ~x = H~x + ~g . Odtud získáme iterační formuli ~x(k+1) = H~x(k) + ~g . Vlastní iterační proces probíhá tak, že 41

1) zvolíme počáteční iteraci ~x(0) , 2) prostřednictvím iterační formule ~x(k+1) = H~x(k) + ~g určíme další vektory řešení ~x(k+1) , k = 0, 1, 2, . . . , 3) proces ukončíme buď po předem stanoveném počtu iterací, nebo když bude splněna zastavovací podmínka |~x(k+1) − ~x(k) | < δ. (Říkáme pak, že řešení ~x je určeno s přesností δ.) 4) Nakonec odhadneme chybu (k + 1)-ní iterace. Skutečnost, že |~x(k+1) − ~x| < ε, vyjadřuje, že iterace ~x(k+1) aproximuje přesné řešení ~x s chybou ε. Pokud se iterační matice H, popř. vektor ~g , v pr˚ uběhu výpočtu mění, hovoříme o nestacionárním iteračním procesu, v opačném případě o iteračním procesu stacionárním. V dalším textu se omezíme na studium stacionárních iteračních proces˚ u.

a) Jacobiova iterační metoda

Při odvozování iterační metody budeme postupovat tak, že z i-té rovnice soustavy a11 x1 + a12 x2 + · · · + a1n xn = b1 , a21 x1 + a22 x2 + · · · + a2n xn = b2 , ... an1 x1 + an2 x2 + · · · + ann xn = bn 42

(2.5)

vyjádříme i-tou neznámou. Iterační rovnice v (k + 1)-ním kroku, kde k = 0, 1, 2, . . . , mají tvar (k+1)

x1

(k+1)

x2

1 (k) (k) (b1 − a12 x2 − a13 x3 − · · · − a1n x(k) n ), a11 1 (k) (k) = (b2 − a21 x1 − a23 x3 − · · · − a2n x(k) n ), a22 =

... x(k+1) = n Pokud označíme    HJ =  

1 (k) (k) (k) (bn − an1 x1 − an2 x2 − · · · − an,n−1 xn−1 ). ann 0 − aa21 22 .. .

− aa12 11 0 .. .

n1 n2 − aann − aann

. . . − aa1n 11 . . . − aa2n 22 .. ... . ... 0





b1 a11

  .   a gJ =  ..  bn

  ,

ann

dostáváme Jacobiovu iterační formuli ve tvaru ~x(k+1) = HJ ~x(k) + ~gJ .

(2.6)

Nadále v tomto odstavci budeme používat označení D = diag(a11 , a22 , . . . , ann ),     0 0 ... 0 0 0 a12 . . . a1,n−1 a1n  a21 0 . . .  0 0 0  0 . . . a2,n−1 a2n         0 0 ,N =  ... ... ... ... ...  M =  a31 a32 . . . .  ... ... ...   ... ... 0 0 ... 0 an−1,n  an1 an2 . . . an,n−1 0 0 0 ... 0 0 Když matici soustavy A = (aij )ni,j=1 vyjádříme ve tvaru A = M + D + N, pak můžeme psát (M + D + N )~x = ~b, D~x = −(M + N )~x + ~b, ~x = −D−1 (M + N )~x + D−1~b. To znamená, že pro matici HJ a vektor ~gJ ve vztahu (2.6) platí HJ = −D−1 (M + N ), gJ = D−1~b. 43

(2.7)

Definice 2.3 Nechť A je čtvercová matice řádu n a A = M + D + N. Jacobiova metoda pro řešení soustavy A~x = ~b je dána iteračními rovnicemi ~x(k+1) = −D−1 (M + N )~x(k) + D−1~b, k = 0, 1, 2, . . . .

Algoritmus (Jacobiova metoda) (0) Vstup: n, A = (aij )ni,j=1 , ~b = (bi )ni=1 , ~x(0) = (xi )ni=1 , m pro k = 0, 1, . . . , m pro i = 1, 2, . . . , n P Pn (k) (k) (k+1) xi = a1ii (bi − i−1 j=1 aij xj − j=i+1 aij xj ) (m)

(m)

(m)

Výstup: ~x(m) = (x1 , x2 , . . . , xn ) Příklad Jacobiovou metodou řešte soustavu 2x1 + x2 = 2, x1 + 4x2 + x3 = 0, x2 + 2x3 = −2. Řešení: Iterační formule (k+1)

x1

(k+1)

x2

(k+1)

x3 Iterační matice

1 (k) = (2 − x2 ), 2 1 (k) (k) = (−x1 − x3 ), 4 1 (k) = (−2 − x2 ). 2



 0 − 12 0 HJ =  − 14 0 − 14  , 0 − 12 0 ~gJ = (1, 0, −1)T . 44

(2.8)

Když zvolíme ~x(0) = (1, 1, 1), pak ~x(1) = (0, 500, −0, 500, −1, 50), ~x(2) = (1, 25, 0, 250, −0, 750), ~x(3) = (0, 875, −0, 125, −1, 13), ~x(4) = (1, 07, 0, 0638, −0, 940), ~x(5) = (0, 970, −0, 0325, −1, 03), ~x(6) = (1, 02, 0, 0150, −0, 985), ~x(7) = (0, 995, −0, 00875, −1, 01), ~x(8) = (1, 01, 0, 00375, −0, 995), ~x(9) = (1, 00, −0, 00375, −1, 00), ~x(10) = (1, 00, 0, −1, 00), ~x(11) = (1, 0, −1). Po deseti iteracích jsme dospěli k přesnému řešení ~x = (1, 0, −1). Příklad Jacobiovou metodou řešte soustavu x1 + 3x2 = 4, −3x1 + x2 = −2. Řešení: Iterační formule (k+1)

(k)

x1 = 4 − 3x2 , (k+1) (k) x2 = −2 + 3x1 . Iterační matice

µ HJ =

0 −3 3 0

¶ , ~gJ = (4, −2)T .

Když zvolíme ~x(0) = (0, 0), pak ~x(1) = (4, −2), ~x(2) = (10, 10), ~x(3) = (−26, 28), ~x(4) = (80, −80), ~x(5) = (−236, −242). 45

Vidíme, že v tomto případě Jacobiova metoda nekonverguje. K otázce konvergence se vrátíme v průběhu následujícího výkladu (viz strana 48).

b) Gaussova-Seidelova metoda Gaussova-Seidelova metoda se od Jacobiovy liší tím, že všechny vypočtené hodnoty okamžitě používáme v dalším iteračním kroku. To znamená, že v (k + 1)-ním kroku (k = 0, 1, 2, . . . ) obdržíme iterační formule 1 (k) (k) (b1 − a12 x2 − a13 x3 − · · · − a1n x(k) n ), a11 1 (k+1) (k+1) (k) x2 = (b2 − a21 x1 − a23 x3 − · · · − a2n x(k) n ), a22 ............... 1 (k+1) (k+1) (k+1) − · · · − an,n−1 xn−1 ). − an2 x2 x(k+1) = (bn − an1 x1 n ann (k+1)

x1

=

Použijme nyní stejného značení jako na str. 43. Když uvážíme, že řešíme soustavu (M + D + N )~x = ~b, pak (M + D)~x = −N~x + ~b. ~x = −(M + D)−1 N~x + (M + D)−1~b.

(2.9)

Dostáváme novou iterační metodu ~x(k+1) = HGS ~x(k) + ~gGS , kde

HGS = −(M + D)−1 N, gGS = (M + D)−1~b.

(2.10)

Definice 2.4 (Gauss 1822) Když A je čtvercová matice řádu n a platí A = M + D + N, pak Gaussova-Seidelova metoda řešení soutavy lineárních rovnic A~x = ~b je určena iteračními rovnicemi ~x(k+1) = −(M + D)−1 N~x(k) + (M + D)−1~b, k = 0, 1, 2, . . . . 46

(2.11)

Algoritmus (Gaussova-Seidelova metoda) (0) Vstup: n, A = (aij )ni,j=1 , ~b = (bi )ni=1 , ~x(0) = (xi )ni=1 , m pro k = 0, 1, . . . , m pro i = 1, 2, . . . , n P P (k+1) (k+1) (k) xi = a1ii (bi − i−1 − nj=i+1 aij xj ) j=1 aij xj (m)

(m)

(m)

Výstup: ~x(m) = (x1 , x2 , . . . , xn ) Předností tohoto algoritmu je, že klade menší nároky na paměť počítače než Jacobiova metoda. Příklad Gaussovou-Seidelovou metodou řešte soustavu 2x1 + x2 = 2, x1 + 4x2 + x3 = 0, x2 + 2x3 = −2. Řešení: Iterační formule (k+1)

x1

(k+1)

x2

(k+1)

x3

1 (k) = (2 − x2 ), 2 (k) 1 1 x (k+1) (k) (k) = (−x1 − x3 ) = (−1 + 2 − x3 ), 4 4 2 1 1 7 1 (k) 1 (k) (k+1) = (−2 − x2 ) = (− − x2 + x3 ). 2 2 4 8 4

Iterační matice a vektor pravých stran pak jsou tvaru 

HGS

 0 0 − 12 1 7 =  0 − 18 − 14  , ~gGS = (1, − , − )T . 4 8 1 1 0 − 16 8 47

Když zvolíme ~x(0) = (1, 1, 1), pak ~x(1) = (0, 500, −0, 375, −0, 815), ~x(2) = (1, 19, −0, 0938, −0, 955), ~x(3) = (1, 05, −0, 0238, −0, 990), ~x(4) = (1, 01, −0, 00500, −1, 00), ~x(5) = (1, 01, −0, 00250, −1, 00), ~x(6) = (1, 00, 0, 00, −1, 00), ~x(7) = (1, 0, −1). V tomto případě jsme k přesnému řešení ~x = (1, 0, −1) dorazili už po šesti iteracích, tedy rychleji než při použití Jacobiovy iterační metody. Oběma výše uvedenými metodami jsme získali posloupnost iterací {~x(k) }. Viděli jsme však také, že tato posloupnost nemusí vždy konvergovat. Jak je to tedy s konvergencí iteračního procesu? Platí následující dvě tvrzení: Věta 2.2 (Nutná a postačující podmínka konvergence) Nechť λi (H) jsou vlastní čísla matice H. Iterace ~x (k+1) = H~x (k) + ~g konvergují pro libovolné ~g a libovolnou počáteční hodnotu ~x(0) , právě když spektrální poloměr iterační matice ρ(H) = max{λi (H)} < 1. i

Důkaz. Iterační matice H je podobná jisté Jordanově matici J. To znamená, že obě matice mají stejná vlastní čísla (viz věta 1.1). Dále existuje regulární matice T tak, že H = T JT −1 a H k = T J k T −1 . Když O je nulová matice a podle předpokladů věty 1 > ρ(H) = λ1 (H) > λ2 (H) > · · · > λn (H), pak lim J k = O, a tedy lim H k = O. k→∞

k→∞

Obráceně, když lim H k = O, je také lim J k = O, a tedy λi (H) < 1 pro všechna i.

k→∞

k→∞

48

Věta 2.3 (Postačující podmínka konvergence) Posloupnost iterací ~x(k+1) = H~x(k) + ~g konverguje pro libovolné ~x(0) a libovolné ~g , když norma iterační matice splňuje nerovnost ||H|| ≤ q < 1.

Důkaz. Protože ~x = H~x + ~g a ~x(k+1) = H~x(k) + ~g , máme pro chybu řešení k~x −~x(k+1) k ≤ kHk k~x −~x(k) k ≤ kH 2 k k~x −~x(k−1) k ≤ · · · ≤ kH k+1 k k~x −~x(0) k. Protože kH k+1 k ≤ kHkk+1 a ||H|| ≤ q < 1, je limk→∞ kH k+1 k = 0. To znamená, že limk→∞ k~x − ~x(k+1) k = 0 a metoda konverguje. Příklad Vhodnou iterační metodou řešte soustavu x1 + 2x2 − 2x3 = 1, x1 + x2 + x3 = 3, 2x1 + 2x2 + x3 = 5. Řešení: O tom, zda zvolíme Jacobiovu nebo Gaussovu-Seidelovu metodu, se rozhodneme podle toho, která metoda bude konvergovat. K tomu využijeme větu 2.2. Matice soustavy   1 2 −2 A =  1 1 1 . 2 2 1 Iterační matice Jacobiovy metody 

0 −2 2  −1 0 −1 HJ = −2 −2 0

 .

Protože 



−λ −2 2  −1 −λ −1 |HJ − λI| = −2 −2 −λ

 = λ3 = 0 ⇒ λ1,2,3 = 0 < 1, 49

bude Jacobiova metoda podle věty 2.2 konvergovat. Iterační matice Gaussovy-Seidelovy metody     1 0 0 0 2 −2 HGS = −(M +D)−1 N =  −1 1 0   0 0 1  =  0 −2 1 0 0 0   −λ −2 2  = (2−λ)2 (−λ) = 0 ⇒ λ1 |HGS −λI| =  0 2 − λ −3 0 0 2−λ

 0 −2 2 0 2 −3  . 0 0 2 = 0, λ2,3 = 2,

tzn. spektrální poloměr ρ(HGS ) = 2 > 1. Gaussova-Seidelova metoda nemůže podle věty 2.2 konvergovat. Úlohu budeme řešit Jacobiovou metodou. Iterační formule mají tvar (k+1)

= 1 − 2x2 + 2x3 ,

(k+1)

= 3 − x1 − x3 ,

(k+1)

= 5 − 2x1 − 2x2 .

x1 x2 x3

(k)

(k)

(k)

(k)

(k)

(k)

Když zvolíme ~x(0) = (0, 0, 0), pak ~x(1) = (1, 3, 5), ~x(2) = (5, −3, −3), ~x(3) = (1, 1, 1), ~x(4) = (1, 1, 1). Vektor ~x = (1, 1, 1) je řešením zadané soustavy. Výše uvedené konvergenční podmínky se však v praxi špatně ověřují. V některých specifických případech je možné rozhodnout o konvergenci už z toho, jaké vlastnosti má matice soustavy A. Tak například platí: Věta 2.4 Jacobiova metoda konverguje pro libovolné ~g a ~x 0 , když q∞ = max i

n X aij | | < 1, aii j=1, j6=i

50

nebo

n X aij qS = max | | < 1, j aii i=1, i6=j

nebo

1/2



n X aij   qE =  | |2  aii  i,j=1,

< 1.

i6=j

Pak pro apriorní odhad chyby dostáváme ||~x(n) − ~x(0) ||i ≤

qin ||~x(1) − ~x(0) ||i , i = ∞, S, E. 1 − qi

Pro aposteriorní odhad chyby je ||~x(n) − ~x(0) ||i ≤

qi ||~x(n) − ~x(n−1) ||i , i = ∞, S, E. 1 − qi

Důkaz. Iterační matice HJ = −D−1 (M +N ) (viz str. 43) má na hlavní diagoa nále nulové prvky, jinak jsou rovny − ai,j . Čísla q∞ , qS , qE představují normy ii Jacobiho iterační matice. Pokud jsou menší než 1, jsou splněny předpoklady věty 1.5. Odtud tvrzení. Poznámka Z podmínky qS < 1 plyne n X

|aij | < |aii |, j = 1, . . . , n.

i=1, i6=j

To znamená, že A je ostře diagonálně dominantní. Jacobiova metoda pak konverguje pro libovolnou počáteční aproximaci ~x(0) . Totéž platí i pro Gaussovu-Seidelovu metodu. Věta 2.5 Pokud matice A je symetrická pozitivně definitní, pak GaussovaSeidelova metoda konverguje pro libovolné ~x(0) . 51

Věta 2.6 (Sassenfeldovo kriterium) Gaussova-Seidelova metoda konverguje pro libovolné ~g a libovolné ~x(0) , když p = max pi < 1, i=1,...,n

kde p1 =

i−1 n n X X X aij aij a1j | |pj + | | < 1, i = 2, . . . , n. | | < 1, pi = a11 aii aii j=2 j=1 j=i+1

Pro apriorní odhad chyby pak platí ||~x(k) − ~x(0) ||∞ ≤

qn ||~x(1) − ~x(0) ||∞ . 1−q

a pro aposteriorní odhad chyby je ||~x(k) − ~x(0) ||∞ ≤

q ||~x(k) − ~x(k−1) ||∞ . 1−q

Důkaz. Uvažujme soustavu (M + D)~x = −N~z, kde ||z||∞ = 1, pak xi = −

n X aij xj − zj , i = 1, . . . , n. aii a j=i+1 ii

i−1 X aij j=1

Indukcí obdržíme |xi | ≤ pi , i = 1, . . . , n, a tedy ||~x||∞ ≤ p. Odtud ||(M + D)−1 N ||∞ ≤ p. Tvrzení plyne z věty 1.5. Příklad Rozhodněte, zda Gaussova-Seidelova metoda bude konvergovat v případě, že matice soustavy má tvar   2 −1 0 0  −1 2 −1 0   A=  0 −1 2 −1  . 0 0 −1 2 52

Řešení: Matice A je sice symetrická, ale není ostře diagonálně dominantní. Nemůžeme tedy použít větu 2.5. Protože však 1 1 1 1 1 3 1 1 7 1 p1 = < 1, p2 = · + < 1, p3 = · + , p4 = · + < 1, 2 2 2 2 2 4 2 2 8 2 jsou splněny předpoklady věty 2.6, a tedy Gaussova-Seidelova metoda bude konvergovat. V případě, že spektrální poloměr ρ(H) je číslo blízké 1, mohou iterace konvergovat velmi pomalu. Východiskem z tohoto problému je použití numerické metody, která proces urychlí. Jacobiovu iterační formuli (2.8) přepišme na tvar ~x(k+1) = ~x(k) − D−1 (M + D + N )~x(k) + D−1~b. To znamená ~x(k+1) = ~x(k) + ~r(k+1) , kde ~r(k+1) = D−1 (~b − A~x(k) ). Když reziduum ~r(k+1) vynásobíme vhodným číslem ω, dostaneme novou numerickou metodu, která může konvergovat rychleji. Definice 2.5 Iterační metoda ~x(k+1) = ~x(k) + ωD−1 (~b − A~x(k) ), k = 0, 1, . . . ,

(2.12)

představuje Jacobiovu metodu s relaxací. Číslo ω je tzv. relaxační parametr. Poznámka Vztah (2.12) lze psát ve složkách i−1

(k+1) xi

=

(k) xi

n

X X ω (k) (k) + (bi − aij xj − aij xj ), i = 1, . . . , n. aii j=1 j=i

Věta 2.7 Nechť HJ je iterační matice Jacobiovy metody, λ1 > · · · > λn jsou vlastní čísla matice HJ , ρ(HJ ) < 1. Spektrální poloměr iterační matice HJω Jacobiovy relaxační metody bude minimální, když zvolíme relaxační parametr 2 . ω= 2 − λ1 − λn 53

Důkaz. Pro ω > 0 jsou rovnice HJ ~v = λ~v [(1 − ω)I + ωHJ ]~v = [1 − ω + ωλ]~v ekvivalentní. Proto také vlastní číslo λ matice HJ odpovídá vlastnímu číslu 1 − ω + ωλ matice (1 − ω)I + ωHJ = HJω . Vlastní čísla matice HJω jsou reálná, 1 − ω + ωλn je nejmenší a 1 − ω + ωλ1 největší vlastní číslo. Spektrální poloměr matice HJω bude nejmenší, když nejmenší a největší vlastní číslo budou mít opačnou hodnotu. Tzn. 1 − ω + ωλn = −(1 − ω + ωλ1 ). Odtud tvrzení. K urychlení iteračního procesu bývá nejčastěji využívána tzv. superrelaxační metoda. Obdobně jako v předchozím případě přepíšeme Gaussovu-Seidelovu iterační formuli (2.11) na tvar i−1

(k+1) xi

=

(k) xi

+

(k+1) ri ,

(k+1)

a i-té reziduum ri iterační formuli (k+1) xi

=

(k) xi

+

kde

(k+1) ri

n

X X 1 (k+1) (k) = (bi − aij xj − aij xj ) aii j=1 j=i

vynásobíme vhodným parametrem ω. Obdržíme novou

(k+1) ωri

= (1 −

(k) ω)xi

i−1 n X X ω (k+1) (k) + (bi − aij xj − aij xj ). aii j=1 j=i+1

Definice 2.6 Iterační metoda ~x(k+1) = ~x(k) + ωD−1 (~b − M~x(k+1) − (D + N )~x(k) ), n = 0, 1, . . . ,

(2.13)

představuje Gaussovu-Seidelovu relaxační metodu. Poznámka Stručně Gaussovu-Seidelovu relaxační metodu budeme označovat jako SOR (successive overrelaxaction). Věta 2.8 K tomu, aby Gaussova-Seidelova relaxační metoda (2.13) konvergovala, je nutné, aby 0 < ω < 2. 54

Důkaz. Vlastní čísla iterační matice HGSω jsou kořeny charakteristického polynomu a platí pro ně n Y

λj = det HGSω = det {(D + ωM )−1 [(1 − ω)D − ωN ]}.

j=1

Protože determinant součinu je roven součinu determinantů a matice (D + ωM )−1 a (1 − ω)D − ωN jsou trojúhelníkové, dostáváme n Y

λj = det (D + ωM )−1 det [(1 − ω)D − ωN ] = (1 − ω)n .

j=1

Odtud λ1 > ω. Z věty 2.2 dostaneme, že SOR bude konvergovat, když 1 > |ρ(Hω )| ≥ |1 − ω|. Odtud tvrzení. Algoritmus (Superrelaxační metoda - SOR) Vstup: n, A = (aij )ni,j=1 , ~b = (bi )ni=1 , ~x (0) = (x1 , . . . , xn )T , ω, m pro k = 0, 1, . . . , n pro i = 1, 2, . . . , m P P (k+1) (k) (k+1) (k) − nj=i+1 aij xj ) xi = (1 − ω)xi + aωii (bi − i−1 j=1 aij xj Výstup: ~x (k+1)

Podmíněnost soustav lineárních rovnic V úvodní kapitole jsme hovořili o tom, že každé číslo se nemusí v počítači zobrazit přesně. Pokud úloha je špatně podmíněná, mohou zaokrouhlovací chyby naprosto znehodnotit celý výpočet. Podívejme se proto nyní, jak je tomu v případě soustav lineárních rovnic. Označme A, ~b — přesné hodnoty matice soustavy a vektoru pravých stran, 55

~x — k nim příslušný vektor řešení. ˜ ~˜b — matice soustavy a vektoru pravých stran, které jsou zatížené chybou, A, ~x˜ — k nim spočtený vektor řešení. ˜ ∆~b = ~b − ~˜b, ∆~x = ~x − ~x˜. Označme ∆A = A − A, Nejprve uvažujme případ, kdy matice A není zadána přesně a vektor ~b ano. To znamená ∆A 6= 0, ∆~b = 0 a řešíme soustavu (A + ∆A)(~x + ∆~x) = ~b. (A + ∆A)∆~x = ~b − A~x − ∆A~x. Protože A~x = ~b, dostaneme A∆~x = −∆A~x − ∆A∆~x, ∆~x = −A−1 ∆A(~x + ∆~x). Pak ale

1

||∆~x|| ||∆A|| ≤ ||A−1 || ||A|| . ||~x + ∆~x|| ||A|| To znamená, že číslo podmíněnosti (viz definice 1.5) cp = ||A−1 || ||A||. Když je matice A zadána přesně a vektor ~b ne, tj. ∆A = 0, ∆~b 6= 0, pak A(~x + ∆~x) = ~b + ∆~b, A∆~x = ~b − A~x + ∆~b, ∆~x = A−1 ∆~b, ||∆~x|| ≤ ||A−1 || ||∆~b||. Zároveň z (1.14) plyne k~bk = kA~xk ≤ kAk k~xk, 1

Definice různých typů norem viz strana 18. Protože normy (1.8) - (1.10) jsou ekvivalentní, nezáleží na tom, se kterým typem normy (nadále pevně zvoleným) se rozhodneme pracovat.

56

a tedy ||~x|| ≥ Dostáváme proto

||~b|| . ||A||

||∆~x|| ||∆~b|| . ≤ ||A−1 || ||A|| ||~x|| ||~b||

I zde je číslo podmíněnosti cp = ||A−1 || ||A||.

(2.14)

Dá se ukázat, že totéž platí i v případě, kdy ∆A 6= 0 a ∆~b 6= 0. Můžeme tedy uzavřít tím, že pro řešení soustavy lineárních rovnic je číslo podmíněnosti definováno vztahem (2.14).

57

58

Kapitola 3 Vlastní čísla a vlastní vektory matice

V předchozím odstavci jsme viděli, jak je při studiu konvergence iteračních metod důležité umět určit vlastní čísla a vlastní vektory regulární matice. Náš cíl však může být různý. • Někdy například při zkoumání konvergence iteračních metod potřebujeme najít pouze největší (tzv. dominantní) vlastní číslo. V takovém případě řešíme tzv. částečný problém vlastních čísel. • Jindy nám p˚ ujde o to, abychom určili všechna vlastní čísla dané matice. Pak říkáme, že řešíme tzv. úplný problém vlastních čísel. Než však začneme vlastní čísla hledat, je dobré mít představu o tom, jaká je jejich přibližná hodnota. Nesmíme ovšem zapomenout, že i když prvky matice A jsou reálná čísla, mohou být vlastní čísla matice A komplexní. K odhadu polohy vlastních čísel lze použít některou z následujících vět: Věta 3.1 (Geršgorin) Všechna vlastní čísla matice A = (aij )ni,j=1 leží v komplexní rovině ve sjednocení kruhů |z − aii | ≤

n X i=1, i6=j

59

|aij |.

Příklad Odhadněte polohu vlastních čísla matice  4  1 A= 0

čísel a hodnotu nejmenšího vlastního  1 0 2 1 . 1 1

Řešení: Všechna vlastní čísla leží ve sjednocení kruhů (v komplexní rovině): |z − 4| ≤ 1, |z − 2| ≤ 2, |z − 1| ≤ 1 → Reλ ∈ h0, 5i.

Im z

1

2

4

Re z

Věta 3.2 Když je λn nejmenší vlastní číslo matice A = (aij )ni,j=1 , pak platí odhad n X √ n |aij |2 )1/2 . |λn | ≤ r1 r2 . . . rn , kde ri = ( j=1



 1 1 0 Příklad Odhadněte polohu vlastních čísel matice  0 2 0  . 0 0 3 √ √ √ √ Řešení: r1 = 12 + 12 = 2, r2 = 22p= 2, r3 = 32 = 3. Pro nejmenší √ . 3 vlastní číslo zadané matice platí |λn | ≤ 2 · 3 · 2 = 2, 0396. 60

Částečný problém vlastních čísel Jak bylo podotknuto výše, když řešíme částečný problém vlastních čísel, chceme najít některé, většinou dominantní, vlastní číslo. Uvedeme si dva zp˚ usoby řešení.

a) Mocninná metoda Věta 3.3 (Mocninná metoda) Nechť A je čtvercová matice řádu n, |λ1 | > |λ2 | ≥ · · · ≥ |λn | její vlastní čísla a ~v1 , ~v2 , . . . , ~vn jsou k nim příslušné lineárně nezávislé vlastní vektory. Položme ~y

(0)

=

n X

ai~vi ,

i=1

~y (k) = A~y (k−1) , k = 1, . . . , n. Když existují limity (k+1)

λ1,j = lim

yj

k→∞

, j = 1, . . . , n,

(k)

yj

pak pro největší vlastní číslo matice A platí n

1X λ1 = λ1,j . n j=1

Důkaz. Položme ~y

(0)

=

n X i=1

61

ai~vi

a pro k ≥ 1 ~y (k) = A~y (k−1) = A2 ~y (k−2) = · · · = Ak ~y (0) . To znamená, že ~y (k) = Ak (a1~v1 + a2~v2 + · · · + an~vn ) = a1 Ak~v1 + · · · + an Ak~vn . Protože λi jsou vlastní čísla matice A a ~vi k nim příslušné vlastní vektory, je A~vi = λi~vi . Po dosazení do předchozího vztahu dostaneme, že ~y (k) = a1 λk1 ~v1 + · · · + an λkn~vn = λk1 (a1~v1 + a2 ( Když označíme ~ε k =

Pn i=1

λ2 k λn ) ~v2 · · · + an ( )k~vn ). λ1 λ1

ai ( λλ1i )k~vi , máme pro j-tou složku vektoru ~y k (k)

yj

= λk1 (a1 v1,j + εkj ).

(3.1)

Protože λ1 je podle předpokladu největší vlastní číslo, platí pro všechna i ≥ 2, že | λλ1i | < 1, a tedy lim εkj = 0. M˚ užeme psát k→∞

(k+1)

λ1,j = lim

yj

k→∞

(k)

yj

k+1 λk+1 ) 1 (a1 v1j + εj = lim , j = 1, 2, . . . , n. k→∞ λk1 (a1 v1j + εkj )

Největší vlastní číslo λ1 ≈

1 n

Pn j=1

λ1,j .

Algoritmus (Mocninná metoda) Vstup: n, A = aij , ~y (0) , ε λ01 = 0 (k) (k−1) pro k = 1, 2, . . . , dokud |λ1 − λ1 |>ε pro i = 1,P 2, . . . , n (k−1) (k) yi = nj=1 aij yj (k)

λ1,i = (k)

λ1 =

1 n

(k) yi (k−1) yi

Pn

i=1

λ1,i

(k)

Výstup: λ1

62



 1 2 0 Příklad Určete hodnotu dominantního vlastního čísla matice  0 0 1  . 3 6 2 Řešení: Položíme ~y (0) = (1, 1, 1)T a mocninnou metodou spočteme ~y (1) = (3, 1, 11)T , ~y (2) = (5, 11, 37)T , ~y (3) = (27, 37, 155)T , ~y (4) = (101, 155, 613)T , ~y (5) = (411, 613, 2459)T , ~y (6) = (1637, 2459, 9829)T , ~y (7) = (6555, 9829, 39323)T , ~y (8) = (26213, 39323, 157285)T , ~y (9) = (104859, 157285, 629147)T , ~y (10) = (419429, 629147, 2516581)T ,

(1)

λ1 (2) λ1 (3) λ1 (4) λ1 (5) λ1 (6) λ1 (7) λ1 (8) λ1 (9) λ1

= 5.3434, = 4.3176, = 3.9616, = 4.0119, = 3.9972, = 4.0007, = 3.9998, = 4.0000, = 4.0000.

Dominantní vlastní číslo λ1 = 4. Nevýhodou metody je, že předem nelze rozhodnout o tom, zda jsou splněny požadované předpoklady, tj. zda |λ1 | > |λ2 | ≥ · · · ≥ |λn |. Také se těžko odhaduje chyba aproximace.

b) Metoda Rayleighova podílu Modifikací mocninné metody je metoda Rayleighova podílu, kterou lze aplikovat tehdy, když jsou splněny stejné předpoklady jako v mocninné metodě, a matice A je navíc symetrická. Věta 3.4 (Rayleigh) Nechť A je symetrická matice řádu n, pro její vlastní čísla platí |λ1 | > |λ2 | ≥ · · · ≥ |λn | 63

a ~v1 , ~v2 , . . . , ~vn jsou odpovídající ortonormální vlastní vektory. Když označíme ~y (0) = ~y

(k)

n X

ai~vi ,

i=1 k (0)

= A ~y

, k = 1, . . . , n,

pak (~y (k) )T ~y (k+1) . k→∞ (~ y (k) )T ~y (k)

λ1 = lim

Důkaz. Použijeme stejných označení jako v důkazu věty 3.3. Především podle (k) (3.1) položíme yj = λk1 (a1 v1,j + εkj ). Když A je symetrická, pak vlastní vektory jsou ortonormální a vektory ~y (k) splňují 2 (~y (k) )T ~y (k) = λ2k v1T + (~ε k )T )(a1~v1 + ~ε k ) = λ2k ε k )T ~ε k ), 1 (a1~ 1 (a1 + (~

(~y (k) )T A~y (k) = λk1 (a1~v1T +(~ε k )T )λk+1 v1 +~ε k+1 ) = (λk1 )2k+1 (a1 2 +(~ε k )T ~ε k+1 ). 1 (a1~ Pro |λ1 | > |λ2 | ≥ · · · ≥ |λn | je limk→∞ (~ε k )T ~ε k = limk→∞ (~ε k )T ~ε k+1 = 0, a tedy (~y (k) )T A~y (k) (~y (k) )T ~y (k+1) λ1 = lim = lim . k→∞ (~ k→∞ (~ y (k) )T ~y (k) y (k) )T ~y (k)

Algoritmus (Metoda Rayleighova podílu) Vstup: n, A = (aij ), ~y (0) , ε (0)

λ1 = 0 (k−1) (k) pro k = 1, 2, . . . , dokud |λ1 − λ1 |>ε (k+1)

λ1

=

(~ y k )T A~ yk k T (~ y ) ~ yk

(k+1)

Výstup: λ1



 4 1 0 Příklad Najděte dominantní vlastní číslo matice A =  1 2 1  . 0 1 1 64

Řešení: Položíme ~y (0) = (1, 1, 1)T . Metodou Rayleighova podílu spočteme (1)

(~ y (0) )T ~ y (1) (0) T (~ y ) ~ y (0) (~ y (1) )T ~ y (2) (~ y (1) )T ~ y (1) (~ y (2) )T ~ y (3) (~ y (2) )T ~ y (2)

~y (1) = A~y (0) = (5, 4, 2)T ,

λ1 =

~y (2) = A~y (1) = (24, 15, 6)T ,

λ1 =

~y (3) = A~y (2) = (111, 60, 21)T ,

λ1 =

~y (4) = A~y (3) = (504, 252, 81)T , ~y (5) = A~y (4) = (2268, 1089, 333)T , ~y (6) = A~y (5) = (10161, 4779, 1422)T , ~y (7) = A~y (6) = (45423, 21141, 6201)T , ~y (8) = A~y (7) = (202833, 93906, 27342)T , ~y (9) = A~y (8) = (905238, 417987, 121248)T , ~y (10) = A~y (9) = (4038939, 1862460, 539235)T ,

(4) λ1 = 4, 4472, (5) λ1 = 4, 4571, (6) λ1 = 4, 4597, (7) λ1 = 4, 4603, (8) λ1 = 4, 4605, (9) λ1 = 4, 4605, (10) λ1 = 4, 4605.

(2)

(3)

= 3, 6667, = 4, 2667, = 4, 4086,

Úplný problém vlastních čísel Když chceme určit všechna vlastní čísla dané matice, m˚ užeme zvolit jeden ze dvou přístup˚ u. Buď vlastní čísla určíme přímo jako kořeny charakteristického polynomu p(λ) (definice charakteristického polynomu viz strana 15), nebo využijeme skutečnosti, že dvě podobné matice mají stejná vlastní čísla.

a) Přímý výpočet vlastních čísel V prvním případě, pokud chceme vlastní čísla počítat přímo jako kořeny charakteristického polynomu a řád matice je vysoký, znamená to počítat také determinant vysokého řádu. To je problém náročný a častokrát nerealizovatelný. Situace se však zjednoduší v případě, když matice A je třídiagonální.

65

Věta 3.5 Když A je třídiagonální  a1 c1  b2 a 2 A= 0 0

matice tvaru 0 ... c2 . . . 0 ...

 0 0 0 0  b n an

a položíme f−1 (λ) = 0, f0 (λ) = 1, fk (λ) = (ak − λ)fk−1 (λ) − bk ck−1 fk−2 (λ) pro k = 1, 2, . . . , n, pak charakteristický polynom matice A je p(λ) = fn (λ).

Poznámka Vlastní čísla matice A jsou kořeny charakteristického polynomu p(λ). Algoritmus (Charakteristický polynom třídiagonální matice) Vstup: n, A, f−1 (λ) = 0, f0 (λ) = 1 pro k = 1, 2, . . . , n fk (λ) = (ak − λ)fk−1 (λ) − bk ck−1 fk−2 (λ) Výstup: p(λ) = fn (λ) Poznámka Tento postup je možné využít i pro obecnou čtvercovou matici po té, co původní matici převedeme na matici třídiagonální, která má stejná vlastní čísla. (Viz [8] strana 542.) 

 2 −1 0 Příklad Určete charakteristický polynom matice  −1 2 −1  . 0 −1 1 66

Řešení: Položíme f−1 (λ) = 0, f0 (λ) = 1 a spočteme f1 = 2 − λ, f2 = (2 − λ)2 − 1, f3 = ((2 − λ)2 − 1)(1 − λ) − (2 − λ) = −λ3 + 5λ2 − 6λ + 1 = P (λ).

b) Určení vlastních čísel metodou LU-rozkladu Tato metoda patří mezi ty, které využívají podobnosti matic. Podle věty 1.1 podobné matice mají stejná vlastní čísla. Věta 3.6 Nechť A je reálná regulární matice a A = LU je její LU-rozklad, pak matice B = U L je podobná matici A. Důkaz. Když A = LU, je U = L−1 A. Odtud B = U L = L−1 AL. To znamená, že matice A a B jsou podobné.

Když budeme chtít spočítat vlastní čísla matice A metodou LU-rozkladu, položíme nejprve A0 = A a provedeme LU-rozklad této matice: A0 = L0 U0 . Záměnou pořadí násobení obou matic dostaneme novou matici A1 = U0 L0 . Provedeme LU-rozklad matice A1 : A1 = L1 U1 , a označíme A2 = U1 L1 atd. Dostaneme posloupnost matic Ak , k = 0, 1, . . . Protože Ak = Uk−1 Lk−1 a zároveň ze vztahu Ak−1 = Lk−1 Uk−1 plyne Uk−1 = L−1 k−1 Ak−1 , dostáváme −1 Ak = Uk−1 Lk−1 = L−1 k−1 Ak−1 Lk−1 = Lk−1 Uk−2 Lk−2 Lk−1 = . . . −1 −1 −1 = L−1 k−1 Lk−2 Ak−2 Lk−2 Lk−1 = · · · = Lk−1 . . . L0 A0 L0 . . . Lk−1 .

Označme Mk = L0 L1 . . . Lk . Platí následující věta Věta 3.7 Nechť A je regulární matice. Pokud posloupnost matic Mk konverguje k regulární matici L, konvergují také matice Ak k horní trojúhelníkové matici s vlastními čísly na hlavní diagonále. 67

Důkaz. Podle předpokladu posloupnost matic Mk konverguje k L. Dále z toho, −1 že Mk = L0 L1 . . . Lk , dostaneme Lk = Mk−1 Mk , a tedy −1 lim Lk = lim Mk−1 Mk = L−1 L = I = lim L−1 k .

k→∞

k→∞

k→∞

Zároveň existuje −1 −1 −1 lim Uk−1 = lim Ak L−1 AM = U. k−1 = lim Mk−1 A0 Mk−1 Lk−1 = M

k→∞

k→∞

k→∞

Pak lim Ak = lim Lk Uk = lim Uk = U.

k→∞

k→∞

k→∞

Algoritmus (Určení vlastních čísel pomocí LU-rozkladu) Vstup: n, A = (ai,j ), m A0 := A, A0 = L0 U0 pro k = 1, 2, . . . , m Ak := Uk−1 Lk−1 LU-rozklad: Ak = Lk Uk m Výstup: Am = (am ij ), λi ≈ aii

Poznámka V některých speciálních případech lze využít vlastností zadané matice a výpočet zjednodušit. Například když A je třídiagonální symetrická pozitivně definitní matice, použijeme ke stanovení počátečního LUrozkladu Choleského algoritmus. Označme Lk = (aij )i,j , Lk+1 = (bij )i,j , Ak+1 = Lk+1 LTk+1 . Prvky matice Ak+1 = B spočteme následovně: Algoritmus (Určení vlastních čísel metodou LU-rozkladu pro symetrickou třídiagonální matici) Vstup: n, aii , ai+1,i , i = 1, 2, . . . n pro i = 1, 2, . . . , n bii = a2ii + a2i+1,i − b2i−1,i b2i+1,i = b12 a2i+1,i a2i+1,i+1 ii

Výstup: bi,i , bi+1,i , i = 1, 2, . . . , n 68

c) Určení vlastních čísel metodou ortogonálních transformací Také zde využíváme podobnosti matic. Platí následující věta: Věta 3.8 Když A je matice řádu n a Q je ortogonální matice stejného řádu, pak matice A a QT AQ jsou podobné. Věta 3.9 (Jacobi 1846) Když A je reálná symetrická matice řádu n a p-tý sloupec 

Qpq

           =           

1

q-tý sloupec

 .. .. . .  .. .. ..  . . .   .. ..  1 . .  . . . . . . cos α . . . . . . . . . − sin α . . . . . . . . . . . .   .. ..   . 1 .  .. .. ...  . .   .. ..  . 1 .  . . . . . . sin α . . . . . . . . . cos α . . . . . . . . . . . .    .. ..  . . 1  .. .. ..  . . .  .. .. . . 1

je ortogonální matice stejného řádu, pak B = QTpq AQpq je reálná symetrická matice, která je tvořena prvky bqq bpp bpq biq bip bil

= app sin2 α − apq sin 2α + aqq cos2 α, = app cos2 α + apq sin 2α + aqq sin2 α, = bqp = apq cos 2α + 21 (aqq − app sin 2α), = bqi = −aip sin α + aiq cos α pro i 6= p, q, = bpi = aip cos α + aiq cos α pro i 6= p, q, = ail , i, l 6= p, q.

69

Poznámka α = π4 .

Prvky bpq = bqp = 0, když tg 2α =

2apq , app −aqq

app 6= aqq , tj. když

Věty 3.8 a 3.9 jsou základem Jacobiovy metody ortogonálních transformací. Tato metoda umožňuje postupně pomocí posloupnosti matic Qpq převést matici A na diagonální s vlastními čísly na hlavní diagonále. Výpočet probíhá následovně: V každém kroku výpočtu zvolíme ortogonální transformaci Qpq tak, aby se v nové matici B vynuloval nediagonální prvek s největší absolutní hodnotou. To znamená, že ve stávající matici A vybereme prvek v tom řádku p a sloupci q, pro který apq = maxi,j {|aij |}. Hodnoty prvk˚ u nově vzniklé matice B spočteme ze vztahů app + aqq + r , 2 app aqq − a2pq , = bpp = bqp = 0, = aip c + aiq s = bpi , i = 1, 2 . . . , n, i 6= p, q, = −aip s + aiq c = bqi , i = 1, 2 . . . , n, i 6= p, q = aij jinak.

bpp = bqq bpq bip biq bij

V předchozích rovnicích jsme použili označení r2 = (app + aqq )2 + 4apq , 1 app − aqq s2 = − , 2 2r 1 app − aqq c2 = + . 2 2r

Vypočtené prvky bij matice B přeznačíme na aij . V nové matici vyhledáme ten prvek, který má největší absolutní hodnotu a znovu podle uvedených vztahů přepočítáme prvky matice B. Věta 3.10 Jacobiova metoda konverguje k diagonální matici s vlastními čísly na hlavní diagonále. 70

Důkaz. V k-tém kroku Jacobiovy metody nejprve označme n X 2 (N (A)) = |aij |2 . i,j=1, i6=j

Když apq = maxi,j a2i,j , je (N (A))2 ≤ (n2 − n)a2pq . Pro největší nediagonální prvek pak dostaneme (N (A))2 a2pq ≥ . n(n − 1) Při této volbě je bpq = bqp = 0 a platí (N (B))2 = (N (A))2 − 2a2pq ≤ q 2 (N (A))2 , kde 1 2 q = (1 − n(n−1) )2 . Když v k-tém kroku přeznačíme B na {A(k) }, pak pro posloupnost podobných matic {A(k) }, k ∈ N 1 , platí [N (A(k) )]2 ≤ q k [N (A(0) )]2 . Odtud a protože q < 1 dostaneme, že limk→∞ N (A(k) ) = 0. Algoritmus: (Jacobiova diagonalizace) Vstup: A = (ai,j )nij=1 , ε P P 2 s = ni=2 i−1 j=1 aij dokud s > ε určíme indexy p, q, p 6= q, pro které ve stávající matici apq = max |ai,j | i,j6=p,q

podle vzahů na straně 70 spočteme prvky bij matice B. A :=P B P 2 s := ni=2 i−1 j=1 aij Výstup: λ1 ≈ a1,1 , . . . , λn ≈ an,n 

 2 −1 0 Příklad Jacobiovou metodou určete vlastní čísla matice A =  −1 2 −1  . 0 −1 2 1

Zde N je množina přirozených čísel

71

Řešení: Postupně Jacobiovou metodou obdržíme  1 0 0, 7071067810  0 3 0, 7071067810 B=  0, 7071067810 0, 7071067810 2 

2, 366025404

0, 6279630301

0

  ,  

  , 0, 6279630301 3 0, 3250575835 B=   0 0, 3250575835 0, 6339745962   3, 386446078 0 0, 2768365599   , 0 1, 979579326 0, 1703641738 B=   0, 2768365599 0, 1703641738 0, 6339745962   3, 414013491 0, 01688138882 0    B=  0, 01688138882 1, 979579326 0, 1695257220  , 0 0, 1695257220 0, 6064071831   3, 414013491 0, 01675788864 −0, 002038248435   , 0, 01675788864 2, 000198603 0 B=   −0, 002038248435 0 0, 5857879069   3, 414212095 0 −0, 002038105311   , 0 1, 999999999 0, 00002415417276 B=   −0, 002038105311 0, 00002415417276 0, 5857879069  0, 585786439 −0, 00002415416649 0  1, 999999999 −0, 00000001740441820 B=  −0, 00002415416649 0 −0, 00000001740441820 3, 414213559 √ (Přesné hodnoty vlastních čísel jsou λ1,2 = 2 ± 2, λ3 = 2.)

72

  . 

Kapitola 4 Řešení nelineárních rovnic Problematiku tohoto odstavce rozdělme do dvou částí. Naším cílem bude najít • řešení nelineární rovnice f (x) = 0 nebo • řešení rovnice Pn (x) = 0, kde Pn (x) je polynom n-tého stupně. Na rozdíl od metod, které jsme používali v odstavci o soustavách lineárních rovnic, máme zde k dispozici pouze iterační metody. Metoda přímé iterace a Newtonova metoda, které si zde uvedeme pro případ jedné nelineární rovnice, se dají ve vektorové podobě použít k nalezení řešení soustavy nelineárních rovnic. Protože pro n > 1 jsou rovnice Pn (x) = 0 rovnicemi nelineárními, můžeme uvedené metody využít také ke stanovení kořenů polynomů. Rovnice Pn (x) = 0 zároveň mají i určité specifické vlastnosti, které daly vzniknout dalším speciálním metodám. Těm se budeme věnovat v druhé části této kapitoly.

Řešení nelineární rovnice f (x) = 0 Na tomto místě se budeme snažit najít takové α ∈ ha, bi, pro které f (α) = 0. Abychom si zajistili existenci alespoň jednoho reálného řešení zadané nelineární rovnice, budeme v celém tomto odstavci předpokládat, že funkce f (x) je spojitá na intervalu ha, bi a že f (a)f (b) < 0. (Viz Bolzanova věta na straně 73

59 ve skriptech [4].) Protože metody řešení nelineárních rovnic jsou iterační, je třeba nějakým způsobem získat počáteční odhad polohy kořene. Jednou z možností je použít grafickou metodu: Když položíme f (x) = f1 (x) − f2 (x) a nakreslíme grafy funkcí y = f1 (x), y = f2 (x), bude x-ová souřadnice průsečíku obou křivek odhadem řešení α rovnice f (x) = 0.

Příklad Odhadněte polohu kořene rovnice x + ln x = 0. Obrázek: Nakreslíme grafy funkcí y = x, y = − ln x. Jejich průsečík je hledaný kořen.

y=x

x0

y = − ln x

Když jsme určili, v jakém intervalu budeme řešení hledat, popř. ze kterého bodu x0 odstartujeme, budeme pokračovat některou vhodnou numerickou metodou. Numerické metody, o kterých dále budeme hovořit, rozdělme do dvou skupin, a to na startovací a zpřesňující. Mezi startovací zařadíme ty metody, které konvergují, avšak rychost konvergence bývá pomalá. Patří sem například metoda bisekce, regula falsi nebo prostá iterace. Většinou je používáme k počátečnímu odhadu řešení. Zpřesňující metody, jako Newtonova nebo Mullerova metoda, nemusejí, když zvolíme počáteční hodnotu nevhodně, konvergovat. Pokud ale konvergují, tak rychleji než metody startovací. Proto slouží hlavně ke zpřesnění řešení.

74

a) Metoda bisekce Postupným půlením intervalu ha, bi sestrojíme posloupnost interval˚ u, které obsahují hledaný kořen. Střed s intervalu ha, bi bereme za první aproximaci kořene. Pokud f (a)f (s) < 0, přeznačíme s na b, jinak přeznačíme s na a. Nový interval ha, bi opět rozp˚ ulíme atd. Chyba metody bisekce po n krocích je pak |xn − α| ≤

b−a . 2n

Obrázek: Metoda bisekce

a

x2

x1

x0

b

Algoritmus: (Metoda bisekce) Vstup: ε > 0, a, b, f (x) pro k=1,2,. . . , dokud 2b−a k−1 > ε b+a s := 2 xk := s když f (a) f (s) < 0, položíme a := s jinak b := s když f (a) f (b) = 0, skončíme Výstup: xk 75

Příklad Řešte nelineární rovnici x arctg x = 1 s přesností 10−2 . Řešení: Nejprve graficky odhadneme polohu kořene. Obrázek: Odhad polohy kořene rovnice x arctg x = 1

y = arctg(x)

1

0.8

y= 0.6

1 x

0.4

0.2

0

1

1.2

1.4

x

1.6

1.8

2

Vidíme, že kořen α leží uvnitř íntervalu h1, 2i. Metodou bisekce budeme počítat tak dlouho, než chyba bude menší než 10−2 . a=1 a=1 a=1 a = 1, 125 a = 1, 125 a = 1, 15625 a = 1, 15625 Chyba |α − s7 | ≤

1 27

b=2 b = 1, 5 b = 1, 25 b = 1, 25 b = 1, 1875 b = 1, 1875 b = 1, 171875

s1 s2 s3 s4 s5 s6 s7

= 1, 5, = 1, 25, = 1, 125, = 1, 1875, . = 1, 15625, = 1, 171875, = 1, 1640625

= 0, 007812500000.

Metoda prosté iterace Spočívá v tom, že levou stranu rovnice f (x) = 0 vyjádříme ve tvaru f (x) = x−ϕ(x). Hledané řešení pak musí splňovat x = ϕ(x). Z tohoto vztahu určíme 76

iterační formuli xk = ϕ(xk−1 ).

(4.1)

Zadáme počáteční hodnotu x0 a přesnost δ, s jakou chceme počítat. Iterujeme do té doby, než |xk − xk−1 | < δ. Zda metoda prosté iterace konverguje, závisí na tom, jak zvolíme iterační funkci y = ϕ(x). Když ji zvolíme nevhodně, nemusí metoda v˚ ubec konvergovat. Uveďme si nyní větu, která nám umožní určit podmínky, za jakých metoda prosté iterace konverguje. Věta 4.1 Nechť ϕ ∈ C 1 ha, bi, sup |ϕ0 (x)| < 1, pak rovnice x = ϕ(x) má x∈ha,bi

jediné řešení α ∈ ha, bi, ke kterému konverguje posloupnost iterací xk+1 = ϕ(xk ), k = 0, 1, . . . , pro libovolné x0 ∈ ha, bi. Pro všechna k ∈ N platí apriorní odhad chyby |xk − α| ≤

qk |x1 − x0 | 1−q

a aposteriorní odhad chyby |xk − α| ≤

q |xk − xk−1 |. 1−q

Důkaz. Když ϕ ∈ C 1 ha, bi, pak podle věty o střední hodnotě (viz str. 112 ve skriptech [4]) pro libovolné dva body x, y ∈ ha, bi, x < y, platí ϕ(x) − ϕ(y) = ϕ0 (ξ)(x − y), ξ ∈ (a, b). Odtud a z předpokladů věty dostaneme |ϕ(x) − ϕ(y)| ≤ sup |ϕ0 (ξ)| |x − y| = q|x − y|. ξ∈(a,b)

Když q < 1, je ϕ kontrakce a podle Banachovy věty 1.3 existuje jediný pevný bod tohoto zobrazení. Tvar apriorního a aposteriorního odhadu plyne z věty 1.4 77

Jako důsledek předchozí věty obdržíme tvrzení Věta 4.2 Nechť x je pevný bod zobrazení ϕ ∈ C 1 ha, bi, |ϕ0 (x)| < 1, pak metoda prostých iterací xk+1 = ϕ(xk ), k = 0, 1, . . . , je lokálně konvergentní (tj. konverguje pro každé x0 z okolí bodu x.) Podmínky, které klademe na funkci ϕ ve větě 4.2 jsou však příliš silné. Stačí, aby funkce ϕ byla spojitá na ha, bi. Pak platí Věta 4.3

Pokud ϕ ∈ C 0 ha, bi a existuje q ∈ h0, 1) tak, že |ϕ(x) − ϕ(y)| ≤ q|x − y|, x, y ∈ ha, bi,

(4.2)

pak existuje právě jeden reálný kořen rovnice x = ϕ(x) a posloupnost {xk } určená vztahem xk = ϕ(xk−1 ) konverguje k tomuto x pro libovolné x0 z intervalu ha, bi. Algoritmus: (Metoda prosté iterace) Vstup: δ, x0 , ϕ(x) pro k = 1, 2, . . . dokud |xk − xk−1 | ≥ δ xk = ϕ(xk−1 ) Výstup: xk Příklad Metodou prosté iterace řešte rovnici x + ln x = 0. Řešení: Kořen rovnice budeme hledat v intervalu h0, 1; 1i. Jednou z možností, jak vybrat iterační funkci, je položit ϕ(x) = − ln(x). Protože 1 q = sup |ϕ0 (x)| = sup | | ≥ 1, x∈h0,1;1i x∈h0,1;1i x nebudou iterace podle věty 4.1 konvergovat. Pokusme se proto vybrat iterační formuli jinak. Přepišme zadanou rovnici na tvar x = e−x . Když ϕ(xk+1 ) = e−xk , k = 0, 1, . . . , 78

bude sup |ϕ0 (x)| = sup e−x = 0.3678794412 = q < 1, x∈h0,1;1i

x∈h0,1;1i

a tato iterační formule bude konvergovat. Zvolme x0 = 1. Spočteme x1 x3 x5 x7

= 0, 5664147332, = 0, 5669089119, = 0, 5670678984, = 0, 5671860501,

x2 x4 x6 x8

= 0, 5675566373, = 0, 5672762322, = 0, 5670678984, = 0, 5671190401.

Chyba aproximace: |x − x8 | ≤

q |x8 − x7 | = 0, 00003899825913. 1−q

V případě iteračních metod není zanedbatelná otázka rychlosti konvergence. Rychlost konvergence iterační metody závisí řádu iterační metody. Metody vyšších řádů konvergují rychleji. Definice 4.1 Nechť α je přesné řešení nelineární rovnice, xk je jeho k-tá iterace, limk→∞ xk = α a c je nenulová konstanta. Řádem uvažované iterační metody budeme rozumět takové reálné číslo r ≥ 1, pro které platí lim

k→∞

|xk − α| = c. |xk−1 − α|r

Poznámka Iterační metoda určená iterační funkcí ϕ(x) je řádu r, právě když ϕ(α) = α, funkce ϕ je dost hladká a její derivace splňují ϕ(j) (α) = 0, když j ∈ h1, r), ϕ(r) (α) 6= 0.

Věta 4.4 Pokud v metodě prosté iterace xk = ϕ(xk−1 ) je ϕ0 (α) 6= 0, pak řád metody je roven jedné. Jestliže ϕ0 (α) = 0, ϕ00 (α) 6= 0, pak řád metody je roven dvěma.

79

Důkaz. Z Taylorova rozvoje funkce ϕ v bodě α máme xk − α = ϕ0 (α)(xk−1 − α) +

ϕ00 (ξ) (xk−1 − α)2 . 2

Když rovnici vydělíme (xk−1 − α) a provedeme limitní přechod pro k → ∞, je zřejmé, že se jedná o metodu prvního řádu. Ve druhém případě obdobným způsobem dokážeme, že se jedná o metodu řádu dva.

c) Metoda regula falsi Je jednou z nejstarších metod řešení nelineárních rovnic. Její princip spočívá v tom, že ke grafu funkce y = f (x) vedeme v bodech (a, f (a)), (b, f (b)) sečnu. Rovnice sečny má tvar y = f (a) +

f (b) − f (a) (x − a). b−a

Spočteme její průsečík s osou x,

f (a) +

f (b) − f (a) b−a (x − a) = 0 ⇒ x = a − f (a) . b−a f (b) − f (a)

Bod, ve kterém sečna protíná osu x, označme p a považujme jej za první aproximaci kořene α. p = a − f (a)

b−a . f (b) − f (a)

(4.3)

Pokud f (a)f (p) < 0, pak p přeznačíme na b, v opačném případě na a. Zkonstruujme novou sečnu. Proces zastavíme, když funkční hodnota v průsečíku bude menší než předem zvolené δ > 0. 80

Obrázek: Metoda regula falsi

p1 p2 a

b

Algoritmus: (Metoda regula falsi) Vstup: δ > 0, a, b, f (x) x1 := a pro k=1,2,. . . dokud f (xk ) > δ, f (a) p := a − f (b)−f (b − a) (a) xk := p když f (a) f (p) < 0, položíme b := p, jinak a := p když f (a) f (b) = 0, skončíme Výstup: xk Poznámka Metoda regula falsi konverguje pro libovolnou spojitou funkci. Ovšem v blízkosti kořene je tato konvergence velmi pomalá. Dá se ukázat, že metoda je řádu 1. Poznámka Pokud je f ∈ C 1 ha, bi, jsou splněny předpoklady Lagrangeovy věty (viz skripta [4] str. 77) a existuje tedy ξ ∈ (xk , α) ⊂ ha, bi, že f 0 (ξ) =

f (xk ) − f (α) . xk − α 81

Odtud, nezávisle na konkrétním tvaru metody, dostaneme odhad chyby |xk − α| ≤

f (xk ) , kde m1 = min |f 0 (x)|. x∈(xk ,α) m1

(4.4)

Příklad Nelineární rovnici x arctg x = 1 řešte metodou regula falsi. Řešení: V našem případě f (x) = x arctg x − 1. Kořen α budeme hledat v intervalu h1, 2i. a=1 a = 1, 150186819 a = 1, 161536517 a = 1, 162287210 a = 1, 162336388 a = 1, 162339607 a = 1, 162339818

b=2 b=2 b=2 b=2 b=2 b=2 b=2

p1 p2 p3 p4 p5 p6 p7

= 1, 150186819 = 1, 161536517 = 1, 162287210 = 1, 162336388 = 1, 162339607 = 1, 162339818 = 1, 162339832

−9

(p7 ) | = | 10 | = 0, 466 · 10−8 . Chyba |α − p7 | ≤ | fm π −1 1 4

d) Newtonova metoda Tato metoda byla objevena roku 1669. Newton ji použil při řešení kubických rovnic a při aproximaci řešení Keplerových rovnic pro pohyb planet. Pokud f ∈ C 1 ha, bi, můžeme levou stranu rovnice f (x) = 0 nahradit Taylorovým polynomem 1. stupně. Obdržíme f (x0 ) + f 0 (x0 )(x − x0 ) = 0. Když odtud spočteme x a připojíme indexy, obdržíme Newtonovu iterační formuli f (xk−1 ) xk = xk−1 − 0 . (4.5) f (xk−1 ) Iterační proces zastavíme podmínkou |xk − xk−1 | ≤ δ pro zvolené δ > 0. 82

Definice 4.2 Nechť f ∈ C 1 ha, bi, f 0 (x) 6= 0, x ∈ ha, bi, pak Newtonova metoda pro řešení rovnice f (x) = 0 je dána iteračním schématem f (xk ) , k = 0, 1, . . . , f 0 (xk )

xk+1 = xk − pro počáteční hodnotu x0 ∈ ha, bi.

Algoritmus: (Newtonova metoda) Vstup: a, b, f (x), x0 , δ > 0 pro k = 1, 2, . . . , dokud |xk − xk−1 | > δ f 0 (xk−1 ) k−1 ) xk = xk−1 − ff0(x (xk−1 ) Výstup: xk Poznámka Newtonova metoda bývá někdy nazývána metodou tečen. To proto, že bod (xk , 0) je průsečíkem tečny sestrojené v bodě (xk−1 , f (xk−1 )) ke křivce y = f (x). Obrázek: Newtonova metoda

a

x2

x1 x0 = b

Poznámka Z Taylorova rozvoje se dá odvodit (viz [3]), že pro odhad chyby 83

Newtonovy metody po n krocích platí 1 k |xn − α| ≤ |c(x0 − α)|2 , c

(4.6)

00

(x1 ) kde c ≥ 12 ff 0 (x pro libovolné x1 , x2 ∈ (a, b). Je vidět, že metoda bude určitě 2) konvergovat, když |c(x0 − α)| < 1, tj. když počáteční aproximace x0 je dost blízká kořenu α. Pokud počáteční aproximace x0 není dost blízká kořenu α rovnice f (x) = 0, nemusí metoda konvergovat.

Uveďme si příklady některých podmínek, za kterých Newtonova metoda konverguje. Věta 4.5 Když f ∈ C 2 (a, b) a kořen α rovnice f (x) = 0 je jednoduchý (tj, f 0 (α) 6= 0), pak Newtonova metoda je lokálně konvergentní. Věta 4.6 Jestliže f 0 (x) 6= 0, x ∈ ha, bi a f 00 (x) v intervalu ha, bi nemění znaménko, f (a)f (b) < 0, ff0(a) < b − a a ff0(b) < b − a, pak Newtonova metoda (a) (b) konverguje pro libovolné x0 . Příklad Newtonovou metodou řešte rovnici x − e−x = 0. (Volte x0 = 1, δ = 10−4 .) Řešení: V našem případě f (x) = x − e−x , f 0 (x) = 1 + e−x . x0 = 1 x1 = 0, 5378828427 x2 = 0, 5669869914 x3 = 0, 5671432860 x4 = 0, 5671432904. Když tyto výsledky porovnáme s výsledky příkladu na str. 78, vidíme, že rychlost Newtonovy metody je podstatně vyšší než u metody prosté iterace. Poznámka Newtonova metoda je metodou řádu 2. Kvadratická konvergence Newtonovy metody znamená, že počet platných číslic v numerické aproximaci se v každém iteračním kroku zdvojnásobí.

84

e) Metoda sečen Iterační formuli odvodíme tak, že v rovnici f (x) = 0 funkci f (x) nahradíme lineárním polynomem v bodech xk−1 , xk . Z nové rovnice f (xk ) − f (xk−1 ) (x − xk ) + f (xk ) = 0 xk − xk−1 spočteme x. Když připojíme indexy, dostaneme iterační formuli xk+1 = xk − f (xk )

xk − xk−1 . f (xk ) − f (xk−1 )

(4.7)

Obdrželi jsme metodu, kterou je také možné získat z Newtonovy metody (4.5) tím, že derivaci f 0 (x) nahradíme poměrnou diferencí. Metoda sečen je příkladem dvoubodové metody, neboť k jejímu odstartování potřebujeme zadat dvě počáteční hodnoty x0 , x1 . Jestliže je zvolíme špatně, nemusí metoda konvergovat. Pokud však metoda sečen konverguje, dá se ukázat, že metoda je přibližně řádu 1,6. Algoritmus: (Metoda sečen) Vstup: a, b, f (x), x0 , x1 , δ > 0 pro k = 1, 2, . . . , dokud |xk+1 − xk | ≥ δ −xk−1 xk+1 = xk − f (xk ) f (xxkk)−f (xk−1 ) Výstup: xk+1 Příklad Metodou sečen řešte rovnici x − e−x = 0. Řešení: Když zvolíme x0 = 2, x1 = 1 a δ = 10−4 , pak x2 = 0, 4871416535, x3 = 0, 5730763105, x4 = 0, 5672301456, x5 = 0, 5671431972.

85

Řešení rovnic Pn(x) = 0 Přímo řešit umíme jen několik málo typ˚ u nelineárních algebraickch rovnic Pn (x) = 0. Patří sem rovnice kvadratické a2 x2 + a1 x + a0 = 0, kubické a3 x3 + a2 x2 + a1 x + a0 = 0, binomické xn − a = 0, trinomické ax2n + bxn + c = 0, reciproké a0 xn + a1 xn−1 + · · · = 0, kde ak = an−k nebo ak = −an−k . U zbývajících typ˚ u určujeme kořeny numericky. K tomu použijeme buď metody, kterými jsme se zabývali v předchozí části odstavce, nebo metody určené speciálně k hledání kořenů polynomu Pn (x). V první fázi bychom se měli pokusit získat předběžné informace o poloze, popř. počtu kořenů polynomu. Věta 4.7 Je dán polynom Pn (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an , kde a0 , . . . an ∈ R. Pokud a0 6= 0, an 6= 0, pak pro jeho kořeny α (reálné i komplexní) platí (1 +

max{|an−1 |, . . . , |a0 |} −1 max{|an |, . . . , |a1 |} ) < |α| < 1 + . |an | |a0 |

Věta 4.8 (Budan-Fourier) Nechť a, b ∈ R, a < b, koeficient an polynomu Pn (x) je kladný a Pn (a) · Pn (b) 6= 0. Když V (a) je počet znaménkových změn (n) v posloupnosti Pn (a), Pn0 (a), . . . , Pn (a) a V (b) je počet znaménkových změn (n) v posloupnosti Pn (b), Pn0 (b), . . . , Pn (b), pak počet reálných kořenů polynomu Pn (x) v intervalu (a, b) (včetně jejich násobnosti) je roven rozdílu V (a)−V (b) popř. o sudý počet méně. 86

Poznámka Následující věta je důsledkem věty 4.8. Počet kladných kořenů resp. záporných kořenů je určen počtem znaménkových změn v posloupnosti koeficientů polynomu Pn (x) resp. Pn (−x). Věta 4.9 (Descartes) Počet kladných reálných kořenů polynomu Pn (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an (včetně jejich násobnosti) je roven počtu znaménkových změn v posloupnosti koeficientů a0 , a1 , . . . , an nebo o sudý počet méně. (Nulové koeficienty neuvažujeme). Počet záporných kořenů je roven počtu znaménkových změn v posloupnosti (−1)n a0 , (−1)n−1 a1 , . . . , an nebo o sudý počet méně.

Příklad Odhadněte počet a polohu kořenů polynomu P (x) = x3 − 2x + 1. Řešení: Podle věty 4.7 je 1 max{1, 2} −1 max{2, 1} = (1 + ) < |α| < 1 + =3 3 1 1 1 1 x ∈ (−3, − ) ∪ ( , 3). 3 3 Podle věty 4.9 má uvažovaný polynom 2 nebo žádný reálný kladný kořen a 1 reálný záporný kořen. Dále P 0 (x) = 3x2 − 2, P 00 (x) = 6x, P 000 (x) = 6, x P (x) P 0 (x) P 00 (x) P 000 (x) V (x)

0 + − 0 + 2

1 0 + + + 0

2 + + + + 0

3 + + + + 0

Protože P (1) = 0, je α = 1 kořenem zadaného polynomu. Odtud a z věty 4.8 je vidět, že polynom má v intervalu (0, 1i dva reálné kořeny. 87

Pokud známe některý kořen polynomu Pn (x) = a0 xn + a1 xn−1 + · · · + an , je možné snížit stupeň polynomu prostřednictvím Hornerova schématu. Jestliže příslušný kořen označíme α1 , dostaneme Pn (x) = (x − α1 )Pn−1 (x), kde Pn−1 (x) = b0 xn−1 + b1 xn−2 + · · · + bn−1 . Pokud α2 je kořenem polynomu Pn−1 (x), můžeme znovu snížit stupeň polynomu pomocí Hornerova schématu atd. Musíme si však uvědomit, že během výpočtu dochází k zaokrouhlování. Pokud hodnoty kořenů jsou přibližné, mohou být koeficienty nových polynomů Pn−1 , Pn−2 , . . . zatížené chybami. V takovém případě nemusíme dostat dobrou aproximaci dalšího kořene polynomu Pn . Dělení lineárními činiteli x − αk , k = 1, . . . , n − 1, přinese dobré výsledky, když |α1 | < |α2 | < . . . |αn−1 |, tj. pokud začínáme dělit kořenovým činitelem příslušným nejmenšímu kořenu. Algoritmus: (Dělení lineárním polynomem - Hornerovo schéma) Vstup: n, α, a0 , . . . , an b0 = a0 pro k = 1, 2, . . . , n bk = αbk−1 + ak Výstup: b1 , b2 , . . . , bn Snížení stupně polynomu je možné dosáhnout také tak, že vydělíme kvadratickým činitelem d(x) = x2 − px − q. Potom Pn (x) = d(x)Qn−2 (x) + ϕ(x), kde Qn−2 (x) = b0 xn−2 + b1 xn−1 + · · · + bn−2 , ϕ(x) = bn−1 x + bn . Algoritmus: (Dělení kvadratickým polynomem) Vstup: n, p, q, a0 , . . . , an b−1 = 0, b0 = a0 pro k = 1, 2, . . . , n − 1 bk = ak + pbk−1 + qbk−2 bn = ak + qbn−2 Výstup: b1 , b2 , . . . , bn

88

a) Bernoulliova metoda Metoda je určena k vyhledání reálného dominantního kořene polynomu. Věta 4.10 Když polynom Pn (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an má reálné různé kořeny α1 > α2 > · · · > αn a yk =

−1 (a1 yk−1 + a2 yk−2 + · · · + an yk−n ) a0

je řešením diferenční rovnice a0 yk + a1 yk−1 + a2 yk−2 + · · · + an yk−n = 0 doplněné podmínkami y0 = y1 = · · · = yn−2 = 0, yn−1 = 1, pak

yk+1 . k→∞ yk

α1 = lim

(4.8)

Důkaz. Algebraická rovnice a0 xn + a1 xn−1 + · · · + an−1 x + an = 0 je charakteristickou rovnicí příslušnou k diferenční rovnici a0 yk + a1 yk−1 + a2 yk−2 + · · · + an yk−n = 0. Když α1 , . . . , αn jsou reálné různé charakteristické rovnice, má řešení Pn kořeny k diferenční rovnice tvar yk = i=0 ci αi . (Konstanty ci vypočteme z počátečních podmínek.) Pro c1 α1 6= 0 je yk =

c1 α1k (1

n X ci αi k + ( ) ). c α1 i=2 1

Za předpokladu, že α1 je největší vlastní číslo a i > 1, je limk→∞ ( αα1i )k = 0. Pak yk+1 c1 α1k+1 lim = lim = α1 . k k→∞ yk k→∞ c1 α1

89

Algoritmus: (Bernoulliova metoda) Vstup: n, a0 , . . . , an , m y0 = y1 = · · · = yn−2 = 0, yn−1 = 1 pro k = n, n +P 1, . . . , m yk = − a10 nj=1 ai yk−i yk α1 = yk−1 Výstup: α1 Příklad Najděte dominantní kořen rovnice x3 − x2 −

x 4

+

1 4

= 0.

Řešení: Položíme y0 = y1 = 0, y2 = 1 a pro k = 4, 5, . . . spočítáme iterace yk = yk−1 + 14 yk−2 − 14 yk−3 . y3 y4 y5 y6 y7 y8 y9

=1 = 1, 250000000, = 1, 250000000, = 1, 312500000, = 1, 312500000, = 1, 328125000, = 1, 328125000,

α11 α12 α13 α14 α15 α16 α17

=1 = 1, 25, = 1, = 1, 05, = 1, = 1, 011904762, = 1.

Dominantní kořen α1 = lim α1n = 1. n→∞

Poznámka Existuje i modifikace Bernoulliovy metody, kterou lze použít k nalezení dominantních komplexně sdružených kořenů.

b) Graefova metoda Naším cílem je najít všechny kořeny polynomu za předpokladu, že tyto kořeny jsou reálné a navzájem různé. Při odvozování metody vyjdeme ze vztahu, který platí mezi kořeny a koeficienty kvadratické rovnice. Když rovnice má tvar a0 x2 + a1 x + a2 = 0, pak pro její kořeny α1 , α2 platí α1 + α2 = −

a2 a1 , α1 α2 = . a0 a0 90

(4.9)

První z těchto rovnic lze přepsat na tvar α1 (1 +

α2 −a1 )= . α1 a0

Když předpokládáme, že | αα21 | << 1, můžeme kořeny zadané kvadratické rovnice přibližně spočítat pomocí vztahů α1 ≈

a2 −a0 −a2 −a1 , α2 ≈ ( )= . a0 a0 a1 a1

Sestavme nyní novou kvadratickou rovnici b0 x2 + b1 x + b2 = 0, jejíž kořeny jsou kvadrátem kořenů α1 , α2 . Platí β1 + β2 = α12 + α22 = (α1 + α2 )2 − 2α1 α2 =

a21 a2 a21 − 2a2 a0 b1 = − 2 =− , 2 2 a0 a0 a0 b0

a22 b2 = . 2 a0 b0 Koeficienty nové kvadratické rovnice proto splňují β1 β2 = α12 α22 =

b0 = a20 , b1 = −(a21 − 2a0 a2 ), b2 = a22 a pro kořeny platí β1 ≈

−b1 b2 , β2 ≈ . b0 b1

Kořeny původní rovnice s

s |α1 | ≈

|b1 | , |α2 | ≈ |b0 |

|b2 | . |b1 |

Této úvahy lze využít při hledání reálných kořenů polynomu Pn (z) = a0 xn + a1 xn−1 + · · · + an , pokud jsou tyto kořeny uspořádány sestupně podle absolutních hodnot jejich velikostí. Definujme posloupnost polynomů P 0 (x), P 1 (x), . . . , P n (x) tak, že P 0 (x) = Pn (x), P k (x) = ak0 xn + ak1 xn−1 + · · · + akn−1 x + akn , 91

kde a0i = ai , i = 0, 1, . . . , n, a pro k = 1, 2, . . . , n položme 2 ak0 = (ak−1 0 ) , 2 k−1 k−1 −ak1 = (ak−1 1 ) − 2a0 a2 , k−1 2 k−1 k−1 + 2ak−1 ak2 = (ak−1 0 a4 , 2 ) − 2a1 a3 2 k−1 k−1 k−1 k−1 −ak3 = (ak−1 + 2ak−1 − 2ak−1 3 ) − 2a2 a4 1 a5 0 a6 , ...... 2 k−1 k−1 (−1)n−1 akn−1 = (ak−1 n−1 ) − 2an−2 an , 2 (−1)n akn = (ak−1 n ) .

Pro kořeny původní rovnice dostáváme s |αj | ≈

2k

akj | k |, j = 1, 2, . . . , n. aj−1

(4.10)

Algoritmus: (Graefova metoda) Vstup: n, a0 , . . . , an , m pro j = 1, 2, . . . , n a0j = aj pro k = 1, 2, . . . , m Pj−1 k−1 akj = (−1)j [(ak−1 )2 + 2 s=1 (−1)s (ak−1 s−1 )as+1 j r ak k |αj | ≈ 2 | ak j | j−1

Výstup: |α1 |, . . . , |αn | Příklad Graefovou metodou najděte kořeny polynomu x2 − 1, 5x + 0, 5. Řešení: Výpočet pro přehlednost zaznamenáme do tabulky 92

k 0 1 2

ak0 ak1 ak2 1 −1, 5 0, 5 1 −1, 25 0, 25 1 −1, 0625 0, 0625

Z hodnot v posledním řádku (viz vztahy (4.10)) dostaneme odhad kořenů s p . 4 a21 |α1 | = | 2 | = 4 1, 0625 = 1, 015915975, a0 s r 0, 0625 . 4 a22 |α2 | = | 2 | = 4 = 0, 4924790605. a1 1, 0625 Poznámka Předpoklady, za kterých jsme konstruovali metodu, připouštějí její aplikaci pouze na případ, že polynom má reálné r˚ uzné kořeny. Uvedená metoda se ale dá modifikovat tak, že ji lze použít i ke stanovení páru kořen˚ u komplexních.

c) Laguerrova metoda Je to hojně používaná iterační metoda, která slouží k výpočtu kořenů algebraické rovnice Pn (x) = 0. Metoda je velice univerzální, protože konverguje jak v případě reálných, tak v případě komplexních kořenů. Konvergence je zvlášť rychlá, když kořeny jsou reálné jednoduché. Předpokládejme, že α1 ≤ α2 ≤ · · · ≤ αn jsou kořeny algebraické rovnice Pn (x) = 0. Počáteční aproximaci x0 ∈ R zvolme libovolně. Další aproximace kořene určíme ze vztah˚ u nP (xk ) p , P 0 (xk ) ± H(xk ) H(xk ) = (n − 1)[(n − 1)(P 0 (xk ))2 − nP (xk )P 00 (xk )]. xk+1 = xk −

(4.11)

Znaménko „+ÿ ve jmenovateli zlomku volíme v případě, že P (xk ) > 0. Jestliže x0 ∈ (αj−1 , αj ), pak hodnoty xk+1 konvergují k jednomu z kořenů αj−1 nebo αj . 93

Algoritmus: (Laguerova metoda) Vstup: n, a0 , . . . , an , x0 , m pro k = 1, 2, . . . , m H(xk ) = (n − 1)[(n − 1)(P 0 (xk ))2 − nP (xk )P 00 (xk )] √k ) xk+1 = xk − 0 nP (x P (xk )±

H(xk )

Výstup: xm Příklad Prostřednictím Laguerrovy metody najděte některý z kořenů polynomu x3 − x2 − x4 + 14 = 0. Řešení: Když položíme x0 = 0, 1, pak x1 = −0, 2333518266, x2 = −0, 4939586061, x3 = −0, 5011847770, x4 = −0, 4997504466, x5 = −0, 5000498708, x6 = −0, 4999900040, x7 = −0, 5000019991.

Obrázek: Graf zadaného polynomu 2 1.8 1.6

p(x) = x3 − x2 −

1.4 1.2 y1 0.8 0.6 0.4 0.2 –1–0.8

–0.4 –0.2

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x

–0.4 –0.6 –0.8 –1 –1.2 –1.4

94

x 4

+

1 4

Kapitola 5 Aproximace funkcí Častokrát máme co do činění s funkcemi, které jsme získali jako výsledky laboratorních měření a u kterých proto máme k dispozici jen jejich hodnoty v diskrétních bodech, popřípadě s funkcemi, které jsou příliš složité a obtížně se s nimi pracuje. V takových případech se můžeme pokusit nahradit funkci f jinou, jednodušší funkcí ϕ. Způsobů, jak zkonstruovat aproximující funkci, je celá řada. Uveďme si některé z nich. Je dáno n + 1 uzlových bodů x0 , x1 , . . . , xn ∈ ha, bi a n + 1 funkčních hodnot f (x0 ), f (x1 ), . . . , f (xn ). Podle toho, jaké požadavky klademe na funkci ϕ, můžeme rozlišit následující typy aproximací • Interpolace: Funkci ϕ vytvoříme tak, aby pro všechna i = 0, 1, . . . , n splňovala interpolační podmínky 1 ϕ(xi ) = f (xi ).

(5.1)

• Aproximace metodou nejmenších čtverců: V uzlových bodech minimalizujeme součet kvadrátů odchylek funkce f od funkce ϕ. Tzn. minimalizujeme n X

(f (xi ) − ϕ(xi ))2 .

(5.2)

i=0

• Čebyševova (stejnoměrná) aproximace: Žádáme, aby max |(f (xi ) − ϕ(xi ))|

xi ∈ha,bi 1

(5.3)

Splnění interpolačních podmínek je možné vyžadovat i od derivací aproximující funkce. Viz Hermitova interpolace na str. 119.

95

bylo co nejmenší. Poznámka V posledních dvou případech se jedná o minimalizaci kvadrátu normy v prostoru l2 ha, bi a normy v prostoru C 0 ha, bi.

Interpolace Už několik století se provádí interpolace pomocí polynomů. Poprvé bylo interplačních polynomů použito v 17. století při tabelování logaritmů. Proč se právě polynomy těší takové oblibě, vyplývá z jejich vlastností. Polynom je funkce, kterou lze získat z konečného počtu vstupních údajů pomocí konečného počtu operací sčítání a násobení. Lehce se derivuje a integruje. Navíc je možné polynomem aproximovat spojité funkce, a to s libovolnou přesností, jak ukazuje následující věta: Věta 5.1 (Weierstrass) Když f (x) ∈ C 0 ha, bi, pak ke každému ε > 0 existuje n a polynom n-tého stupně p(x), tak, že max |f (x) − p(x)| < ε.

x∈ha,bi

Nechť tedy funkce, kterou chceme aproximovat, je zadaná funkčními hodnotami v n + 1 uzlových bodech: x0 x1 . . . xn−1 xn f0 f1 . . . fn−1 fn Funkční hodnota v bodě xi je v tabulce označena jako fi . Dále označme Pn množinu všech polynomů n-tého stupně. Bázi tohoto prostoru tvoří polynomy {1, x, x2 , . . . , xn }. Naším cílem bude zkonstruovat interpolační polynom p(x) ∈ Pn .

96

a) Lagrangeův interpolační polynom Nejprve si uveďme postup, který při konstrukci interpolačního polynomu v roce 1794 použil francouzský matematik J.L.Lagrange. Věta 5.2 (Lagrangeova konstrukce interpolačního polynomu) Když je dáno n+1 uzlových bodů x0 , x1 , . . . , xn ∈ ha, bi a n+1 funkčních hodnot f0 , f1 , . . . , fn−1 , fn ∈ R, pak existuje právě jeden polynom p ∈ Pn tak, že p(xi ) = fi . V Lagrangeově reprezentaci má tento polynom tvar L(x) =

n X

li (x)fi ,

(5.4)

i=0

kde

n Y (x − xj ) li (x) = . (xi − xj ) j=0,

(5.5)

j6=i

Důkaz. Existence: Především je vidět, že pro všechna i = 1, . . . ,P n jsou li (x) ∈ Pn . Dále pro i, j = 1, . . . , n platí li (xj ) = δij . Pak ale L(xj ) = i δij fi = fj . Tzn. polynom L splňuje interpolační podmínky (5.1). Jednoznačnost: Předpokládejme, že existují dva různé interpolační polynomy L1 , L2 ∈ Pn , L1 6= L2 . Pak polynom L = L1 − L2 splňuje v uzlových bodech podmínky L(xj ) = 0, j = 0, . . . , n. Je tedy L polynom stupně n, který má n + 1 kořenů. To může nastat pouze tehdy, když všechny jeho koeficienty jsou nulové. To znamená, že L = L1 − L2 = 0, a tedy L1 = L2 . Algoritmus: (Lagrange˚ uv interpolační polynom) Vstup: n, x0 , x1 , . . . , xn , f0 , f1 , . . . , fn pro i = 0, 1, Q ...,n (x−x ) li (x) = nj=0,i6=j (xi −xjj ) Pn+1 L(x) = i=0 li (x)fi Výstup: L(x) 97

Příklad Najděte Lagrangeův interpolační polynom funkce zadané tabulkou xi −1 0 1 3 fi 2 1 2 0 Řešení: l0 (x) = l1 (x) =

x(x − 1)(x − 3) x3 − 4x2 + 3x = , −1(−2)(−4) −8

(x + 1)(x − 1)(x − 3) x3 − 3x2 − x + 3 = , 1(−1)(−3) 3

l2 (x) =

(x + 1)x(x − 3) x3 − 2x2 − 3x = , 2.1(−2) −4

(x + 1)(x − 1)x x3 − x = . 4.3.2 24 1 L(x) = 2l0 (x) + l1 (x) + 2l2 (x) = (−5x3 + 12x2 + 5x + 12). 12 l3 (x) =

Poznámka

Při označení qn+1 (x) =

n Y

(x − xj )

j=0

lze psát Lagrangeovy polynomy li ve tvaru li (x) =

qn+1 (x) . 0 (x − xi )qn+1 (xi )

Věta 5.3 (Odhad chyby) Když f : ha, bi → R, f ∈ C n+1 (a, b) a L je interpolační polynom zkonstruovaný v bodech x0 , . . . , xn ∈ ha, bi, pak pro všechna x ∈ ha, bi je n f (n+1) (ξ) Y f (x) − L(x) = (x − xj ), kde ξ ∈ ha, bi. (n + 1)! j=0

98

(5.6)

Důkaz. Je vidět, že v uzlových bodech x0 , . . . , xn je f (xi ) − L(xi ) = 0, a tedy uvedená formule zde platí. Označme qn+1 (x) =

n Y

(x − xj ).

j=0

Pro pevné x, x 6= xi , i = 0, 1, . . . , n, uvažujme funkci g : ha, bi → R takovou, že f (x) − L(x) g(y) = f (y) − L(y) − qn+1 (y) , y ∈ ha, bi. qn+1 (x) Protože podle předpokladů věty f ∈ C n+1 (a, b), je také g ∈ C n+1 (a, b). To ale znamená, že g musí mít alespoň n + 1 kořenů, g 0 má alespoň n kořenů a g (n+1) alespoň jeden kořen ξ ∈ ha, bi. Pro tento kořen platí 0 = f (n+1) (ξ) − (n + 1)!

f (x) − L(x) . qn+1 (x)

Odtud plyne tvrzení. Poznámka Protože bod ξ ve vzorci (5.6) nelze obecně stanovit, odhadujeme chybu interpolačního vzorce vztahem ||f (x) − L(x)||∞ ≤

1 ||qn+1 (x)||∞ ||f (n+1) (x)||∞ . (n + 1)!

(5.7)

Z tohoto vztahu je vidět, že chyba interpolace závisí na vlastnostech interpolované funkce. Protože vztah (5.7) obsahuje výraz qn+1 , je chyba interpolace ovlivněna také rozložením uzlů. Poznámka Kromě chyby metody může být výpočet ovlivněn ještě zaokrouhlovacími chybami. Pokud jsou hodnoty v uzlech xi zatíženy chybami, pak je hodnota aproximace f˜i = fi + ∆fi , i = 0, . . . , n, a Lagrangeův interpolační polynom n X i=0

li (x)f˜i =

n X

li (x)fi +

i=0

n X i=0

99

li (x)∆fi .

Dále |

n X i=0

li (x)∆fi | ≤

n X

|li (x)| |∆fi | ≤ c

i=0

n X

|li (x)|,

i=0

P P kde c ≥ |∆fi |, i = 0, 1, . . . , n. Protože ni=0 li (x) = 1, je ni=0 |li (x)| ≥ 1 a můžeme učinit Pzávěr, že s rostoucím počtem uzlů také pomalu roste koeficient šíření chyby ni=0 |li (x)|.

Poznámka Předpokládejme nyní, že uzly interpolace jsou v intervalu ha, bi rozloženy ekvidistantně (tj. jsou od sebe stejně vzdálené) s krokem h. Pro uzlové body pak máme xi = x0 + ih, i = 0, 1, . . . , n. Pokud označíme relativní vzdálenost bodu x od bodu x0 symbolem t, x − x0 t= , h pak x = x0 + th a pro koeficienty Lagrangeova polynomu (viz vzorec (5.5)) dostaneme n n n Y Y (x − xj ) (x0 + th − x0 − jh) Y (t − j) li (x0 + th) = = = = li (t). (x − x ) (x + ih − x − jh) (i − j) i j 0 0 j=0, j=0, j=0, i6=j

i6=j

i6=j

To znamená, že v případě ekvidistantních uzlů polynom li závisí na relativní vzdálenosti t, nikoliv na skutečné velikosti kroku h. Lagrangeův interpolační polynom (5.4) L(x0 + th) =

n X

li (t)fi = L(t)

i=0

je pak také funkcí závisou na t. Podle (5.6) je rozdíl n n f (n+1) (ξ) n+1 Y f (n+1) (ξ) Y (t − j)h = h (t − j). f (x0 + th) − L(x0 + th) = (n + 1)! j=0 (n + 1)! j=0

(5.8) Vidíme, že čím je krok h větší, tím je větší chyba metody. Zaokrouhlovací chyby ovšem nejsou podle předchozí poznámky na kroku h závislé.

100

b) Newtonův interpolační polynom Lagrangeův interpolační vzorec má význam především při teoretickém zkoumání vlastností interpolačních polynomů. Pro praktické počítání ale není příliš výhodný. Například když se rozhodneme přidat další uzlový bod, musíme znovu přepočítat všechny polynomy li (x). V takovém případě je výhodnější použít Newtonův interpolační vzorec. Newtonův interpolační polynom předpokládejme ve tvaru N (x) = a0 + a1 (x − x0 ) + · · · + an (x − x0 ) . . . (x − xn−1 ).

(5.9)

Když požadujeme, aby pro všechny uzlové body platilo N (xi ) = fi , obdržíme soustavu lineárních rovnic N (x0 ) = a0 = f0 , N (x1 ) = a0 + a1 (x1 − x0 ) = f1 , N (x2 ) = a0 + a1 (x2 − x0 ) + a2 (x2 − x0 )(x2 − x1 ) = f2 , ...... Z ní určíme koeficienty a 0 = f0 , f1 − a0 f1 − f0 a1 = = , x1 − x0 x1 − x0 [f2 − a0 − a1 (x1 − x0 )] + [a1 (x1 − x0 ) − a1 (x2 − x0 )] f2 − a0 − a1 (x2 − x0 ) = = a2 = (x2 − x0 )(x2 − x1 ) (x2 − x0 )(x2 − x1 ) f1 −f0 (x1 − x2 ) 1 f2 − f1 1 f2 − f1 f1 − f0 x1 −x0 = ( + )= ( − ), x2 − x0 x2 − x1 x2 − x1 x2 − x0 x 2 − x1 x1 − x0 ...... Pokud označíme Di0 = fi , i = 0, . . . , n, a Dik =

k−1 Dik−1 − Di−1 , k = 1, 2, . . . n, i = k, . . . , n, xi − xi−k

(5.10)

pak ak = Dkk a Newtonův interpolační polynom můžeme psát ve tvaru N (x) =

D00

+

n X

Dkk

k=1

101

k−1 Y

(x − xi ).

i=0

(5.11)

Výpočet koeficientů při ručním počítání provádíme v tabulce: x0 x1 .. .

D00 = f (x0 ) D10 = f (x1 ) .. .

D11 .. .

0 1 xn−1 Dn−1 = f (xn−1 ) Dn−1

xn

Dn0 = f (xn )

Dn1

.. . .. . .. .

.. .

.. .

n−1 Dn−1

Dnn−1 Dnn

Věta 5.4 V Newtonově reprezentaci má interpolační polynom z věty 5.2 tvar (5.11). Poznámka Výhodou Newtonovy konstrukce je, že když přidáme další uzel, stačí dopočítat pouze jeden řádek tabulky. Poznamenejme ještě, že Lagrangeova i Newtonova konstrukce vedou po úpravě na jeden a tentýž interpolační polynom. Algoritmus (Výpočet koeficientů Newtonova polynomu) Vstup: n, x0 , . . . , xn , f0 , . . . , fn pro i = 0, 1, . . . , n Di0 = fi pro k = 1, 2, . . . , n pro i = k, k + 1, . . . , n Dik =

k−1 Dik−1 −Di−1 xi −xi−k

Výstup: a0 = D00 , a1 = D11 , . . . , an = Dnn Poznámka Newtonův polynom lze psát N (x) = (. . . (an (x − xn−1 ) + an−1 )(x − xn−2 ) + · · · + a1 )(x − x0 ) + a0 . Je to stejný tvar, jaký se používá při odvozování Hornerova schématu. Při určování Newtonova interpolačního polynomu je proto, stejně jako v Hornerově algoritmu, zapotřebí O(n2 ) operací.

102

Příklad Určete Newton˚ uv interpolační polynom, když je funkce zadána tabulkou xi −1 0 1 3 fi 2 1 2 0 Řešení:

−1 0 1 3

2 1 −1 2 1 1 5 0 −1 − 23 − 12

N (x) = 2 − (x + 1) + (x − 0)(x + 1) − 5 5 3 x + x2 + 12 x + 1. = − 12

5 (x 12

− 0)(x + 1)(x − 1) =

Poznámka Jistou analogií uvedené Newtonovy konstrukce je tzv. Nevill˚ uv algoritmus, který umožňuje stanovit funkční hodnotu interpolačního polynomu v bodě α, aniž bychom znali konkrétní tvar interpolačního polynomu. Označme pki hodnoty Lagrangeova interpolačního polynomu k-tého stupně v bodě α. (Polynom jsme získali interpolací funkce f v uzlech xi−k , xi−k+1 , . . . , xi , i = 0, 1, . . . , n, k ≤ i.) Pak p0i (x) = f (xi ), i = 0, . . . , n, (x − xi )p0i+1 (x) − (x − xi+1 )p0i (x) 1 = pi (x) = xi+1 − xi [(x − xi+1 ) + (xi+1 − xi )]p0i+1 (x) − (x − xi+1 )p0i (x) = = xi+1 − xi p0 (x) − p0i−1 (x) = p0i + (x − xi ) i , i = 0, . . . , n − 1, xi − xi−1 (x − xi )p1i+1 (x) − (x − xi+2 )p1i (x) p2i (x) = = xi+2 − xi p1 (x) − p1i−1 (x) , i = 0, . . . , n − 2, = p1i (x) + (x − xi ) i xi − xi−1 ............ (x) + (x − xi ) pki (x) = pk−1 i

pk−1 (x) − pk−1 i i−1 (x) , i = k, . . . , n. xi − xi−1 103

Pro hodnotu interpolačního polynomu v bodě α platí p(α) = pnn (α). Algoritmus (Neville - výpočet hodnoty interpolačního polynomu) Vstup: n, x0 , . . . , xn , f0 , . . . , fn , α pro i = 0, 1, . . . , n p0i = fi pro k = 1, 2, . . . , n pro i = k, k + 1, . . . , n pki = pk−1 + (α − xi ) i

k−1 pk−1 −pi−1 i xi −xi−k

Výstup: p(α) = pnn Příklad Určete hodnotu interpolačního polynomu v bodě α = 2, když funkce je zadaná tabulkou xi −1 0 1 3 fi 2 1 2 0 Řešení: Výpočet provedeme v tabulce xi α − xi −1 3 0 2 1 1 3 −1

fi p1i p2i p3i 2 1 −1 2 3 5 0 1 5/3 5/2

Zde 1−2 = −1, 0 − (−1) 2−1 = 3, p12 = 2 + 1 1 − (0) 0−2 p13 = 0 − 1 = 1, 3−1 3+1 p22 = 3 + 1 = 5, atd. 1 − (−1) p11 = 1 + 2

104

Hodnota interpolačního polynomu p(2) = 52 . Při ekvidistantních uzlech můžeme Newtonův interpolační polynom zapsat ve tvaru, v němž vystupují diference vpřed a vzad. V uzlech xi = x0 + ih, i = 0, . . . , n, definujeme diference vpřed následovně ∆0 fi = fi , ∆1 fi = fi+1 − fi , ∆2 fi = ∆(fi+1 − fi ) = fi+2 − 2fi+1 + fi , ......... ¶ k µ X k k k−1 ∆ fi = ∆(∆ fi ) = fi+k−j . j j=0

Všimněme si, že k-tá diference vpřed v bodě xi závisí na funkčních hodnotách fi , fi+1 , . . . , fi+k , které leží napravo od bodu xi . Protože Dik =

∆k fi , k!hk

můžeme psát, že Newtonův interpolační polynom µ ¶ µ ¶ n µ ¶ X t t t n ∆i f0 . ∆ f0 = ∆f0 + · · · + N (x0 + th) = f0 + i n 1 i=0 +

(5.12)

V případě, že vycházíme při konstrukci interpolačního polynomu z uzlů na pravém konci interpolačního intervalu, můžeme využít diferencí vzad. Tyto diference definujeme vztahy ∇0 f i = f i , ∇1 fi = fi − fi−1 , ∇2 fi = ∇(fi − fi−1 ) = fi − 2fi−1 + fi−2 , ......... ¶ k µ X k k k−1 ∇ fi = ∇(∇ fi ) = fi−k+j . j j+0

105

Tentokrát ale k-tá diference vzad v bodě xi závisí na hodnotách fi , fi−1 , . . . , fi−k , které leží vlevo od bodu xi . Protože mezi diferencemi vpřed a vzad platí: ∆k fi−k = ∇k fi , můžeme Newtonův interpolační polynom psát ve tvaru: µ

¶ µ ¶ µ ¶ n X −t −t n i −t N (x) = fn + ∇fn + · · · + ∇ f0 = (−1) ∇i f 0 . 1 n i i=0 −

Extrapolace O interpolaci hovoříme, pokud určujeme funkční hodnoty uvnitř interpolačního intervalu. V případě, že počítáme funkční hodnoty vně tohoto intervalu, hovoříme o extrapolaci. Extrapolace však dává uspokojivé výsledky pouze v těch případech, kdy nejsme příliš vzdáleni od koncových bod˚ u zadaného intervalu. Příklad Když znáte hodnoty funkce sin x pro x = 20◦ a 21◦ , určete (interpolací) sin 20◦ 180 a (extrapolací) sin 21◦ 180 . xi 20 21 sin xi 0, 34202 0, 35837 Řešení: Při interpolaci máme f (x) ∼ L1 (x) = sin 20◦ 180 ∼

x − x1 x − x0 f0 + f1 , x0 − x1 x1 − x 0

20, 3 − 20 20, 3 − 21 0, 34202 + 0, 35837 = 0, 34693. 20 − 21 21 − 20

Chyba E(x) = |f (x) − L1 (x)| ≤

M2 |(x − x0 )(x − x1 )|. (n + 1)!

106

V našem případě je n = 1, f 00 (x) = − sin x a M2 = max | − sin x| = 0, 35837. x∈h20,21i

To znamená, že E(x) ≤

0, 35837 0, 3π 0, 7π ( ) = 1, 1 · 10−5 . 2 180 180

Extrapolací dostaneme 21, 3 − 21 21, 3 − 20 0, 34202 + 0, 35837 = 0, 36328, 20 − 21 21 − 20 0, 35837 0, 3π 1, 7π E(x) ≤ ( ) = 3 · 10−5 . 2 180 180

sin 21◦ 180 ∼

Poznámka Ke zpřesnění hodnoty, kterou jsme před tím numericky spočítali (třeba hodnoty integrálu nebo řešení diferenciální rovnice), je možné využít tzv. Richardsonovy extrapolace. Označme veličinu, kterou počítáme, F a T (h) výsledek numerického výpočtu nějakou metodou. Předpokládejme, že při zjemňování kroku h dostáváme přesnější řešení, lim+ T (h) = F, h→0

a že T (h) lze vyjádřit ve formě polynomu prvního stupně v proměnné hp : T (h) = F + chp + O(hp ). Funkci T (h) nyní nahraďme Lagrangeovým interpolačním polynomem zkonstruovaným v bodech H p , (qH)p , L1 (h) =

h − (H)p h − (qH)p T (H) − T (qH). H p − (qH)p (qH)p − H p

Richardsonovou extrapolací budeme rozumět hodnotu tohoto polynomu v bodě h = 0. 2 (H)p q p T (H) − T (qH) −(qH)p T (H) + T (qH) = = L1 (0) = p H − (qH)p (qH)p − H p qp − 1 T (H) − T (qH) = T (H) + . qp − 1 2

O extrapolaci zde hovoříme proto, že 0 ∈ / (H, qH).

107

Hodnota L1 (0) reprezentuje novou formuli T1 (H) = T (H) +

T (H) − T (qH) , qp − 1

která je lepší aproximací hodnoty F, než bylo T (H) nebo T (qH). Řád chyby aproximace je v tomto případě roven hr , r > p. Zobecněním tohoto postupu je tzv. opakovaná Richardsonova extrapolace. Je možné ji realizovat, pokud numerická metoda T (h) umožňuje vyjádřit chybu aproximace ve formě polynomu s-tého stupně v mocninách hp , tj. ve tvaru T (h) − F ≈ c1 hp + c2 h2p + · · · + cs hsp . Pak můžeme ke konstrukci nové formule použít interpolační polynom s-tého stupně. Jeho hodnota v bodě h = 0 bude představovat přesnější aproximaci veličiny F. K výpočtu hodnoty interpolačního polynomu Ls (0) je možné použít Nevillův algoritmus, který jsme si uvedli na straně 104.

Splajny Mohli bychom se ptát, zda když budeme zvyšovat počet uzlů interpolace, bude interpolační polynom lépe aproximovat funkci f. Odpověď je záporná. Pro jakoukoliv volbu uzlů interpolace vždy existuje spojitá funkce, pro kterou interpolační proces nekonverguje stejnoměrně. Speciálně problematická je interpolace spojitých funkcí při použití ekvidistantních uzlů. (Např. pro funkci y = |x| konverguje interpolační proces s ekvidistantními uzly pouze v bodech x = −1, 0, 1 intervalu h−1, 1i.) Navíc polynomy vyšších stupňů se na koncích více „vlníÿ a to nemusí být v řadě případů pro aproximaci zrovna optimální. Pokud chceme dosáhnout přesnějších výsledků, nepomůže nám tedy vzít větší počet uzlů. Jednou z možností, jak zlepšit interpolaci, je použít k aproximaci po částech polynomiální funkci - tzv. splajn. Splajny byly studovány už začátkem 20. století. Označení „splineÿ použil poprvé v roce 1946 Shoenberg. Původně byla tímto názvem označována speciální pružná dřevěná nebo kovová pravítka, kterých používali konstruktéři, když při vytváření návrhů trupů lodí potřebovali zadanými body proložit hladkou křivku.

108

Definice 5.1 Nechť a = x0 < x1 < x2 < · · · < xn = b, m ∈ N. Funkci s(x) : ha, bi → R, pro kterou platí s ∈ C m−1 ha, bi, s |<xi−1 ,xi > ∈ Pn , n ≤ m, i = 1, 2, . . . , n, budeme nazývat splajnem m-tého stupně. n Množinu všech těchto splajnů označíme Sm . Nyní si ukážeme způsob, jak je možné zkonstruovat kubický splajn. V uzlových bodech x0 , x1 , . . . , xn ∈ ha, bi jsou dány funkční hodnoty f0 , f1 , . . . , fn interpolované funkce. Požadujeme, aby kubický splajn s(x), který chceme najít, splňoval s ∈ C 2 ha, bi, s(x) ∈ Pn , n ≤ 3 na každém subintervalu < xi , xi+1 >, i = 0, 1, . . . , n − 1, s(xi ) = f (xi ), i = 0, 1, . . . , n. Označme hi = xi+1 − xi , si (x) = s(x)|<xi ,xi+1 > . Funkce si (x) uvažujme ve tvaru polynomu 3. stupně si (x) = ai + bi (x − xi ) + ci (x − xi )2 + di (x − xi )3 .

(5.13)

Jeho druhá derivace je lineární funkce s00i (x) = 2ci + 6di (x − xi ).

(5.14)

Předpokládejme na chvíli, že známe s00 (xi ) = Mi . (Mi jsou tzv. momenty.) Funkci s00i (x) pak můžeme na intervalu < xi , xi+1 > nahradit lineárním polynomem x − xi xi+1 − x + Mi+1 . (5.15) s00i (x) = Mi hi hi Když tento vztah integrujme, dostaneme s0i (x) = −Mi si (x) = Mi

(xi+1 − x)2 (x − xi )2 + Mi+1 + Ai , 2hi 2hi

(xi+1 − x)3 (x − xi )3 + Mi+1 + Ai (x − xi ) + Bi . 6hi 6hi 109

(5.16) (5.17)

Protože splajn splňuje v uzlových bodech interpolační podmínky s(xi ) = fi , dostaneme pro x = xi ze vztahu (5.17) Bi = f i −

Mi h2i . 6

(5.18)

Pro x = xi+1 z (5.17) je h3i h2i fi+1 fi − Mi+1 2 − + Mi = Ai = hi 6hi hi 6hi

(5.19)

fi+1 − fi (Mi+1 − Mi )hi − . hi 6 To znamená, že pro x ∈ hxi , xi+1 i, i = 0, . . . , n − 1, je s(x) = si (x) a =

si (x) = Mi

(xi+1 − x)3 (x − xi )3 + Mi+1 + Ai (x − xi ) + Bi . 6hi 6hi

Koeficienty Ai , Bi spočteme pomocí hodnot fi a Mi . Hodnoty n+1 momentů Mi však ve skutečnosti neznáme. Pokusme se je proto nyní spočítat. M0 , Mn si zvolme libovolně např. M0 = f000 , Mn = fn00 . K určení zbývajících neznámých M1 , M2 , . . . , Mn−1 využijeme skutečnosti, že hledaný splajn má mít v uzlových bodech, které leží uvnitř intervalu ha, bi, spojitou první derivaci tj. − 0 s0i−1 (x+ i ) = si (xi ), i = 1, . . . , n − 1.

Pak podle (5.16), (5.19) fi − fi−1 hi−1 hi−1 fi+1 − fi hi hi + Mi + Mi−1 = − Mi − Mi+1 , hi−1 3 6 hi 3 6 hi 6 fi−1 − fi fi − fi−1 hi−1 Mi−1 + 2Mi + Mi+1 = ( − ). hi−1 + hi hi−1 + hi hi−1 + hi hi hi−1 Označíme-li αi =

hi−1 , hi−1 + hi 110

βi =

γi =

hi , hi−1 + hi

fi−1 − fi fi − fi−1 6 ( − ), hi−1 + hi hi hi−1

budou momenty M1 , . . . , Mn řešením soustavy (5.20)      

2 α2 ... 0 0

β1 2 ... 0 0

0 β2 ... 0 0

0 0 ... 0 0

... 0 0 ... 0 0 ... ... ... αn−2 2 βn−2 0 αn−1 2

     

M1 M2 ... Mn−2 Mn−1





    =    

Algoritmus (Konstrukce kubického splajnu) Vstup: n, x0 , . . . , xn , f0 , . . . , fn , f000 , fn00 pro i = 0, . . . , n − 1 hi = xi+1 − xi pro i = 1, . . . , n − 1 hi−1 αi = hi−1 +hi βi = γi =

hi = 1 − αi hi−1 +hi 6 i−1 ( fi−1hi−fi − fih−f ) hi−1 +hi i−1

M0 = f000 , Mn = fn00 M1 , . . . , Mn−1 — řešení soustavy (5.20) pro i = 1, . . . , n − 1 spočteme Mi h2i 6 fi+1 −fi i )hi − (Mi+1 −M hi 6

Bi = f i − Ai =

3

3

−x) i) Výstup: si (x) = Mi (xi+1 + Mi+1 (x−x + Ai (x − xi ) + Bi 6hi 6hi

111

γ1 − α1 f000 γ2 ... γn−2 γn−1 − βn−1 fn00

   .  

Příklad Funkce je zadána tabulkou xi −1 0 1 3 fi 3 3 1 2 Najděte přirozený kubický splajn (tj. takový, pro který jsou momenty M0 = Mn = 0.)

i xi fi hi αi βi γi

0 −1 3 1 − − −

1 0 3 1

2 1 1 2

3 3 2 − 1 1 − 2 3 1 2 − 2 3 0 6 −

M0 = M3 = 0 1 2M1 + M2 = 0 2 1 M1 + 2M2 = 6 3 Odtud: M1 = − Dále

72 18 , M2 = . 23 23

3 61 71 , A1 = − , A2 = , 23 23 46 72 25 B0 = 3, B1 = , B2 = − . 23 23

A0 =

Dostáváme tedy s0 (x) = s1 (x) = −

72 3 3 − (x + 1)3 + x pro x ∈ h−1, 0i 23 23 23

12 61 72 3 (1 − x)3 + x3 − x + pro x ∈ h0, 1i 23 23 23 23 112

s2 (x) =

6 121 71 (3 − x)3 − − x pro x ∈ h1, 3i. 23 46 46

Obrázek: Kubický splajn (Na obrázku je znázorněn kubický splajn plnou čarou a interpolační polynom čárkovaně.)

3.4 3.2 3 2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 –1

0 –0.2

1

x

2

3

n Obecně lze ke generování prostoru Sm použít funkce

3

uk (x) = (x − x0 )k , k = 0, . . . m, vk (x) =(x − xk )m + , k = 1, . . . , n − 1.

(5.21)

Na jednotlivých subintervalech intervalu ha, bi potom máme s(x) =

m X k=0

αk uk (x) +

i−1 X

βk vk (x), x ∈ hx0 , xi i, i = 1, . . . , n.

(5.22)

k=1

Konstruované splajny musí především vyhovovat n+1 interpolačním podmínkám. Ve vztahu (5.22) však potřebujeme určit m + n neznámých koeficientů 3

Zde značíme

½ (x − xk )+ =

x − xk , 0

113

když x − xk ≥ 0, jinak.

α0 , . . . , αm , β0 , . . . , βn−1 . Přidáme proto ještě dalších m − 1 podmínek pro derivace v krajních bodech intervalu. Abychom je mohli do krajních bodů rozdělit rovnoměrně, předpokládejme, že m = 2l − 1. Neznámé koeficienty pak určíme ze soustavy m X k=0 m X

αk uk (xi ) + (i)

αk uk (a) +

k=0 m X

(i)

αk uk (b) +

k=0

n−1 X k=1 n−1 X k=1 n−1 X

βk vk (xi ) = fi , i = 0, . . . , n,

(5.23)

(i)

βk vk (a) = ai , j = 1, . . . , l − 1, (i)

βk vk (b) = bi , i = 1, . . . , l − 1.

k=1

Tato soustava je sice jednoznačně řešitelná, ale kvůli zvoleným bázovým funkcím (5.21), které mají globální nosič, je špatně podmíněná. Z tohoto důvodu bývá zvykem za bázové funkce volit B-splajny ½ 1 pro |x| ≤ 0, 5, B0 (x) = 0 pro |x| > 0, 5. Z Bm+1 (x) =

x+ 21 x− 12

Bm (y) dy, x ∈ R, m = 0, 1, . . .

Např. pro m = 1 máme ½ B1 (x) =

1 − |x| pro |x| ≤ 1, 0 pro |x| > 1.

Pro m = 2 je   2 − (|x| − 12 )2 − (|x| + 12 )3 pro |x| ≤ 12 , (|x| − 32 )2 pro 12 < |x| < 32 , B2 (x) =  0 pro |x| ≥ 32 . Když m = 3, dostaneme  (2 − |x|)2 − 4(1 − |x|)3 pro |x| ≤ 1 1 (2 − |x|)3 pro 1 < |x| < 2 B3 (x) = 6 0 pro |x| ≥ 2 114

Tyto funkce už mají lokální nosič. Dá se ukázat, že když uzly jsou v intervalu , n ≥ 2, m = 2l − 1, pak za bázi ha, bi rozloženy ekvidistantně, h = b−a n n prostoru Sm můžeme vzít B-splajny Bm (

x − xk ), x ∈ ha, bi, k = −l + 1, . . . , n + l − 1. h

Aproximace trigonometrickými polynomy Při každé aproximaci je vhodné pokusit se zachovat vlastnosti interpolované funkce. Z tohoto důvodu nejsou pro aproximaci periodických funkcí výše uvedené interpolační polynomy právě nejvhodnější. Místo prostoru polynomů Pn budeme v takovém případě raději pracovat s prostorem trigonometrických polynomů. Poprvé tento typ interpolace použili Clairout v roce 1759 a Lagrange v roce 1782. Definice 5.2 Trigonometrickým polynomem rozumíme funkci g(x) =

n X

ak cos kx +

n X

bk sin kx,

k=1

k=0

kde n ∈ N, an , bn ∈ C, |an | + |bn | > 0. Množinu všech takových trigonometrických polynomů budeme značit Tn . Poznámka

Prvky prostoru Tn jsou periodické funkce.

Příklad Určete trigonometrický interpolační polynom, když funkce je zadána tabulkou 4π xi 0 2π 3 3 fi 1 0 −1 Řešení: V našem případě interpolační polynom g(x) = a0 + a1 cos x + b1 sin x 115

musí vyhovovat interpolačním podmínkám a0 + a1 = 1, √ 1 3 a0 − a1 + b1 = 0, 2 √2 1 3 a0 − a1 − b1 = −1. 2 2 Odtud a0 = 0, a1 = 1, b1 =

√ 3 , 3

g(x) = cos x +

√ 3 3

sin x.

Předpokládejme, že funkce, kterou aproximujeme, je periodická s periodou 2π, a že interpolujeme na intervalu h0, 2πi. Bez újmy na obecnosti budeme nadále předpokládat, že uzlové body jsou v intervalu h0, 2πi rozloženy ekvidistantně. Věta 5.5 Když jsou dány ekvidistantní uzlové body x0 , . . . , x2n ∈ h0, 2π) a v nich jsou zadány funkční hodnoty f0 , . . . , f2n ∈ R, pak existuje právě jeden interpolační polynom g ∈ Tn tak, že n

a0 X [ak cos kx + bk sin kx], g(x) = + 2 k=1

(5.24)

který vyhovuje interpolačním podmínkám g(

2πj ) = fj , j = 0, . . . , 2n. 2n + 1

Jeho koeficienty jsou určeny vztahy 2n

2 X 2πj ak = fj cos k, k = 0, . . . , n, 2n + 1 j=0 2n + 1

(5.25)

2n

bk =

2 X 2πj k, k = 1, . . . , n. fj sin 2n + 1 j=0 2n + 1

(5.26)

V případě, že počet uzlových bodů je sudý, chybí nám jedna interpolační podmínka k jednoznačnému určení koeficientů trigonometrického polynomu. 116

Ze součtu (5.24) je však možné vyjmout jeden sčítanec, ve kterém figuruje funkce sin nx, protože tato funkce je v uzlových bodech nulová. Věta 5.6 Když je dáno 2n ekvidistantních uzlových bodů x0 , . . . , x2n−1 ∈ h0, 2π) a hodnoty f0 , . . . , f2n−1 ∈ R, pak existuje právě jeden interpolační polynom g ∈ Tn , n−1

a0 X an g(x) = + [ak cos kx + bk sin kx] + cos nx, 2 2 k=1

(5.27)

který vyhovuje interpolačním podmínkám g(

πj ) = fj , j = 0, . . . , 2n − 1. n

Jeho koeficienty jsou určeny vztahy 2n−1 πj 1 X fj cos k, k = 0, . . . , n, ak = n j=0 n

(5.28)

2n−1 1 X πj bk = fj sin k, k = 1, . . . , n − 1. n j=0 n

(5.29)

Příklad Určete trigonometrický interpolační polynom, když je funkce zadána tabulkou 4π xi 0 2π 3 3 fi 1 0 −1 Řešení: Počet ekvidistantních uzlů je lichý. Podle věty 5.5 obdržíme trigonometrický polynom g(x) = a20 + a1 cos x + b1 sin x, kde 2 a0 = (1 · cos 0 + 0 · cos 0 − 1 · cos 0) = 0, 3 2π 4π 2 − 1 · cos ) = 1, a1 = (1 · cos 0 + 0 · cos 3 3 3 √ 2 2π 4π 3 b1 = (1 · sin 0 + 0 · sin − 1 · sin ) = . 3 3 3 3 117

Odtud g(x) = cos x +

√ 3 3

sin x.

Poznámka Uvedená aproximace periodických funkcí trigonometrickými polynomy je diskrétní analogií Fourierových řad. Ve smyslu metody nejmenších čtverců (viz dále str. 124) je tato aproximace nejlepší aproximací na třídě všech polynomů. Poznámka Kromě uvedeného tvaru se často používá vyjádření trigonometrického polynomu jako funkce komplexní proměnné. Díky Eulerovým vzorcům je n

n

a0 X a0 X 1 ikx g(x) = + [ak cos kx + bk sin kx] = + [ak (e + e−ikx )+ 2 2 2 k=1 k=1 n

1 a0 1 X [(ak + ibk )e−ikx + (ak − ibk )eikx ]. + bk (eikx − e−ikx )] = + 2i 2 2 k=1 Když označíme c0 = a20 , ck = 12 (ak − ibk ) a c−k = 12 (ak + ibk ), k = 1, . . . n, lze trigonometrický polynom psát ve tvaru g(x) =

n X

ck eikx .

(5.30)

k=−n

Pokud známe hodnoty funkce f v 2n+1 ekvidistantních uzlech x0 , x1 , . . . x2n , pak pro koeficienty ck platí 2n

1 X −ikxj ck = fj e , k = −n, . . . , n. 2n + 1 j=0

Poznámka Když položíme w = e−ixj , pak ck = P (wk ). Vidíme, že koeficienty ck jsou rovny hodnotám polynomů (2n + 1)-ního stupně v bodech wk , k = 0, . . . , n. To znamená, že jejich výpočet můžeme provést pomocí Hornerova schématu. Koeficienty ck lze také spočítat jiným způsobem, a to pomocí tzv. rychlé Fourierovy transformace (FFT). Princip diskrétní Fourierovy transformace 118

objevil už Gauss. Rychlou realizaci metody zveřejnili v roce 1965 Cooley a Tukey. Předpokládejme, že k dispozici máme N = 2n ekvidistantních uzlů. Označme 2π fj w = e−i N , Fj = . N Pak N X ck = Fj wkj . j=0

Rozdělme tento součet na sudé a liché sčítance tak, že N

ck =

−1 2 X

N

F2j1 (w2 )kj1 + wk

j1 =0

−1 2 X

F2j1 +1 (w2 )kj1 , kde j = 2j1 .

j1 =0

Když k = r N2 + k1 a uvědomíme si, že wN = e−2πi = 1, pak N

(w2 )kj1 = (w2 )r 2 j1 (w2 )k1 j1 = (wN )rj1 (w2 )k1 j1 = (w2 )k1 j1 . Vidíme, že N

N

ck =

−1 2 X j1 =0

F2j1 (w2 )k1 j1 + wk

−1 2 X

F2j1 +1 (w2 )k1 j1 = Pk1 (wk ) + wk P˜k1 (wk ).

j1 =0

Hodnoty Pk1 (wk ), Pk1 (wk ) obsahují výrazy, které lze opět rozdělit na sudé a liché sčítance, atd. O rychlé Fourierově transformaci zde hovoříme proto, že k realizaci algoritmu je řádově zapotřebí pouze O(N log N ) operací, na rozdíl od Hornerova schématu, kde by bylo třeba uskutečnit O(N 2 ) operací.

Hermitova interpolace Pokud kromě funkčních hodnot známe i hodnoty derivací funkce v uzlových bodech, je možné konstruovat interpolační polynom stupně 2n + 1. Věta 5.7 (Hermitova interplace) Když jsou dány body x0 , . . . , xn ∈ ha, bi a hodnoty f0 , . . . , fn , f00 , . . . , fn0 ∈ R, pak existuje právě jeden interpolační 119

polynom p ∈ P2n+1 , který vyhovuje interpolačním podmínkám p(xi ) = fi , p0 (xi ) = fi0 , i = 0, . . . , n. Tento tzv. Hermitův polynom je definován vztahem n X H(x) = [fi Hi0 + fi0 Hi1 ], i=0

kde Hi0 (x) = [1 − 2li0 (xi )(x − xi )][li (x)]2 , Hi1 (x) = (x − xi )[li (x)]2 a li (x) jsou Lagrangeovy polynomy (5.5). Důkaz. Protože Hi0 , Hi1 ∈ P2n+1 , i = 0, . . . , n, je Hermitův polynom H ∈ P2n+1 . Protože 0 Hi0 (xj ) = Hi1 (xj ) = δi,j , 0

Hi0 (xj ) = Hi1 (xj ) = 0, i, j = 0, . . . , n, splňuje Hermitův polynom H interpolační podmínky. Jednoznačnost: Předpokládejme, že existují dva různé interpolační polynomy H1 , H2 ∈ P2n+1 , H1 6= H2 . Pak polynom H = H1 − H2 splňuje v uzlových bodech vztahy H(xj ) = H 0 (xj ) = 0, j = 0, . . . , n. H je polynom stupně 2n + 1 a má n + 1 dvojnásobných kořenů. To ale znamená, že pak je H identicky roven nule a tedy H1 = H2 .

Bézierovy křivky Z hlediska počítačové grafiky není interpolace příliš efektivním nástrojem. Při modelování křivek a ploch je výhodnější, když nepožadujeme, aby křivka či plocha procházela určitými body. Zadané body budeme v takovém případě chápat jako vrcholy polygonu, od kterého se pak odvíjí tvar příslušné křivky popř. plochy. Tento způsob konstrukce poprvé navrhli pro účely automobilového průmyslu v 60. létech B. Bézier (Renault) a P. de Casteljau (Citroen). Za základ aproximace si vzali tzv. Bernsteinovy polynomy,4 které 4

Bernsteinovy polynomy hrály důležitou roli v Bernsteinově důkazu Weierstrassovy věty (viz věta 5.1). Podstata důkazu spočívala v konstrukci posloupnosti funkcí, které stejnoměrně konvergují k dané spojité funkci.

120

tvoří v prostoru Pn bázi. To nám dává možnost dospět k nové reprezentaci křivek a ploch prostřednictvím Bézierových polynomů. Oproti polynomům, které jsme dosud používali, mají Bézierovy polynomy tu přednost, že jejich koeficienty lze geometricky interpretovat, určují hodnotu interpolované křivky v počátečním i koncovém bodě a už dopředu umožňují odhadnout tvar křivky. V dalším výkladu se pro jednoduchost omezíme pouze na interval h0, 1i. To však nebude újmou na obecnosti, protože k obecnému intervalu ha, bi lze . dospět prostřednictvím transformace t = x−a b−a Definice 5.3 Když pro t ∈ h0, 1i označíme µ ¶ n i n t (1 − t)n−i , i = 0, . . . , n, Bi (t) = i

(5.31)

Bernsteinův polynom n-tého stupně, pak p(t) =

n X

bi Bin (t)

(5.32)

i=0

je tzv. Bézierův polynom. Grafem polynomu p(t) je tzv. Bézierova křivka. Parametry bi , i = 0, . . . , n, nazýváme kontrolními (Bézierovými) body. Jejich konvexní obal conv {b1 , b2 , . . . , bn } =

n X i=0

αi bi , αi ≥ 0,

n X

αi = 1,

i=1

je tzv. Bézierův polygon. Příklad Nakreslete Bézierovu křivku, když jsou dány řídicí body a) b0 = (1, 3), b1 = (4, 4), b2 = (5, 2), b) b0 = (1, 3), b1 = (1, 1), b2 = (4, 4), b3 = (5, 2). Řešení: a) p(t) = (1, 3)(1 − t)2 + 2(4, 4)t(1 − t) + (5, 2)t2 , t ∈ h0, 1i. Parametrické rovnice křivky: x(t) = (1 − t)2 + 8t(1 − t) + 5t2 , y(t) = 3(1 − t)2 + 8t(1 − t) + 2t2 , t ∈ h0, 1i. 121

b) p(t) = (1, 3)(1 − t)3 + 3(1, 1)t(1 − t)2 + 3(4, 4)t2 (1 − t) + (5, 2)t3 , t ∈ h0, 1i. Parametrické rovnice křivky: x(t) = (1 − t)3 + 3t(1 − t)2 + 12t2 (1 − t) + 5t3 , y(t) = 3(1 − t)3 + 3t(1 − t)2 + 12t2 (1 − t) + 2t3 , t ∈ h0, 1i.

Poznámka Důležité vlastnosti Bézierových polynomů: 1. V krajních bodech intervalu ha, bi je hodnota Bézierova polynomu rovna hodnotě odpovídajícího kontrolního bodu. 2. Tečny Bézierova polynomu a hrany Bézierova polygonu v krajních bodech splývají. Viz následující obrázek (a). 3. Jestliže Bézierův polygon je konvexní (konkávní), je také konvexní (konkávní) Bézierova křivka. Viz následující obrázek (a). 4. Změna polohy kontrolního bodu nemá jen lokální charakter, ale odrazí se ve tvaru celé křivky. 5. Jednotlivé Bézierovy křivky lze na sebe napojovat. Koncový bod jedné křivky je pak počátečním bodem druhé křivky. Tečny v bodech, kde se dva Bézierovy polynomy stýkají, pak splynou. Viz následující obrázek (b). Obrázek: Vlastnosti Bézierových polynom˚ u 2.stupně.

(b)

(a)

b01

b1

b0 b0

b2 = b00

b2 b1 122

b02

Určení hodnot Bézierova polynomu p(t) pro t ∈ h0, 1i se provádí Casteljauovým algoritmem, který je založen na tom, že Bézierovu křivku je možné zkonstruovat jako konvexní kombinaci řídících bodů. Např. když v rovině jsou dány tři kontrolní body b0 , b1 , b2 , pak p(t) =

2 X

bi Bi2 (t) = (1 − t2 )2 b0 + 2t(1 − t)b1 + t2 b2 =

i=0

= (1 − t) [(1 − t)b0 + tb1 ] + [(1 − t)b1 + tb2 ] . | {z } | {z } b10

b11

To znamená, že když zvolíme t0 ∈ h0, 1i pevně, můžeme určit bod p(t0 ), který leží na Bézierově křivce tak, že položíme b10 = (1 − t0 )b0 + t0 b1 a b11 = (1 − t0 )b1 + t0 b2 a spočteme hodnotu p(t0 ) = (1 − t0 )b10 + t0 b11 . Obrázek: Konstrukce bodu na Bézierově křivce n = 2, t = 21 . b1

b10

b20

b0

b11

b2

Obecně, když je dáno n+1 řídících bodů, položíme b00 = b0 , b01 = b1 , . . . , b0n = bn a postupně pro k = 0, 1, . . . , n a j = 0, 1, . . . , n − k počítáme subpolynomy (t) + tbk−1 bkj (t) = (1 − t)bk−1 j j+1 (t). Pro Bézierův polynom pak platí p(t) = bn0 (t). 123

Poznamenejme ještě, že výpočet koeficientů bkj (t) můžeme zachytit v tabulce: b0 b1 .. .

b10 .. .

bn−1 b1n−2 bn

b1n−1

.. . .. . .. .

.. .

.. .

b0n−1 b1n−1

bn0

Algoritmus (Casteljauův algoritmus) Vstup: n, b0 , . . . , bn , t pro i = 0, 1, . . . , n b0i = bi pro k = 1, . . . , n pro i = 0, 1, . . . , n − k bki = (1 − t)bk−1 + tbk−1 i i+1 Výstup: p(t) = bn0 (t) Bézierovy polynomy rychle konvergují k hledané křivce, a proto se využívají při kreslení křivek na počítači.

Metoda nejmenších čtverců Jinou možností, když netrváme na tom, aby aproximující funkce procházela zadanými body, je pokusit se najít funkci, která je těmto hodnotám v jistém smyslu co nejblíž. Pokud velikost odchylky aproximující a dané funkce budeme měřit v normě l2 (a, b), pak hovoříme o aproximaci metodou nejmenších čtverců. Základní myšlenka metody spočívá v tom, že si zvolíme základní 124

funkce ϕ1 , ϕ2 , . . . , ϕm ,

5

aproximující funkci hledáme ve tvaru ϕ(x) =

m X

cj ϕj (x)

(5.33)

j=0

a snažíme se najít konstanty c1 , . . . , cm tak, aby číslo n n m X X X 2 S(c1 , c2 , . . . , cm ) = [f (xi ) − ϕ(xi )] = [f (xi ) − cj ϕj (xi )]2 i=0

i=0

j=0

bylo minimální. Vektor c1 , . . . , cn je stacionárním bodem S, když n

m

X X ∂S =2 [f (xi ) − cj ϕj (xi )]ϕk (xi ) = 0, k = 1, . . . , m. ∂ck i=0 j=0 Dostáváme tak soustavu tzv. normálních rovnic m X j=0

cj

n X

ϕj (xi )ϕk (xi ) =

i=0

n X

f (xi )ϕk (xi ), k = 1, . . . , m.

(5.34)

i=0

Když její řešení c1 , . . . , cn dosadíme do vztahu (5.33) obdržíme funkci ϕ(x). Poznámka

Soustavu (5.34) lze pomocí skalárního součinu v l2 psát jako m X (ϕj , ϕk )cj = (f, ϕk ), k = 0, . . . , n.

(5.35)

j=0

Poznámka Pokud zvolíme ϕj (xi ) = xji , j = 1, . . . , m, znamená to, Pmže funkci f aproximujeme polynomem m-tého stupně. Pak funkce ϕ(x) = j=0 cj xji a soustava normálních rovnic (5.34) přejde na tvar n n m X X X j+k xi )cj = f (xi )xki , k = 0, . . . , n. ( j=0 i=0

5

Za základní funkce lze {1, cos x, sin x, cos 2x, sin 2x, . . . , sin mx}

i=0

zvolit

125

např.

{1, x, x2 , . . . , xm }

nebo

Algoritmus (Metoda nejmenších čtverců) Vstup: n, x0 , . . . , xn , f0 , . . . , fn , m, ϕ1 (x), . . . , ϕm (x) Pro j = 0, . . . , m pro k = 0, . . . , m − k P ϕj,k = ni=0 ϕj (xi )ϕk (xi ) ϕk,j = ϕj,k P fk = ni=0 fi ϕj (xi ). c0 , c1 , . . . , cm — řešení soustavy 

ϕ0,0 ϕ0,1  ϕ1,0 ϕ1,1   ... ... ϕm,0 ϕm,1 Výstup: ϕ(x) =

P j

 . . . ϕ0,m c0  c1 . . . ϕ1,m   ... ...  ... . . . ϕm,m cm





 f0   f1  =    ... . fm

cj ϕj (x).

Příklad Funkce je zadána tabulkou xi 1 2 3 5 fi 3 3 1 2 Pomocí metody nejmenších čtverců aproximujte tuto funkci kvadratickým polynomem. Řešení: Funkci aproximujeme polynomem ϕ(x) = c0 + c1 x + c2 x2 , jehož koeficienty získáme jako řešení normální soustavy c0 (ϕ0 , ϕ0 ) + c1 (ϕ1 , ϕ0 ) + c1 (ϕ2 , ϕ0 ) = (f, ϕ0 ), c0 (ϕ0 , ϕ1 ) + c1 (ϕ1 , ϕ1 ) + c1 (ϕ2 , ϕ1 ) = (f, ϕ1 ), c0 (ϕ0 , ϕ2 ) + c1 (ϕ1 , ϕ2 ) + c1 (ϕ2 , ϕ2 ) = (f, ϕ2 ), 126

kde ϕ0 = 1, ϕ1 = x, ϕ2 = x2 . Je tedy (ϕ0 , ϕ0 ) =

3 X

ϕ0 (xi )ϕ0 (xi ) =

i=0

3 X

1 = 4,

i=0

(ϕ0 , ϕ1 ) = (ϕ1 , ϕ0 ) = (ϕ0 , ϕ2 ) = (ϕ2 , ϕ0 ) =

3 X i=0 3 X

ϕ1 (xi )ϕ0 (xi ) = ϕ2 (xi )ϕ0 (xi ) =

3 X i=0 3 X

i=0

(ϕ1 , ϕ1 ) =

3 X

xi = 11, x2i = 39,

i=0

ϕ1 (xi )ϕ1 (xi ) =

i=0

3 X

x2i = 39,

i=0

(ϕ1 , ϕ2 ) = (ϕ2 , ϕ1 ) =

3 X

ϕ2 (xi ), ϕ1 (xi ) =

i=0

(ϕ2 , ϕ2 ) = (f, ϕ0 ) = (f, ϕ1 ) = (f, ϕ1 ) =

3 X i=0 3 X i=0 3 X i=0 3 X

3 X

x3i = 161,

i=0

ϕ2 (xi )ϕ2 (xi ) =

3 X

x4i = 723,

i=0

f (xi )ϕ0 (xi ) = 9, f (xi )ϕ1 (xi ) = 22, f (xi )ϕ1 (xi ) = 74.

i=0

Řešíme tedy soustavu 4c0 + 11c1 + 39c2 = 9, 11c0 + 39c1 + 161c2 = 22, 39c0 + 161c1 + 723c2 = 74. 49 Odtud c2 = 14 , c1 = − 37 a c0 = 10 . Metodou nejmenších čtverců jsme získali 20 aproximaci zadané funkce ve tvaru

37 49 1 ϕ(x) = x2 − x + . 4 20 10 127

Obrázek: Aproximace metodou nejmenších čtverců 6

ϕ(x) = 14 x2 −

5

4

y3

2

1

0

1

2

3 x

128

4

5

6

37 x 20

+

49 . 10

Kapitola 6 Numerická kvadratura a derivace

Numerickou metodu pro integrování neboli numerickou kvadraturu použijeme v případě, když integrál nelze určit analyticky prostřednictvím elementárních funkcí nebo když integrovaná funkce je zadaná tabulkou. Princip numerické kvadratury spočívá v tom, že funkci f (x) nahradíme jinou (jednodušší) funkcí ϕ a tu pak analyticky integrujeme. Pokud je ϕ dobrou aproRb ximací funkce f na nějakém intervalu ha, bi, je také a ϕ(x) dx dobrou aproRb ximací a f (x) dx. Když x0 , x1 , . . . , xn ∈ ha, bi a známe hodnoty integrované funkce v uzlových bodech, můžeme ji nahradit třeba Lagrangeovým polynomem P f (x) ≈ Ln (x) = ni=0 f (xi )li (x). Integrací obdržíme Rb Rb RbP Rb P f (x) dx ≈ a Ln (x) dx = a ni=0 f (xi )li (x) dx = ni=0 f (xi )( a li (x) dx). a Rb Když označíme wi = a li (x) dx, obdržíme kvadraturní formuli Z b n X f (x) dx ≈ wi f (xi ). (6.1) a

i=0

Koeficienty wi v kvadraturní formuli budeme nazývat váhy. Dále bude pro 129

chybu aproximace platit Z

Z

b

|

f (x) dx − a

Z

b

b

Ln (x) dx| ≤ a

|f (x) − Ln (x)| dx.

(6.2)

a

Vidíme, že chyba numerické kvadratury závisí na tom, jak velká je chyba aproximace funkce. Poznámka V dalším textu budeme používat pro zadaný integrál označení Z

b

I=

f (x) dx

(6.3)

a

a pro kvadraturní formuli Qn =

n X

wi fi .

(6.4)

i=0

Poznamenejme ještě, že integrál i kvadratura představují lineární operátory I, Qn : C 0 ha, bi → R.

Definice 6.1 Pokud uzlové body volíme tak, že a = x0 a b = xn , hovoříme o uzavřených kvadraturních formulích, když a < x0 < x1 < · · · < xn < b, x0 − a = b − xn , pak konstruujeme otevřené kvadraturní fomule. Budeme se zabývat vesměs uzavřenými formulemi.

Newtonovy-Cotesovy kvadraturní vzorce Nejjednodušším typem kvadraturních formulí jsou Newtonovy-Cotesovy vzorce. Předpokládejme, že známe hodnoty integrované funkce v ekvidistantně rozložených uzlech. Označme tyto uzly x0 , x0 + h, . . . , x0 + nh ∈ ha, bi, h = 130

b−a . n

Příklad Pro n = 1 zkonstruujte Newtonovu-Cotesovu kvadraturní formuli a odhadněte chybu kvadraturního vzorce. , x0 = a, x1 = b. Integrovanou funkci f aproximuŘešení: Položíme h = b−a 2 jeme lineárním polynomem L1 (x) = x−b f (a) + x−a f (b). Pak a−b b−a Z

Z

b

b

f (x) dx ≈ a

L1 (x) dx = a

b−a (f (a) + f (b)). 2

Nyní odhadneme chybu. Lze psát Z b Z b f (x) − L1 (x) b−a (f (a) + f (b)) = (x − a)(x − b) dx. f (x) dx − 2 (x − a)(x − b) a a Pokud je funkce f ∈ C 2 (a, b), pak z věty o střední hodnotě máme: existuje takové y ∈ ha, bi, že Z b Z b b−a f (y) − L1 (y) f (x) dx − (f (a) + f (b)) = (x − a)(x − b) dx. 2 (y − a)(y − b) a a Odtud a z odhadu chyby (5.8) Z b b−a f 00 (ξ) b − a f (x) dx − (f (a) + f (b)) = (− ). 2 2 6 a Poznámka

(6.5)

Kvadraturní formule Q1 =

b−a (f (a) + f (b)), 2

(6.6)

kterou jsme odvodili v předchozím příkladu, je vzhledem ke své geometrické interpretaci nazývána lichoběžníkové pravidlo. Definice 6.2 Řekneme, že kvadraturní vzorec je řádu n, když kvadraturní formule je přesná pro polynomy p ∈ Pk , k ≤ n. Příklad Za předpokladu, že uzly jsou rozloženy ekvidistantně, najděte na intervalu ha, bi kvadraturní formuli prvního a druhého řádu. Řešení: Podle vztahu (6.4) je kvadraturní formule 1. řádu Q1 = w0 f0 + w1 f1 . Stejně jako v předchozím příkladu položme x0 = a, x1 = b. Kvadraturní 131

formule bude řádu 1, když bude dávat přesné hodnoty pro polynomy p0 (x) = 1 a p1 (x) = x. To znamená (viz (6.1)) Z b w0 + w1 = 1 dx = b − a, a Z b b2 − a 2 −w0 a + w1 b = . x dx = 2 a Rb Odtud w0 = w1 = b−a a a f (x) dx ≈ b−a (f0 + f1 ) = Q1 (x). 2 2 Ukázali jsme, že lichoběžníkové pravidlo představuje kvadraturní formuli 1. řádu. Podle (6.4) kvadraturní formule 2. řádu Q2 = w0 f0 + w1 f1 + w2 f2 . Zvolme x0 = a, x1 = a+b , x2 = b. Z požadavku, aby kvadraturní formule byla přesná 2 pro polynomy p0 (x) = 1, p1 (x) = x, p2 (x) = x2 , získáme soustavu Z b w0 + w1 + w2 = 1 dx = b − a, a Z b b2 − a2 a+b x dx = w0 a + w1 + w2 b = , 2 2 a Z b b3 − a3 (a + b)2 2 2 x2 dx = + w2 b = . w0 a + w1 4 3 a Odtud a0 = a2 =

b−a , a1 3

=

b−a . 3

Obdrželi jsme integrální formuli 2. řádu

b−a a+b (f (a) + 4f ( ) + f (b)) = Q2 (x), 3 2 tzv. Simpsonovo pravidlo. Q2 (x) =

(6.7)

Poznámka Obdobně lze odvodit hodnoty wi i pro ostatní NewtonovyCotesovy vzorce. Když interval, přes který integrujeme, je h−1, 1i, pak pro koeficienty v Newtonových-Cotesových vzorcích dostaneme n 1 Lichoběžníkové pravidlo 2 Simpsonovo pravidlo 3 Newtonovo 3/8 pravidlo 4 Milnovo pravidlo

132

wi 1 2 1 3 3 8 14 45

1 2 4 3 9 8 64 45

1 3 9 8 24 45

3 8 64 45

14 45

Složené kvadraturní vzorce Položme si nyní otázku, zda lze nějakým způsobem dosáhnout zvýšení přesnosti kvadraturního vzorce. Přidáním uzlových bod˚ u sice zvýšíme algebraický řád kvadraturní formule, ale ani pak ještě nemusíme dosáhnout požadované přesnosti. Východiskem z tohoto problému je použití složených kvadraturních vzorc˚ u. Postupujeme tak, že interval ha, bi rozdělíme na subintervaly < aj , aj+1 >, j = 0, . . . , m − 1, a na každém z nich aplikujeme kvadraturní formuli (6.4). Když j-tý interval rozdělíme na nj dílčích intervalů, dostaneme Z

b

f (x) dx = a

Z

m−1 X

aj+1

(

f (x) dx) ≈ aj

j=0

kde

Z

nj m−1 XX

wij fi ,

(6.8)

j=0 i=0

aj+1

wij =

lij (x) dx. aj

Tímto způsobem můžeme získat následující složené Newtonovy-Cotesovy vzorce: a) Obdélníkové pravidlo. Když na j-tém intervalu hxj , xj+1 i nahradíme funkci f (x) hodnotou f (xj + h2 ), dospějeme ke kvadraturní formuli R(f, h) =

n−1 Z X j=0

xj+1

f (x) dx ≈ h xj

n−1 X

h f (xj + ). 2 j=0

Chyba aproximace je Z b Z b h2 | f (x) dx − R(f, h) dx| ≤ (b − a) f 00 (ξ), ξ ∈ ha, bi. 24 a a Algoritmus: (Obdélníkové pravidlo) Vstup: a, b, f (x), n h = b−a n pro j = 0, 1, . . . , n xj := a + jh P h R(f, h) = h n−1 j=0 f (xj + 2 ) Výstup: R(f, h) 133

(6.9)

(6.10)

b) Lichoběžníkové pravidlo. Pokud na intervalu hxj , xj+1 i aproximujeme funkci f (x) aritmetickým pr˚ uměrem 12 (f (xj ) + f (xj+1 )), pak získáme kvadraturní formuli T (f, h) =

n−1 Z X

f (x) dx ≈ h xj

j=1

=

xj+1

n−1 X 1 j=0

2

(f (xj ) + f (xj+1 )) =

h (f0 + 2f1 + 2f2 + · · · + 2fn−1 + fn ). 2

(6.11)

Chyba aproximace je zde Z

Z

b

|

b

f (x) dx − a

T (f, h) dx| ≤ (b − a) a

h2 00 f (ξ), ξ ∈ ha, bi. 12

(6.12)

Algoritmus: (Lichoběžníkové pravidlo) Vstup: a, b, f (x), n h = b−a n pro j = 0, 1, . . . , n xj := a + jh P T (f, h) = h2 n−1 j=0 (f (xj ) + f (xj+1 )) Výstup: T (f, h) c) Simpsonovo pravidlo. Další kvadraturní vzorec obdržíme, když na intervalu hxj , xj+2 i funkci f (x) nahradíme kvadratickým polynomem f (xj )

(x − xj+1 )(x − xj+2 ) (x − xj )(x − xj+2 ) + f (xj+1 ) + (xj − xj+1 )(xj − xj+2 ) (xj+1 − xj )(xj+1 − xj+2 ) +f (xj+2 )

(x − xj )(x − xj+1 ) . (xj+2 − xj )(xj+2 − xj+1 )

Dostaneme S(f, h) =

m−1 X Z xj+1 j=0

xj−1

f (x) dx ≈

h (f0 + 4f1 + 2f2 + 4f3 + · · · + 2fn−2 + 4fn−1 + fn ). 3 (6.13) 134

Chyba aproximace Z

Z

b

|

b

f (x) dx −

S(f, h) dx| ≤ (b − a)

a

a

h4 (4) f (ξ), ξ ∈ ha, bi. 180

(6.14)

Algoritmus: (Simpsonovo pravidlo) Vstup: a, b, f (x), n h = b−a n pro j = 0, 1, . . . , n xj := a + jh S(f, h) =

h 3

n −1 2

P

(f (x2j ) + 4f (x2j+1 ) + f (x2j+2 ))

j=0

Výstup: S(f, h) Příklad Simpsonovým pravidlem spočítejte

R2

ex 1 x

dx.

Řešení: Pro n = 6 obdržíme Z 2 x e 4 7 1 4 4 3 1 5 4 11 1 . 1 dx = e+ e 6 + e 3 + e 2 + e 3 + e 6 + e2 = 3, 059142296, 18 21 12 27 15 33 36 1 x Z

2

| 1

ex 1 1 dx − S(f, h)| ≤ 15 = 0, 0052083. x 180 6

Ve vzorcích (6.10), (6.12) a (6.14) se vyskytují hodnoty 2. a 4. derivace funkce f. Když však je integrovaná funkce zadaná tabulkou, nemůžeme hodnoty jejích derivací stanovit přesně. Můžeme jen získat jejich hrubý odhad. To znamená, že také dostáváme jen hrubý odhad chyby. Z tohoto důvodu se k odhadu velikosti chyby používá Rungeova metoda dvojnásobného kroku. Spočívá v tom, že z hodnot kvadraturních formulí pro n a 2n intervalů provedeme odhad prostřednicvím vztahů Z

Z

b

|

b

f (x) dx − a

a

1 R2n (f, h) dx| ≤ |Rn (f, h) − R2n (f, h)|, 3 135

Z | Z

Z

b

b

f (x) dx − a b

|

Z

a b

f (x) dx − a

1 T2n (f, h) dx| ≤ |Tn (f, h) − T2n (f, h)|, 3

S2n (f, h) dx| ≤ a

1 |Sn (f, h) − S2n (f, h)|. 15

R2 x Příklad Lichoběžníkovým pravidlem spočítejte 1 ex dx. (Volte n = 6.) Odhad chyby kvadratury proveďte pomocí vztahu (6.12) i Rungeovou metodou dvojnásobného kroku. Řešení: Pro n = 6 je Z 2 x e 1 1 1 7 1 4 1 3 1 5 1 11 dx ≈ T6 = e+ e2 + e 6 + e 3 + e 2 + e 3 + e 6 = 3, 063385879, 12 24 7 8 9 10 11 1 x Z 2 x e h3 1 1 | dx − T6 (f, h)| ≤ M2 = ( )3 2, 718 = 0, 10487 · 10−2 . 12 6 12 1 x Pro n = 3 dostaneme Z 2 x e 1 1 4 1 5 1 dx ≈ T3 = e + e 3 + e 3 + e2 = 3, 076116630, 6 4 5 12 1 x Z 2 x e 1 | dx − T6 (f, h)| ≤ |T3 (f, h) − T6 (f, h)| = 0, 4243583667 · 10−2 . 3 1 x Poznámka V roce 1955 navrhl Romberg způsob, který vede ke zpřesnění kvadraturních formulí nižších řádů. Na tuto myšlenku ho přivedla následující skutečnost: Pokud funkce f má derivace až do (m + 1)-ního řádu 1 1 a B2 = 16 , B4 = − 30 , B6 = 42 , . . . , Bm jsou Bernoulliova čísla, lze chybu lichoběžníkového pravidla psát ve tvaru Z b m X Bs 2s (2s−1) T (f, h) − f (x) dx = h [f (b) − f (2s−1) (a)] 2s! a s=0 (tzv. Eulerův-Maclaurinův vzorec). Chybu složené lichoběžníkové formule je tedy možné rozvinout do řady sudých mocnin integračního kroku h Z b T (f, h) − f (x) dx = a1 h2 + a2 h4 + · · · + am h2m + O(h2m+2 ). (6.15) a

136

To znamená (viz poznámku na str. 107), že je možné aplikovat Richardsonovu extrapolaci. , lichoběžníkové Odvození nových formulí: Když funkce f ∈ C 2 ha, bi, h = b−a n pravidlo pro krok h označíme T00 = T (h, f ) a pro krok h2 označíme T10 = T ( h2 , f ), pak podle vztahu (6.15) platí Z Z

b a b a

f (x) dx = T00 + a1 h2 + O(h4 ), f (x) dx = T10 + a1

h2 + O(h4 ). 4

Když z předchozích dvou výrazů vyloučíme a1 , dostaneme Z

b a

1 f (x) dx = (4T10 − T00 ) + O(h4 ). 3

Dále tuto lineární kombinaci složených lichoběžníkových pravidel budeme značit 1 T11 = (4T10 − T00 ). 3 Obdobně pro f ∈ C 4 ha, bi obdržíme Z Z

b a

f (x) dx = T11 + a2 h4 + O(h6 ),

b

f (x) dx = a

T21

h4 + b2 2 + O(h6 ). 4

Opět eliminujeme a2 a dostaneme novou kvadraturní formuli T22 =

1 (16T21 − T11 ), 15

která je ale už přesnější než předchozí formule T11 , protože aproximuje hodnotu uvažovaného integrálu s chybou řádu h6 . A tak bychom mohli při dostatečné hladkosti funkce f pokračovat dále. Obecně je pro m = 1, 2, . . . k Tkm+1 =

m 4m Tk+1 − Tkm−1 . 4m − 1

137

Pokud f ∈ C 2m ha, bi, platí pro odhad chyby Z

b

| a

f (x) dx − Tkm | ≤ c(m)kf (2m) k∞ (

h 2m ) , k = 0, 1, . . . 2k

Protože Rombergova kvadratura představuje v podstatě aplikaci Richardsonovy extrapolace na lichoběžníkové pravidlo, lze výpočet hodnoty Tkm při ručním počítání provést v tabulce k 2k m = 0 m = 1 m = 2 0 1 T00 1 2

T10

T11

2 4

T20

T21

T22

Hodnoty T11 , T22 , . . . na diagonále tabulky představují aproximace hodnoty integrálu s chybami řádu O(h4 ), O(h6 ), . . . Podle nich hned zjistíme, kdy jsme dosáhli požadované přesnosti. Tak se během výpočtu můžeme rozhodnout, zda už výpočet zastavíme nebo zda ještě přidáme a dopočteme další řádky tabulky. Pro výpočet použijeme hodnoty už předtím spočtené. Algoritmus (Rombergova kvadratura) Vstup: f (x), h, m pro i = 1, 2, . . . , m Ti0 = T (f, 2hi ) pro j = i, i + 1, . . . , m Tij = Tij−1 +

j−1 Tij−1 −Ti−1 j 4 −1

Výstup: I ≈ Tmm Příklad Spočítejte

R1

ln(x+1) 0 (x2 +1)

dx. 138

Řešení: i 0 1 2 3 4 5

Ti0 0.173287 0.248829 0.256458 0.270769 0.271841 0, 272109

Ti1 0.274010 0.272334 0.272206 0.272199 0.272198

Ti2

Ti3

Ti4

Ti5

0.272222 0.272197 0.272197 0.272198 0.272198 0.272198 0.272198 0.272198 0.272198 0.272198

Přibližná hodnota integrálu je 0,272198.

Poznámka Rombergova metoda je účinným nástrojem v případě, když integrovaná funkce je dost hladká. Není to však metoda univerzální. Pro integrály z periodických funkcí se ukazuje jako postačující aplikovat lichoběžníkové pravidlo a podobná situace je i při stanovení hodnoty nevlastního integrálu.

Gaussovy kvadraturní vzorce Rb Uvažujme nyní integrál v obecnějším tvaru a ω(x)f (x) dx, kde ω je tzv. vá1 hová P funkce . Opět budeme konstruovat kvadraturní formuli ve formě součtu i wi f (xi ). Tentokrát ale uzlové body xi nebudou zvoleny ekvidistantně, ale (stejně jako i koeficienty wi ) budou závislé na tvaru váhové funkce ω. Přesněji řečeno, za uzlové body zvolíme kořeny polynomů, které jsou ortogonální s vahou ω na intervalu ha, bi. 2 Výhodou Gaussových kvadratur je, že při stejném počtu uzlů dávají přesnější výsledky než Newtonovy-Cotesovy vzorce. 1 intervalu spojitá a kladná. Např. ω(x) = 1, ω(x) = √ Tato funkce je na uvažovaném 1 1 − x2 nebo ω(x) = √1−x . 2 2 Viz definice (1.13). Příslušný skalární součin polynomů pi (x), pj (x) má tvar

Z (pi , pj ) =

b

ω(x)pi (x)pj (x)dx. a

139

Pro ω(x) = 1 a ha, bi = h−1, 1i dostáváme tzv. Gaussův-Legendrův kvadraturní vzorec Z 1 n X f (x) dx = wi f (xi ). (6.16) −1

i=0

Uzlové body xi zde volíme jako kořeny Legendrových polynomů Pn (x) =

1 dn (x2 − 1)n 2n · n! dxn

a koeficienty wi jsou definovány vztahy wi =

−2 . 0 (n + 2)Pn+2 (xi )Pn+1 (xi )

Hodnoty uzlových bodů a koeficientů v Gaussově-Legendrově vzorci nemusíme počítat, protože jsou tabelovány v literatuře. Uveďme si alespoň některé z nich: n xi wi 1 √ 1 ± 3 1 8 2 0 9 q 5 ± 35 9 3 ±0, 861136311 0, 347854845 ±0, 339981043 0, 652145154

Poznámka Gaussovy-Legendrovy kvadraturní vzorce pro R1 n = 0 : −1 f (x) dx = f (0) + 13 f 00 (ξ), ξ ∈ (−1, 1). R1 1 n = 1 : −1 f (x) dx = f (− √13 ) + f ( √13 ) + 135 f (4) (ξ), ξ ∈ (−1, 1). q q R1 3 8 5 1 5 f (6) (ξ), n = 2 : −1 f (x) dx = 9 f (− 5 ) + 9 f (0) + 9 f ( 35 ) + 15750 ξ ∈ (−1, 1). Poznámka Když pro x ∈ h−1, 1i je |f (2n+2) (x)| ≤ M a označíme dn+1 =

22n+3 [(n + 1)!]4 , (2n + 3)[2(n + 1)!]3 140

lze chybu v případě jednoduché kvadratury (6.16) odhadnout vztahem Z

1

|

f (x) dx − −1

k X

wi f (xi )| ≤ M dn+1 .

i=0

Poznámka Jestliže chceme integrovat přes interval ha, bi a stanovit hodnotu Rb integrálu a f (x) dx, převedeme jej substitucí x = 21 (t(b − a) + a + b) (dx = R 2x−a−b b−a 1 1 (b − a) dt, t = ), na integrál f ( 12 ((b − a)t + a + b)) dt. Pro jeho 2 b−a 2 −1 výpočet lze použít buď jednoduchý Gaussův-Legendrův vzorec Z b Z b−a 1 1 f (x) dx = f ( (b − a)t + a + b)) dt ≈ 2 2 a −1 n

b−aX 1 ≈ wi f ( (xi (b − a) + a + b)), 2 i=0 2 kde chyba Z

b

f (x) dx −

X

a

f (2n) (ξ) wi f (xi ) = (2n + 2)!

Z

1

qn+1 (x)2 dx.

−1

Eventuálně je možné použít vzorec složený Z

b

f (x) dx = a

ve kterém h = O(h2m ).

m−1 X Z aj+1 aj

j=0 b−a , m

m−1 n hX X 1 1 f (x) dx ≈ ( wi f ( xi h + a + h(j + )), 2 j=0 i=0 2 2

xi = a + hi. Chyba tohoto kvadraturního vzorce je řádu

Algoritmus: (Složený Gaussův-Legendrův kvadraturní vzorec) Vstup: a, b, f (x), m, n, w0 , . . . , wn h = b−a m pro j = 0, 1, . . . , n xj = a + hj. P Pn 1 1 G(f, h) = h2 m i=0 wi f ( 2 xi h + a + h(j + 2 )). j=0 Výstup: I ≈ G(f, h) 141

Příklad Pomocí kvadraturního vzorce stanovte hodR 2 eGaussova-Legendrova x notu integrálu 1 x dx. (Volte m = 2, n = 1.) Řešení: Integrační interval rozdělíme na dva subintervaly a na každém z nich integrál nahradíme Gaussovým-Legendrovým kvadraturním vzorcem 1. řádu. Z

2 1

1

1

1

1

5 + √3 5 − √3 7 + √3 7 − √3 ex . 1 dx = (f ( )+f ( )+f ( )+f ( )) = 3, 059035425. x 4 4 4 4 4

Numerická derivace Cílem numerického derivování je určení hodnoty derivace funkce v určitém bodě z hodnot funkce v bodech z jeho okolí. Nejprve se pokusme najít vzorec pro numerickou derivaci za předpokladu, že funkce f ∈ C (5) hx − h, x + hi. Protože jsou splněny předpoklady Taylorovy věty, můžeme psát f (x + h) = f (x) + hf 0 (x) +

h2 00 h3 h4 f (x) + f 000 (x) + f (4) (x) + O(h5 ), 2! 3! 4!

f (x − h) = f (x) − hf 0 (x) +

h3 h4 h2 00 f (x) − f 000 (x) + f (4) (x) + O(h5 ). 2! 3! 4!

Odtud

f (x + h) − f (x − h) h2 000 0 = f (x) + f (x) + O(h4 ). 2h 3! Derivaci funkce f můžeme tedy přibližně spočítat pomocí vzorce f 0 (x) ≈

f (x + h) − f (x − h) . 2h

Obecný tvar vzorců pro numerickou derivaci můžeme odvodit také tak, že funkci aproximujeme interpolačním polynomem a pak derivujeme. Nevýhodou je, že i když aproximace uvažované funkce je dobrá, může být chyba 142

spočítané numerické derivace příliš velká.

Když vyjdeme z Lagrangeova interpolačního polynomu pro ekvidistantní uzly (viz vztahy (5.4), (5.8) a označení na str. 98), pak f (x) =

n X i=0

Zde t =

x−xi ,i h

f

(k)

li (t)fi +

hn+1 qn+1 (t) (n+1) f (ξ). (n + 1)!

= 1, . . . , n. Pro k-tou derivaci dostaneme

m 1 X (k) hn+1−k qn+1−k (t) (n+1) (x) = k li (t)fi + f (ξ), h i=0 (n + 1)!

(6.17)

Zatímco u numerické interpolace je chyba metody úměrná kroku h a zaokrouhlovací chyba na něm nezávisí (viz poznámku na str. 100), je u numerické derivace chyba metody úměrná kroku h a zaokrouhlovací chyba je nepřímo úměrná kroku h. To znamená, že když krok h je příliš velký, je chyba metody pro numerickou derivaci velká. Když je h malé, je velká zaokrouhlovací chyba. Je však možné spočítat optimální krok hopt (viz. kniha [8]) tak, aby zaokrouhlovací chyba a chyba metody byly co nejmenší. Poznámka Na závěr si uveďme některé vzorce, které bývají nejčastěji používány pro výpočet 1. a 2. derivace. Získáme je derivováním interpolačních polynomů 1. nebo 2. stupně. Když známe hodnoty funkce f v ekvidistantních uzlech x − h, x, x + h, pak je možné spočítat f 0 (x) =

f (x + h) − f (x) h 00 − f (ξ), h 2

f (x) − f (x − h) h 00 + f (ξ), h 2 f (x + h) − f (x − h) h2 000 f 0 (x) = − f (ξ), 2h 6 f (x + h) − 2f (x) + f (x − h) − h2 f 000 (ξ). f 00 (x) = h2 f 0 (x) =

143

Příklad Stanovte přibližně f 0 (0, 5), když f (x) = Řešení:

ex x

a h = 0, 1.

xi 0.4 0.5 0.6 fi 3, 729561745 3, 297442542 3, 036864667

Nejprve spočteme derivaci z hodnoty funkce v bodě x a jednom jeho sousedním bodě f (x + h) − f (x) , h . f (0, 6) − f (0, 5) f 0 (0, 5) = = −2, 60577875. 0.1 f 0 (x) =

Když derivaci určíme z hodnot funkce ve dvou bodech, které sousedí s bodem x, máme f (x + h) − f (x − h) , 2h . f (0.6) − f (0.4) f 0 (0.5) = = −3, 46348539. 0.2 f 0 (x) =

Skutečná hodnota derivace f 0 (0.5) = −3, 297442542. Druhý vzorec nám dal přesnější výsledek.

144

Literatura [1] Buchanan J.I.,Turner P.R.: Numerical Methods and Analysis, New York 1992 [2] Dont M.: Numerické metody - cvičení, ČVUT Praha 1990 [3] Míka S.: Numerické metody algebry, SNTL Praha 1985 [4] Mošová V.: Matematická analýza I, UP Olomouc 2002 [5] Kobza J.: Numerické metody, UP Olomouc 1984, 1993 [6] Kress R.: Numerical Analysis, Springer - Verlag New York 1998 [7] Přikryl P.: Numerické metody matematické analýzy, SNTL Praha 1988 [8] Ralston A.: Základy numerické matematiky, Academia Praha 1973 [9] Rektorys K.: Aplikovaná matematika I a II, SNTL Praha 1977 [10] Segeth K.: Numerický software I, Karolinum Praha 1998 [11] Vitásek E.: Numerické metody, SNTL Praha 1987

145

Rejstřík algoritmus, 8 aproximace funkcí, 95

chyba řešení, 42 chyby v aritmetických operacích, 10 chyby vzniklé při zobrazení čísla v počítači, 9

Banachův prostor, 21 Bernoulliova metoda, 89 Bernsteinův polynom, 121 Budanova-Fourierova věta, 86 Bézierova křivka, 121 Bézierův polygon, 121 Bézierův polynom, 121

interpolace na ekvidistantních uzlech, 105 iterační metody, 41 Jacobiova iterační metoda, 42 Jacobiova metoda, 44

částečný problém vlastních čísel, 59, 61

korektní úloha, 12

diference vpřed, 105 diference vzad, 105 dobře podmíněná úloha, 12 extrapolace, 106 Gaussova eliminační metoda, 28 Gaussova eliminační metoda s výběrem hlavního prvku, 31 Gaussova-Seidelova iterační metoda, 46 Gaussovy kvadraturní vzorce, 139 Gaussův-Legendr˚ uv kvadraturní vzorec, 141 Graefova metoda, 90 Hilbertův prostor, 19

Lagrangeův interpolační polynom, 97 Laguerrova metoda, 93 Legenderovy polynomy, 21 lichoběžníkové pravidlo, 131, 134 lineární operátor, 22 LU-rozklad matice, 35 metoda bisekce, 75 metoda LU-rozkladu, 33 metoda nejmenších čtverců, 124 metoda prosté iterace, 76 metoda Rayleighova podílu, 63 metoda regula falsi, 80 metoda sečen, 85 mocninná metoda, 61 nejlepší aproximace vzhledem k množině, 26

charakteristický polynom, 15 146

Nevillův algoritmus, 103 řešení soustav nelineárních rovnic, 85 Newtonova metoda, 82 Newtonovy-Cotesovy kvadraturní vzorce, Simpsonovo pravidlo, 132, 134 130 skalární součin, 19 Newtonův interpolační polynom, 101 složené kvadraturní vzorce, 133 norma, 18 splajny, 108 norma matice, 18 spojitý lineární operátor, 22 norma operátoru, 22 stabilita algoritmu, 13 normovaný prostor, 18 stanovení hodnoty determinantu, 39 numerická derivace, 142 stanovení inverzní matice, 40 numerická úloha, 8 superrelaxační metoda, 54, 55 symetrická matice, 14 obdélníkové pravidlo, 133 odhad polohy vlastních čísel, 59 určení vlastních čísel metodou LUortogonální matice, 14 rozkladu, 67 ortogonální prvky, 20 určení vlastních čísel metodou orortonormální prvky, 20 togonálních transformací, 69 ostře diagonálně dominantní matice, úplný normovaný prostor, 21 14 úplný problém vlastních čísel, 59, 65 platné číslice, 11 podmíněnost soustav lineárních rovvlastní číslo matice, 15 nic, 55 podobné matice, 17 Weierstrassova věta, 96 pozitivně definitní matice, 14 pásová matice, 15 přesnost řešení, 42 přímé metody řešení soustav lineárních rovnic, 28 přímý výpočet vlastních čísel, 65 Rayleighův podíl, 17 Richardsonova extrapolace, 107 Rombergova kvadratura, 136 Rungeova metoda dvojnásobného kroku, 135 řešení nelineární rovnice, 73 řešení rovnic Pn (x) = 0, 86 147