Numerické metodyVratislava Mošová

1

Předmluva

S rozvojem počítačů vzrostl význam numerických metod. Řada výpočetníchpostupů, které sem řadíme, vznikla sice už dávno předtím, ale teprve nástupvý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 me-tod jako jsou Gaussova a Jakobiova metoda pro soustavy lineárních rovnic,Newtonova metoda pro nelineární rovnice, Lagrangeova interpolace nebo ob-dé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 polynomynebo 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ěminejzákladnějšími. K pochopení textu je nutná znalost základních pojmů z li-neární algebry a matematické analýzy. Samotný text je rozdělen do šestikapitol. 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 tvr-zení 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 alge-bry, jako je řešení soustav lineárních rovnic, určování vlastních čísel matice avýpočet kořenů nelineárních rovnic. Poslední dvě kapitoly budou zasvěcenynumerickým metodám řešení úloh matematické analýzy, zejména problema-tice aproximace funkcí a otázkám numerické kvadratury.

Na konci skript je seznam doporučené literatury. K připomenutí si po-tř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é proble-matiky. 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 skrip-

3

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 numeric-ké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 dopodstaty 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žívatsoftwarové produkty, které jsou v současné době k dispozici.

4

Obsah

Předmluva 3

1 Úvodní pojmy 7Numerické metody a matematické modelování . . . . . . . . . . . . 7O nepřesnostech při řešení problému . . . . . . . . . . . . . . . . . 9Některá fakta z lineární algebry a funkcionální analýzy . . . . . . . 14

2 Řešení soustav lineárních rovnic 27Přímé metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Gaussova eliminační metoda . . . . . . . . . . . . . . . . . . . 28Metoda LU-rozkladu . . . . . . . . . . . . . . . . . . . . . . . 33

Iterační metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41Jacobiova iterační metoda . . . . . . . . . . . . . . . . . . . . 42Gaussova-Seidelova iterační metoda . . . . . . . . . . . . . . . 46Podmíněnost soustav lineárních rovnic . . . . . . . . . . . . . 55

3 Vlastní čísla a vlastní vektory matice 59Částečný problém vlastních čísel . . . . . . . . . . . . . . . . . . . . 61

Mocninná metoda . . . . . . . . . . . . . . . . . . . . . . . . . 61Metoda Rayleighova podílu . . . . . . . . . . . . . . . . . . . 63

Úplný problém vlastních čísel . . . . . . . . . . . . . . . . . . . . . 65Přímý výpočet vlastních čísel . . . . . . . . . . . . . . . . . . 65Určení vlastních čísel metodou LU-rozkladu . . . . . . . . . . 67Určení vlastních čísel metodou ortogonálních transformací . . 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 . . . . . . . . . . . . . . . . . . . . . . . 76Metoda regula falsi . . . . . . . . . . . . . . . . . . . . . . . . 80Newtonova metoda . . . . . . . . . . . . . . . . . . . . . . . . 82Metoda sečen . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

Řešení rovnic Pn(x) = 0 . . . . . . . . . . . . . . . . . . . . . . . . 86Bernoulliova metoda . . . . . . . . . . . . . . . . . . . . . . . 89Graefova metoda . . . . . . . . . . . . . . . . . . . . . . . . . 90Laguerrova metoda . . . . . . . . . . . . . . . . . . . . . . . . 93

5 Aproximace funkcí 95Interpolace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

Lagrangeův interpolační polynom . . . . . . . . . . . . . . . . 97Newtonův interpolační polynom . . . . . . . . . . . . . . . . . 101Extrapolace . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106Splajny . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108Aproximace trigonometrickými polynomy . . . . . . . . . . . . 115Hermitova interpolace . . . . . . . . . . . . . . . . . . . . . . 119

Bézierovy křivky . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120Metoda nejmenších čtverců . . . . . . . . . . . . . . . . . . . . . . 124

6 Numerická kvadratura a derivace 129Newtonovy-Cotesovy kvadraturní vzorce . . . . . . . . . . . . . . . 130

Složené kvadraturní vzorce . . . . . . . . . . . . . . . . . . . . 133Gaussovy kvadraturní vzorce . . . . . . . . . . . . . . . . . . . . . 139Numerická derivace . . . . . . . . . . . . . . . . . . . . . . . . . . . 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ému. Aby ses nimi vypořádal, vytvořil si řadu nástroju, které mu to umožní. Matematikaje jedním z nich. Tak například úlohu týkající se bezpečné konstrukce mostulze modelovat pomocí diferenciálních rovnic. Jinou úlohu muže představovatproblém, jak optimálně řídit stavbu nějakého areálu, což se dá s určitýmizjednodušeními zachytit prostřednictvím soustavy lineárních rovnic. Řešittakové úlohy přímo však není možné. To je jeden z duvodu, proč je zapo-třebí rozvíjet numerickou matematiku. Dalším důvodem je fakt, že algebraa 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ýchpřípadech dávají konkrétní návod, jak tyto problémy řešit. (Např. v důkazuvěty o existenci a jednoznačnosti řešení Cauchyova počátečního problému sekonstruuje 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 Cramero-vým pravidlem, bude výpočet hodnot determinantů vyšších řádů působitproblé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čivyžaduje, abychom zadali na vstup konečný počet číselných údaju a postup,

7

prostřednictvím kterého po konečném počtu kroku dojdeme k výsledku. Tenje dán opět konečným počtem výstupních číselných údaju. 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 mezikoneč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čitvlastní čí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 ope-rací, 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 prox 6= 0 máme 2x = x+ a

xa x = 1

2(x+ ax). Odtud už obdržíme iterační formuli,

tzv. Heronuv vzorec,

xn =12

(xn−1 +a

xn−1), n = 1, 2, . . .

Algoritmus (Heronuv vzorec)

Vstup: a, x0, npro k = 1, 2, . . . , n

xk = 12(xk−1 + a

xk−1)

Výstup: xn - přibližná hodnota odmocniny√

a

8

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 elimi-nace), 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 sta-bilita 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šíchnepřesností.

Zdrojem chyb mohou být1. 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íhomodelu 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á mantisase 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ž muže zpusobit 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šim-ně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

q+

x2

q2+ . . . ) qb, x1 6= 0,

9

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

x = sgn x (x1

q+

x2

q2+ · · ·+ xt

qt) qb,

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

∆x = |x− x| = qb|∞∑i=1

xi

qi−

t∑i=1

xi

qi| ≤

{qbq−t při řezání,12q

bq−t při zaokrouhlování.

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

|x− x

x| =

qb|∑∞i=1

xi

qi −∑t

i=1xi

qi |qb|∑∞

i=1xi

qi | ≤{

q1−t při řezání,12q

1−t při zaokrouhlování,

neboť |∑∞i=1

xi

qi | ≥ |∑∞i=1

1qi | = 1

q.

Díky zaokrouhlovacím chybám vznikají také chyby ve výsledcích aritmetic-ký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í exponentu 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čemobecně neplatí ani asociativní ani komutativní zákon.

Příklad. V M(10, 3) pomocí částečných součtů odhadněte součet řady∑∞

n=0 0, 9n.

10

Řešení: Protože se jedná o součet geometrické řady s kvocientem q = 0, 9,je přesná hodnota součtu s = 1

1−0,9 = 10. Částečné součty spočteme dvěmazpů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 sk rk s− sk s− rk

50 9.98 9.97 0.02 0.03100 9.98 10.0 0.02 0150 9.98 10.0 0.02 0

V prvním případě, kdy sčítáme od největších členů, se nejdřív rychle blížímek 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. Oprotitomu přičítání odzadu umožňuje do celkového součtu zahrnout příspěvekvšech sčítanců.

Definice 1.3 Řekneme, že v normalizovaném tvaru aproximace

x = sgn x (t∑

i=1

xi

10i) · 10b

čísla

x = sgn x (∞∑i=1

xi

10i) · 10b

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íchdatech. 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žinavstupních dat a Y množina výstupních dat. Řekneme, že úloha U je korektní,když

∀x ∈ X ∃! y ∈ Y : y = Ux,

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

cp =relativní chyba na výstupurelativní 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 = 1003 .

1Banachův prostor viz definice 1.16.

12

Řešení: Protože

p(33) = −28, p(1003

) = −509

, cp =|−28− 50

9−28 |

|33− 509

33 |=

22,4281/333

.= 79,

je úloha špatně podmíněná.

Úspěch numerického řešení problému také závisí na volbě algoritmu. Pokudmalá nepřesnost ve vstupních datech začne generovat velké chyby běhemvýpočtu, nepovede zvolený postup ke správnému výsledku. Zda určitý algo-ritmus je vhodný pro daný problém či ne, souvisí s otázkou stability algo-ritmu. Algoritmus je stabilní, když je málo citlivý na poruchy v datech a navliv zaokrouhlovacích chyb během výpočtu. Stabilitu algoritmů lze studovatpomocí numerických experimentů.

Příklad Je dána rovnice x2−2bx+c = 0. Její kořeny lze určit prostřednictvímalgoritmu

A1: x1 = b +√

b2 − c, x2 = b−√b2 − cnebo algoritmu

A2: x1 = b +√

b2 − c, x2 = (b−√b2−c)(b+√

b2−c)b+√

b2−c= c

b+√

b2−c. (viz [3])

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 = x1x2.

Řešení: Pro rovnici x2 − 105x + 1 = 0 je v M(10, 8) : x1 = 100000, 00.Když spočteme druhý kořen pomocí algoritmu A1, dostaneme x2 = 0 ax1x2 = 0 6= c.Když využijeme algoritmus A2, je x2 = 0, 0000100000000 a x1x2 = 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ý mo-del. To znamená zadat jednoznačný a srozumitelný funkční vztah mezi zada-nými a hledanými matematickými objekty (čísly, maticemi, funkcemi, křiv-kami).

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ýsledku.

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 kurzechalgebry 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 budouzejmé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∑

k=1,k 6=i

|aik|. (1.1)

Symetrická, když pro všechna i, j = 1, . . . , n platí

aij = aji.

Ortogonální, když pro všechna i, j = 1, . . . , n platí 2

n∑

k=1

aikakj = δij. (1.2)

Pozitivně definitní, když

n∑i=1

n∑j=1

aijxixj > 0 pro libovolný vektor ~x = (x1, . . . , xn). (1.3)

14

Čtvercová matice je (p, q)-pásová, když existují čísla p, q < n tak, že proprvky 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

A =

2 −1 0−1 5 2

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 řádun. Ř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řenu λ1, λ2, . . . , λn, které nemusejí být navzájem různé. Pak se ale můžestát, že počet lineárně nezávislých vlastních vektorů je menší než je řád ma-tice (viz [3] strana 107).

2

δij =

{0 pro i 6= j1 jinak

15

Příklad Stanovte vlastní čísla a vlastní vektory matice A =

1 1 00 2 00 0 3

.

Řešení:

|A− λI| =

∣∣∣∣∣∣

1− λ 1 00 2− λ 00 0 3− λ

∣∣∣∣∣∣= (2− λ)(λ− 3)(λ− 1) = 0.

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 = 0tzn.

0 1 0 | 00 1 0 | 00 0 2 | 0

−→

(0 1 0 | 00 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 | 00 0 1 | 0

−→

( −1 1 0 | 00 0 1 | 0

),

~v2 = (r, r, 0)T , r ∈ R.

Pro λ3 = 3 je−2 1 0 | 00 −1 0 | 00 0 0 | 0

−→

( −2 1 0 | 00 −1 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í

λ =~vT A~v

~vT~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 po-jmem podobné matice.

Definice 1.9 Řekneme, že matice A je podobná matici B, když existujeregulární matice T tak, že

T−1AT = B.

Píšeme pak A ∼ B.

Při konstrukci numerických metod pro výpočet vlastních čísel matice se vy-uží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

A =

1 0 00 2 00 0 2

∣∣∣∣∣∣a B =

1 1 | 00 2 | 00 0 2

sice mají stejná vlastní čísla, ale nejsou podobné, protože jsou tvořeny na-vzájem různými Jordanovými bloky

Následující systém definic a vět představuje obecný aparát, který využijemepř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žkyMn — množina čtvercových matic řádu n,Pn — množina polynomů n-tého stupně,Cn(a, b) — množina funkcí n-krát spojitě diferencovatelných na (a, b),Lp(a, b) — množina funkcí, pro které je (Lebesgueův) integrál

∫ b

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ýmlineárním prostorem.

Příklad Pro ~x ∈ Rn představuje

‖~x‖∞ = maxi=1,...,n

{|xi|} nebo ‖~x‖E =

√√√√n∑

i=1

x2i

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

||A||∞ = maxi=1,...,n

n∑j=1

|aij|. (1.8)

Sloupcová norma

||A||S = maxj=1,...,n

n∑i=1

|aij|. (1.9)

18

Eukleidovská norma

||A||E =

√√√√n∑

i=1

(n∑

j=1

a2ij). (1.10)

Pro f ∈ C0(a, b) lze vzít

||f(x)||∞ = maxx∈<a,b>

|f(x)|. (1.11)

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

||f(x)||p = (∫ b

a

|f(x)|p dx)1p . (1.12)

Poznámka Říkáme, že dvě normy ‖.‖1 a ‖.‖2 v lineárním prostoru X jsouekvivalentní, právě když existují kladné konstanty c1, c2 ∈ R tak, že provšechna x ∈ X platí

c1‖x‖1 ≤ ‖x‖2 ≤ c2‖x‖1.

Protože v konečnědimenzionálních prostorech jsou všechny normy ekviva-lentní, 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é provš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é definovatskalární součin3

(f, g) =∫ b

a

f(x)g(x) dx. (1.15)

(Diskrétní případ) Když funkce f a g jsou definovány v uzlových bodechx0, x1, . . . , xn ∈ 〈a, b〉, je skalární součin

(f, g) =n∑

j=0

f(xj)g(xj) dx. (1.16)

Poznámka V Hilbertově prostoru X lze definovat normu vztahem

||x|| =√

(x, x), x ∈ X.

Příklad Nechť i =√−1 je imaginární jednotka. V L2(−π, π) stanovte

normu funkce f(x) = eix.

Řešení: ‖eix‖2 =∫ π

−πeixe−ix dx = 2π, ‖eix‖ =

√2π.

Definice 1.13 Nechť X je Hilbertův prostor. Řekneme, že prvky xi, xj ∈ Xjsou ortogonální, když

(xi, xj)

{= 0 pro i 6= j,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í sys-té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 〈−π, π〉.

3Funkce 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〈−π, π〉.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

Pn(x) =1

2n n!dn(x2 − 1)n

dxn,

jsou ortogonální na intervalu 〈−1, 1〉.Ukažme si nyní, že funkce ϕj = eijx (j = 0,±1, . . . ) jsou opravdu na intervalu〈−π, π〉 ortogonální.Spojitý případ:

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

−π

eijxe−ikx dx =

{0 pro j 6= k,

2π pro j = k.

Diskrétní případ: Když xs = 2πsn+1 , s = 0, 1, . . . n, jsou uzlové body z intervalu

〈−π, π〉, pak

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

s=0

ei(j−k) 2πsn+1

je n-tým částečným součtem geometrické řady, která má kvocient q = ei(j−k) 2πn+1 .

Pokud j−kn+1 je celé, je q = 1 a (ϕj(x), ϕk(x)) = n + 1. V opačném případě je

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

Definice 1.14 Nechť X je normovaný prostor. Řekneme, že posloupnostprvků xn ∈ X konverguje k x ∈ X, když

||xn − x|| → 0, xn →∞.

Definice 1.15 Nechť X je normovaný prostor. Posloupnost {xn} je cauchy-ovská v X, když

||xm − xn|| → 0, n, m →∞.

Definice 1.16 Normovaný prostor, ve kterém je každá cauchyovská posloup-nost konvergentní, je úplný normovaný neboli Banachův prostor.

21

Příklad Prostor C(0)(a, b) s normou ||.||∞ je úplný. Konvergence posloup-nosti spojitých funkcí v normě ||.||∞ je totiž ekvivalentní se stejnoměrnoukonvergencí těchto funkcí. Limita posloupnosti spojitých funkcí je tedy opětspojitá funkce.Prostor Lp(a, b) s normou ||.||p je úplný.

Definice 1.17 Nechť X, Y jsou lineární prostory. Zobrazení A : X → Y jelineá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 ope-rátoru A. Platí

||A|| = sup||x||=1

||Ax||.

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

Věta 1.2 V normovaném prostoru je lineární operátor spojitý, když jeomezený.

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átorA : U → X nazveme kontrakcí, když existuje kontanta q ∈ 〈0, 1) tak, že

||Ax− Ay|| ≤ q||x− y|| (1.17)

22

pro všechna x, y ∈ U.

Definice 1.21 Nechť X je Banachův prostor. Řekneme, že x ∈ X je pevnýmbodem operátoru A : X → X, když

Ax = x.

Věta 1.3 (Banach) Nechť U je úplný podprostor normovaného prostoru Xa 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 jepevně 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|| ≤≤ q2 ||xn−1 − xn−2|| ≤ · · · ≤ qn ||x1 − x0||.

Pro m > n dostaneme

||xn − xm|| ≤ ||xn − xn+1||+ ||xn+1 − xn+2||+ · · ·+ ||xm−1 − xm|| ≤≤ (qn + qn+1 + · · ·+ qm−1)||x1 − x0||. (1.18)

Protože A je kontrakce, máme pro n →∞

||xn − xm|| ≤ qn

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

To znamená, že posloupnost {xn} je cauchyovská v úplném prostoru U, atedy i konvergentní. Ze spojitosti operátoru A máme

x = limn→∞

xn+1 = limn→∞

Axn = Ax.

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éhopevného bodu operátoru A, ale její důkaz poskytuje i návod, jak pevný bodnají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

||xn − x|| ≤ qn

1− q||x1 − x0|| (1.19)

a aposteriorní odhad

||xn − x|| ≤ q

1− q||xn − xn−1||. (1.20)

Důkaz. Existence pevného bodu plyne z Banachovy věty 1.3. Apriorní od-had obdržíme ze vztahu (1.18), když m → ∞. Aposteriorní odhad plynez 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

1− q||x1 − x0|| < ε,

qn <1− q

||x1 − x0|| ε,

n ≥ ln ε

ln q, kde ε =

(1− 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ž od-had 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||n1− ||H|| ||x1 − x0||

a aposteriorní odhad

||xn − x|| ≤ ||H||1− ||H|| ||xn − xn−1||.

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 propostupné aproximace máme

||xn|| ≤n∑

i=0

||H iz|| ≤n∑

i=0

||H i|| ||z|| ≤ 11− ||H|| ||z||.

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

||(I −H)−1z|| ≤ ||z||1− ||H|| .

25

Definice 1.22 Nechť X je normovaný prostor, U ⊂ X a x ∈ X. Prvek v ∈ Unazveme nejlepší aproximací x vzhledem k U, když

||x− v|| = infu∈U

||x− u||. (1.21)

Poznámka Prvek v ∈ U je nejlepší aproximací x ∈ X, právě když provšechna u ∈ U je

(x− v, u) = 0. (1.22)

Poznámka Ke každému prvku normovaného prostoru existuje nejlepší apro-ximace 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 prostoruX a u1, . . . , un jeho báze. Prvek

v =n∑

i=1

αiui (1.23)

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

vyhovují soustavěn∑

i=1

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

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 jehobáze. Nejlepší aproximaci v prvku u pak lze psát

v =n∑

i=1

(x, ui)ui. (1.25)

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

26

Kapitola 2

Řešení soustav lineárníchrovnic

Řešení soustav lineárních rovnic je velice frekventovaný problém. Na úlohunajít řešení soustavy lineárních rovnic vedou nejruznější optimalizační úlohynebo 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 sepouze 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 metodyposkytují přesné řešení v konečném počtu kroků, pokud během výpočtu ne-zaokrouhlujeme.

• Metody nepřímé tzv. iterační (Gaussova-Seidelova metoda, Jacobiova me-toda). Těmito metodami získáme pouze aproximativní řešení.

Obecně nelze říci, který z uvedených typů pro který okruh řešených úlohje nejvýhodnější použít. Pro první orientaci lze zhruba uvést, že soustavys 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

a011x1 + a0

12x2 + · · ·+ a01nxn = a0

1,n+1, (2.2)

a021x1 + a0

22x2 + · · ·+ a02nxn = a0

2,n+1,

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

a0n1x1 + a0

n2x2 + · · ·+ a0nnxn = a0

n,n+1.

Nyní budeme postupně eliminovat prvky pod hlavní diagonálou matice sou-stavy.Nejprve 1. rovnici využijeme k eliminaci 1. neznámé ze zbývajících n − 1rovnic: 1. řádek vynásobíme čıslem m0

21 = −a021

a011

a přičteme k 2. řádku,

1. řádek vynásobený číslem m031 = −a0

31a0

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

získáme soustavu

a011x1+ a0

12x2+ . . . + a01nxn = a0

1,n+1,a1

22x2+ . . . + a12nxn = a1

2,n+1,. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .a1

n2x2+ . . . + a1nnxn = a1

n,n+1.

Obecně v k-tém kroku (k = 1, 2, . . . , n−1) vynulujeme prvky v k-tém sloupci

tím, že k-tý řádek vynásobený multiplikátorem mk−1ik = −ak−1

ik

ak−1kk

přičteme k i-

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

a011x1+ a0

12x2+ . . . +a01nxn = a0

1,n+1,a1

22x2+ . . . +a12nxn = a1

2,n+1,a2

33x3+ . . . +a23nxn = a2

3,n+1,. . . . . . . . . . . . . . . . . .an−1

nn xn = an−1n,n+1,

28

z níž zpětnou eliminací získáme hodnoty neznámých xn, xn−1, . . . , x1.

Algoritmus (GEM)

Vstup: n, A = (a0ij), ~b = (a0

i,n+1)pro k = 1, 2, . . . , n− 1

pro i = k + 1, k + 2, . . . , n

mk−1ik = −ak−1

ik

ak−1kk

(pokud ak−1kk 6= 0)

pro j = k + 1, k + 2, . . . , n + 1ak

ij = ak−1ij + mk−1

ik ak−1kj

pro i = n, n− 1, . . . , 1xi = 1

ai−1ii

(ai−1i,n+1 −

∑nj=i+1 ai−1

ij xj)

Výstup: ~x = (x1, x2, . . . , xn)

Poznámka Protože platín∑

k=1

k =n(n + 1)

2,

n∑

k=1

k2 =n(n + 1)(2n + 1)

6,

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∑

k=1

(n− k) = (n− 1)n− (n− 1)n2

=n(n− 1)

2

dělení a k výpočtu prvků akij potřebujeme

n−1∑

k=1

(n− k)(n− k + 1) =n−1∑

k=1

n(n + 1)− k(2n + 1) + k2 =n3 − n

3

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

n2 + n

2

29

násobení an2 − n

2sčítání. Celkem v našem případě tedy potřebujeme n3

3 + n2 − n3 násobení a

dělení a n3

3 + n2

2 − 5n6 sčítání. Z hlediska doby potřebné k realizaci algoritmu

jsou operace násobení a dělení v podstatě totožné. Protože časová náročnostná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 funkceg, když existuje konečná limita

limx→a

f(x)g(x)

.

Píšeme pak f(x) = O(g(x)), x → a.

Poznámka K realizaci algoritmu GEM je potřeba n3

3 + O(n2) násobení.To znamená, že když v soustavě počet rovnic zdvojnásobíme, vzroste dobavýpočtu osmkrát.

Příklad Gaussovou eliminační metodou řešte soustavu

x1 + 2x2 − 2x3 = 1,x1 + x2 + x3 = 3,2x1 + 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 m0

21 = −1, m031 =

−2; m132 = −2. Dostaneme

1 2 −2 | 11 1 1 | 32 2 1 | 5

−→

1 2 −2 | 10 −1 3 | 20 −2 5 | 3

−→

1 2 −2 | 10 −1 3 | 20 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 krokuje hlavní prvek ak−1

kk = 0. Pokud taková situace nastane a některé ak−1kk se v po-

čítači zobrazí jako nula, můžeme zkusit použít Gaussovu eliminační metodus 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 beremenejvětší prvek v k-tém sloupci. Vybíráme mezi řádky, z nichž jsme dosudnevzali vedoucí prvek. Hlavní prvek je v řádku p a platí pro něj

|apk| = maxi

|ak−1ik |.

ii) GEM s řádkovým výběrem — za hlavní prvek v k-tém kroku beremenejvětší prvek v k-tém řádku. Vybíráme mezi sloupci, z nichž jsme dosudnevzali vedoucí prvek. Hlavní prvek je ve sloupci q a platí pro něj

|akq| = maxj

|ak−1kj |.

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| = maxij

|ak−1ij |.

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 | 11 2 1 | 42 2 1 | 5

−→

1 2 −2 | 10 0 3 | 30 −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. sloupcipuvodní matice má největší absolutní hodnotu prvek a0

31. Proto přesuneme

31

3. řádek na začátek a dále už pokračujeme standardní Gaussovou eliminačnímetodou:

1 2 −2 | 11 2 1 | 42 2 1 | 5

→

2 2 1 | 51 2 −2 | 11 2 1 | 4

,

2 2 1 | 51 2 −2 | 11 2 1 | 4

−→

2 2 1 | 50 1 −5

2 | −32

0 1 12 | 3

2

−→

2 2 1 | 50 1 −5

2 | −32

0 0 62 | 6

2

,

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ěremhlavního prvku může přispět ke snížení kumulování zaokrouhlovacích chybv 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 | 1001 1 | 1

)−→

(1 200 | 1000 −199 | −99

),

x2 =99199

= 0, 4974 . . . ,

x1 = 100− 20099199

=100199

= 0, 5025 . . .

Když během výpočtu zaokrouhlujeme, pak v M(10, 2) bude(

1 200 | 1001 1 | 1

)−→

(1 200 | 1000 −200 | −99

),

x2 =99200

.= 0, 50,

x1 = 100− 200 · 0, 50 = 0.

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 eli-minač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ímev původní matici první a druhý sloupec a pokračujeme eliminační metodou(

1 200 | 1001 1 | 1

)→

(200 1 | 1001 1 | 1

)−→

(200 1 | 1000 0, 995 | 0, 455

).

Při zaokrouhlení v M(10, 2)

x1 =0, 4550, 995

.= 0, 46,

x2 =100− 0, 46

200.= 0, 50.

Při použití Gaussovy eliminační metody s řádkovým výběrem jsme tak do-cílili lepších výsledku 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í maticesoustavy prostřednictvím lineární transformace na horní trojúhelníkovou ma-tici, ze které pak zpětnou eliminací stanovíme řešení. Když A je puvodní ma-tice soustavy, U je příslušná trojúhelníková matice a T = Tn−1Tn−2 . . . T2T1,kde pro k = 1, 2, . . . , n− 1 je

Tk =

1 0 . . . 0 0 0 . . . 00 1 . . . 0 0 0 . . . 0

. . . . . . . . . . . . . . . . . . . . . . . .0 0 mk−1

k+1,k 1 0 0 . . . 00 0 mk−1

k+2,k 0 1 0 . . . 0. . . . . . . . . . . . . . . . . . . . . . . .0 0 mk−1

n,k 0 0 0 . . . 1

,

33

pak proces převedení matice A na trojúhelníkovou matici lze zapsat jako

TA = U.

Matice T je regulární, a proto k ní existuje matice inverzní L = T−1. Platítedy

A = T−1U = LU.

Tím jsme vyjádřili matici A jako součin dolní trojúhelníkové matice

L =

l11 0 0 . . . 0l21 l22 0 . . . 0l31 l32 l33 . . . 0. . . . . . . . . . . . . . .ln1 ln2 ln3 . . . lnn

a horní trojúhelníkové matice

U =

u11 u12 u22 . . . u1n

0 u22 u23 . . . u2n

0 0 u33 . . . u3n

. . . . . . . . . . . . . . .0 0 0 . . . unn

.

Proto zde také hovoříme o tzv. LU-rozkladu.

Uvedený rozklad je ekvivalentní s Gaussovou eliminační metodou, a proto jemožné ho realizovat, právě když lze realizovat GEM.

Příklad LU-rozklad nelze aplikovat přímo, když A =

(0 −53 0

). V této

matici pivot a011 = 0, a proto nelze stanovit multiplikátor m0

21.

Příklad Pro matici A =

(2 54 1

)LU-rozklad existuje. Lehce ověříme, že

A =

(1 02 1

)(2 50 −9

)nebo také A =

(2 04 −1

)(1 5

20 9

).

34

Poznámka Z posledního příkladu je vidět, že rozklad A = LU není určenjednoznač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. Jedno-značnost rozkladu si zajistíme třeba tak, že zadáme prvky na hlavní diagonálematice 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

a21 a22 a23

a31 a32 a33

=

1 0 0l21 1 0l31 l32 1

u11 u12 u13

0 u22 u23

0 0 u33

.

u11 = a11, l21u11 = a21 =⇒ l21 = a21u11

,l31u11 = a31 =⇒ l31 = a31

u11,

u12 = a12, l21u12 + u22 = a22 =⇒ u22 = a22 − l21u12,l31u12 + l32u22 = a32 =⇒ l32 = a32−l31u12

u22,

u13 = a13, l21u13 + u23 = a23 =⇒ u23 = a32 − l21u13,l31u13 + l32u23 + u33 = a33 =⇒ u33 = a33 − l31u13 − l32u23.

Algoritmus (LU-rozklad)

Vstup: n, A = (aij)pro j = 1, 2, . . . , n

pro i = 1, 2, . . . , juij = aij −

∑i−1r=1 lirurj

pro i = j + 1, j + 2, . . . , nlij = 1

ujj(aij −

∑j−1s=1 lisusj) (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ímsoustavy (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

A =

1 2 −21 1 12 2 1

, L =

1 0 01 1 02 2 1

, U =

1 2 −20 −1 30 0 −1

.

Eliminace ze soustavy L~y = ~b :

1 0 0 | 11 1 0 | 32 2 1 | 5

,

y1 = 1, y2 = 2, y3 = −1.

Zpětná eliminace z U~x = ~y :

1 2 −2 | 10 −1 3 | 20 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šitví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 dife-renciálních rovnic metodou konečných prvku, se setkáváme s maticemi, kterémají speciální tvar. Na závěr si proto uveďme dva příklady algoritmu LU-rozkladu pro takové matice.

i) Pokud je matice A symetrická pozitivně definitní (viz definice 1.6), dá seukázat, že existuje horní trojúhelníková matice U s kladnými prvky na hlavnídiagonále tak, že A = UT 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, 2, . . . , n

ujj =√

ajj −∑j−1

r=1 u2jr

pro i = j + 1, j + 2, . . . , nuij = 1

ujj(aij −

∑j−1s=1 uisujs)

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.

Řešení: Označme A =

1 1 11 5 51 5 14

, ~b = (2, 6, 6)T .

37

Platí A = UT U, když UT =

1 0 01 2 01 2 3

, U =

1 1 10 2 20 0 3

.

Řešení soustavy UT~y = ~b :

1 0 0 | 21 2 0 | 61 2 3 | 6

, ~y = (2, 2, 0)T .

Ze soustavy U~x = ~y máme

1 1 1 | 20 2 2 | 20 0 3 | 0

, ~x = (1, 1, 0)T .

ii) Dalším typem speciálních matic, které se v praktických problémech častovyskytují, jsou pásové matice. (Pásová matice viz definice 1.6.) V takovémpřípadě stačí pracovat pouze s nenulovými prvky matice. Sníží se pak početoperací 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, j − p)uij = aij −

∑i−1r=α lirujr

pro j = i + 1, i + 2, . . . , γα = max(1, j − p)lji = 1

uii(aji −

∑i−1r=α ljruri)

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 + x2 = 0,2x1 + x2 + 2x3 = 5,

3x2 − x3 + 3x4 = 5,4x3 + x4 = 5.

Řešení: Označme A =

−1 1 0 02 1 2 00 3 −1 30 0 4 1

, ~b = (0, 5, 5, 5)T .

Platí, že A = LU, když L =

1 0 0 0−2 1 0 00 1 1 00 0 −4/3 1

a U =

−1 1 0 00 3 2 00 0 −3 30 0 0 5

.

Řešení soustavy L~y = ~b :

1 0 0 0 | 0−2 1 0 0 | 5

0 1 1 0 | 50 0 −4/3 1 | 5

, ~y = (0, 5, 0, 5)T .

Ze soustavy U~x = ~y obdržíme

−1 1 0 0 | 00 3 2 0 | 50 0 −3 3 | 00 0 0 5 | 5

, ~x = (1, 1, 1, 1)T .

Gaussovu eliminační metodu a LU-rozklad lze s výhodou použít i k jinýmúčelum než jen k řešení soustav lineárních rovnic. Například mužeme vypo-čítat hodnotu determinantu nebo určit tvar inverzní matice.

i) Stanovení hodnoty determinantu: Pomocí GEM převedeme matici A nahorní trojúhelníkovou matici. Pak

det A = a011a

122 . . . an−1

nn .

39

Pokud jsme provedli LU-rozklad matice A, je

det A = detL detU = 1 · u11u22 . . . unn.

Počet násobení při Gaussově eliminaci je řádově n3

3 . To je při výpočtu deter-minantů vyšších řádů o hodně méně než n! násobení, která bychom muselirealizovat při výpočtu determinantu rozvojem.

Příklad Spočítejte determinant matice A =

1 −1 12 −1 23 −1 1

pomocí Gaus-

sovy metody a metody LU-rozkladu.

Řešení:

Gaussova metoda: A =

1 −1 12 −1 23 −1 1

−→

1 −1 10 1 00 0 −2

, det A = −2.

LU-rozklad: A = LU =

1 0 02 1 03 2 1

1 −1 10 1 00 0 −2

, det A = −2.

ii) Obě metody je možné také využít při stanovení inverzní matice: ModifikacíGEM je Gaussova-Jordanova metoda. Matici tvořenou puvodní maticí Aa 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žijemeskutečnosti, že když A = LU, pak

A−1 = U−1L−1.

K realizaci Gaussovy-Jordanovy metody potřebujeme řádově n3

2 +O(n2) ope-rací, 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

stanovte inverzní matici k matici A =

1 −1 12 −1 23 −1 1

.

40

Řešení:Gaussova-Jordanova metoda:

(A|I) =

1 −1 1 | 1 0 02 −1 2 | 0 1 03 −1 1 | 0 0 1

−→

1 0 0 | −1/2 0 1/20 1 0 | −2 1 00 0 1 | −1/2 1 −1/2

= (I|A−1).

LU-rozklad:

A =

1 −1 12 −1 23 −1 1

=

1 0 02 1 03 2 1

1 −1 10 1 00 0 −2

,

A−1 = U−1L−1 =

1 0 0−2 1 01 −2 1

1 −1 1/20 1 00 0 −1/2

=

−1/2 0 1/2−2 1 0−1/2 0 1/2

.

Iterační metody

Předností iteračních metod je, že nekladou takové nároky na paměť počítačejako metody přímé a že jejich použití je univerzální. Můžeme pomocí nichnejen získat řešení soustavy lineárních rovnic, ale také je můžeme využít kezpř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é rovnicesoustavy A~x = ~b vyjádříme právě jednu neznámou v závislosti na zbývajícíchnezná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 pruběhu výpočtu mění, ho-voříme o nestacionárním iteračním procesu, v opačném případě o iteračnímprocesu stacionárním. V dalším textu se omezíme na studium stacionárníchiteračních procesu.

a) Jacobiova iterační metoda

Při odvozování iterační metody budeme postupovat tak, že z i-té rovnicesoustavy

a11x1 + a12x2 + · · ·+ a1nxn = b1,

a21x1 + a22x2 + · · ·+ a2nxn = b2, (2.5)

. . .

an1x1 + an2x2 + · · ·+ annxn = bn

42

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

x(k+1)1 =

1a11

(b1 − a12x(k)2 − a13x

(k)3 − · · · − a1nx(k)

n ),

x(k+1)2 =

1a22

(b2 − a21x(k)1 − a23x

(k)3 − · · · − a2nx(k)

n ),

. . .

x(k+1)n =

1ann

(bn − an1x(k)1 − an2x

(k)2 − · · · − an,n−1x

(k)n−1).

Pokud označíme

HJ =

0 −a12a11

. . . −a1n

a11−a21a22

0 . . . −a2n

a22...

.... . .

...− an1

ann− an2

ann. . . 0

a gJ =

b1a11...

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),

M =

0 0 . . . 0 0a21 0 . . . 0 0a31 a32 . . . 0 0. . . . . . . . . . . . . . .an1 an2 . . . an,n−1 0

, N =

0 a12 . . . a1,n−1 a1n

0 0 . . . a2,n−1 a2n

. . . . . . . . . . . . . . .0 0 . . . 0 an−1,n

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. (2.7)

43

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, . . . . (2.8)

Algoritmus (Jacobiova metoda)

Vstup: n, A = (aij)ni,j=1,

~b = (bi)ni=1, ~x

(0) = (x(0)i )n

i=1, m

pro k = 0, 1, . . . ,m

pro i = 1, 2, . . . , n

x(k+1)i = 1

aii(bi −

∑i−1j=1 aijx

(k)j −∑n

j=i+1 aijx(k)j )

Výstup: ~x(m) = (x(m)1 , x

(m)2 , . . . , x

(m)n )

Příklad Jacobiovou metodou řešte soustavu

2x1 + x2 = 2,

x1 + 4x2 + x3 = 0,

x2 + 2x3 = −2.

Řešení: Iterační formule

x(k+1)1 =

12

(2− x(k)2 ),

x(k+1)2 =

14

(−x(k)1 − x

(k)3 ),

x(k+1)3 =

12

(−2− x(k)2 ).

Iterační matice

HJ =

0 −12 0

−14 0 −1

40 −1

2 0

,

~gJ = (1, 0,−1)T .

44

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

x(k+1)1 = 4 − 3x

(k)2 ,

x(k+1)2 = −2 + 3x

(k)1 .

Iterační matice

HJ =

(0 −33 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 kon-vergence 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á, žev (k + 1)-ním kroku (k = 0, 1, 2, . . . ) obdržíme iterační formule

x(k+1)1 =

1a11

(b1 − a12x(k)2 − a13x

(k)3 − · · · − a1nx(k)

n ),

x(k+1)2 =

1a22

(b2 − a21x(k+1)1 − a23x

(k)3 − · · · − a2nx

(k)n ),

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

x(k+1)n =

1ann

(bn − an1x(k+1)1 − an2x

(k+1)2 − · · · − an,n−1x

(k+1)n−1 ).

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

(M + D)~x = −N~x +~b.

~x = −(M + D)−1N~x + (M + D)−1~b. (2.9)

Dostáváme novou iterační metodu

~x(k+1) = HGS~x(k) + ~gGS,

kdeHGS = −(M + D)−1N, 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 rovnicA~x = ~b je určena iteračními rovnicemi

~x(k+1) = −(M + D)−1N~x(k) + (M + D)−1~b, k = 0, 1, 2, . . . . (2.11)

46

Algoritmus (Gaussova-Seidelova metoda)

Vstup: n, A = (aij)ni,j=1,

~b = (bi)ni=1, ~x

(0) = (x(0)i )n

i=1, m

pro k = 0, 1, . . . , m

pro i = 1, 2, . . . , n

x(k+1)i = 1

aii(bi −

∑i−1j=1 aijx

(k+1)j −∑n

j=i+1 aijx(k)j )

Výstup: ~x(m) = (x(m)1 , x

(m)2 , . . . , x

(m)n )

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

x(k+1)1 =

12

(2− x(k)2 ),

x(k+1)2 =

14

(−x(k+1)1 − x

(k)3 ) =

14

(−1 +x

(k)2

2− x

(k)3 ),

x(k+1)3 =

12

(−2− x(k+1)2 ) =

12

(−74− 1

8x

(k)2 +

14x

(k)3 ).

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

HGS =

0 −12 0

0 −18 −1

40 − 1

1618

, ~gGS = (1,−1

4,−7

8)T .

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 šestiiterací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 jeto 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) jsouvlastní čí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) = maxi{λi(H)} < 1.

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 = TJT−1 a Hk = TJkT−1.Když O je nulová matice a podle předpokladů věty

1 > ρ(H) = λ1(H) > λ2(H) > · · · > λn(H),

pak limk→∞

Jk = O, a tedy limk→∞

Hk = O.

Obráceně, když limk→∞

Hk = O, je také limk→∞

Jk = O, a tedy λi(H) < 1 pro

všechna i.

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ňujenerovnost

||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í

‖~x−~x(k+1)‖ ≤ ‖H‖ ‖~x−~x(k)‖ ≤ ‖H2‖ ‖~x−~x(k−1)‖ ≤ · · · ≤ ‖Hk+1‖ ‖~x−~x(0)‖.

Protože ‖Hk+1‖ ≤ ‖H‖k+1 a ||H|| ≤ q < 1, je limk→∞ ‖Hk+1‖ = 0. Toznamená, že limk→∞ ‖~x− ~x(k+1)‖ = 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, serozhodneme podle toho, která metoda bude konvergovat. K tomu využijemevětu 2.2. Matice soustavy

A =

1 2 −21 1 12 2 1

.

Iterační matice Jacobiovy metody

HJ =

0 −2 2−1 0 −1−2 −2 0

.

Protože

|HJ − λI| =

−λ −2 2−1 −λ −1−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

HGS = −(M+D)−1N =

1 0 0−1 1 0

0 −2 1

0 2 −20 0 10 0 0

=

0 −2 20 2 −30 0 2

.

|HGS−λI| =

−λ −2 2

0 2− λ −30 0 2− λ

= (2−λ)2(−λ) = 0 ⇒ λ1 = 0, λ2,3 = 2,

tzn. spektrální poloměr ρ(HGS) = 2 > 1. Gaussova-Seidelova metoda nemůžepodle věty 2.2 konvergovat.

Úlohu budeme řešit Jacobiovou metodou. Iterační formule mají tvar

x(k+1)1 = 1− 2x

(k)2 + 2x

(k)3 ,

x(k+1)2 = 3− x

(k)1 − x

(k)3 ,

x(k+1)3 = 5− 2x

(k)1 − 2x

(k)2 .

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∞ = maxi

n∑j=1,j 6=i

|aij

aii

| < 1,

50

nebo

qS = maxj

n∑i=1,i6=j

|aij

aii

| < 1,

nebo

qE =

n∑i,j=1,i 6=j

|aij

aii

|2

1/2

< 1.

Pak pro apriorní odhad chyby dostáváme

||~x(n) − ~x(0)||i ≤ qni

1− qi

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

Pro aposteriorní odhad chyby je

||~x(n) − ~x(0)||i ≤ qi

1− qi

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

Důkaz. Iterační matice HJ = −D−1(M +N) (viz str. 43) má na hlavní diago-nále nulové prvky, jinak jsou rovny −ai,j

aii. Čísla q∞, qS, qE představují normy

Jacobiho iterační matice. Pokud jsou menší než 1, jsou splněny předpokladyvěty 1.5. Odtud tvrzení.

Poznámka Z podmínky qS < 1 plyne

n∑i=1,i6=j

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

To znamená, že A je ostře diagonálně dominantní. Jacobiova metoda pakkonverguje 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 Gaussova-Seidelova metoda konverguje pro libovolné ~x(0).

51

Věta 2.6 (Sassenfeldovo kriterium) Gaussova-Seidelova metoda konvergujepro libovolné ~g a libovolné ~x(0), když

p = maxi=1,...,n

pi < 1,

kde

p1 =n∑

j=2

|a1j

a11| < 1, pi =

i−1∑j=1

|aij

aii

|pj +n∑

j=i+1

|aij

aii

| < 1, i = 2, . . . , n.

Pro apriorní odhad chyby pak platí

||~x(k) − ~x(0)||∞ ≤ qn

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

a pro aposteriorní odhad chyby je

||~x(k) − ~x(0)||∞ ≤ q

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

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

xi = −i−1∑j=1

aij

aii

xj −n∑

j=i+1

aij

aii

zj, i = 1, . . . , n.

Indukcí obdržíme |xi| ≤ pi, i = 1, . . . , n, a tedy ||~x||∞ ≤ p. Odtud

||(M + D)−1N ||∞ ≤ 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

A =

2 −1 0 0−1 2 −1 0

0 −1 2 −10 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

p1 =12

< 1, p2 =12· 1

2+

12

< 1, p3 =12· 3

4+

12, p4 =

12· 7

8+

12

< 1,

jsou splněny předpoklady věty 2.6, a tedy Gaussova-Seidelova metoda budekonvergovat.

V případě, že spektrální poloměr ρ(H) je číslo blízké 1, mohou iterace kon-vergovat 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 nu-merickou 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

x(k+1)i = x

(k)i +

ω

aii

(bi −i−1∑j=1

aijx(k)j −

n∑j=i

aijx(k)j ), i = 1, . . . , n.

Věta 2.7 Nechť HJ je iterační matice Jacobiovy metody, λ1 > · · · > λn jsouvlastní čí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 číslu1 − ω + ωλ matice (1 − ω)I + ωHJ = HJω. Vlastní čísla matice HJω jsoureá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í číslobudou 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-Seidelovuiterační formuli (2.11) na tvar

x(k+1)i = x

(k)i + r

(k+1)i , kde r

(k+1)i =

1aii

(bi −i−1∑j=1

aijx(k+1)j −

n∑j=i

aijx(k)j )

a i-té reziduum r(k+1)i vynásobíme vhodným parametrem ω. Obdržíme novou

iterační formuli

x(k+1)i = x

(k)i + ωr

(k+1)i = (1− ω)x(k)

i +ω

aii

(bi −i−1∑j=1

aijx(k+1)j −

n∑j=i+1

aijx(k)j ).

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čo-vat jako SOR (successive overrelaxaction).

Věta 2.8 K tomu, aby Gaussova-Seidelova relaxační metoda (2.13) kon-vergovala, je nutné, aby 0 < ω < 2.

54

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

n∏j=1

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

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∏j=1

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

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, . . . , npro i = 1, 2, . . . , m

x(k+1)i = (1− ω)x(k)

i + ωaii

(bi −∑i−1

j=1 aijx(k+1)j −∑n

j=i+1 aijx(k)j )

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čizobrazit přesně. Pokud úloha je špatně podmíněná, mohou zaokrouhlovacíchyby naprosto znehodnotit celý výpočet. Podívejme se proto nyní, jak jetomu v případě soustav lineárních rovnic.

OznačmeA, ~b — přesné hodnoty matice soustavy a vektoru pravých stran,

55

~x — k nim příslušný vektor řešení.

A, ~b — matice soustavy a vektoru pravých stran, které jsou zatížené chybou,~x — k nim spočtený vektor řešení.

Označme ∆A = A− A, ∆~b = ~b− ~b, ∆~x = ~x− ~x.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||||~x + ∆~x|| ≤ ||A−1|| ||A|| ||∆A||

||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

‖~b‖ = ‖A~x‖ ≤ ‖A‖ ‖~x‖,1Definice různých typů norem viz strana 18. Protože normy (1.8) - (1.10) jsou ekviva-

lentní, nezáleží na tom, se kterým typem normy (nadále pevně zvoleným) se rozhodnemepracovat.

56

a tedy

||~x|| ≥ ||~b||||A|| .

Dostáváme proto||∆~x||||~x|| ≤ ||A−1|| ||A|| ||∆

~b||||~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 tedyuzavřít tím, že pro řešení soustavy lineárních rovnic je číslo podmíněnostidefinováno vztahem (2.14).

57

58

Kapitola 3

Vlastní čísla a vlastní vektorymatice

V předchozím odstavci jsme viděli, jak je při studiu konvergence iteračníchmetod 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řebu-jeme 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 pujde 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á jejejich přibližná hodnota. Nesmíme ovšem zapomenout, že i když prvky maticeA jsou reálná čísla, mohou být vlastní čísla matice A komplexní. K odhadupolohy 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 kom-

plexní rovině ve sjednocení kruhů

|z − aii| ≤n∑

i=1,i6=j

|aij|.

59

Příklad Odhadněte polohu vlastních čísel a hodnotu nejmenšího vlastníhočísla matice

A =

4 1 01 2 10 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λ ∈ 〈0, 5〉.

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| ≤ n√

r1r2 . . . rn, kde ri = (n∑

j=1

|aij|2)1/2.

Příklad Odhadněte polohu vlastních čísel matice

1 1 00 2 00 0 3

.

Řešení: r1 =√

12 + 12 =√

2, r2 =√

22 = 2, r3 =√

32 = 3. Pro nejmenšívlastní číslo zadané matice platí |λn| ≤ 3

√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 dvazpusoby ř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∑

i=1

ai~vi,

~y (k) = A~y (k−1), k = 1, . . . , n.

Když existují limity

λ1,j = limk→∞

y(k+1)j

y(k)j

, j = 1, . . . , n,

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

λ1 =1n

n∑j=1

λ1,j.

Důkaz. Položme

~y (0) =n∑

i=1

ai~vi

61

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) = a1Ak~v1 + · · ·+ anAk~vn.

Protože λi jsou vlastní čísla matice A a ~vi k nim příslušné vlastní vektory, jeA~vi = λi~vi. Po dosazení do předchozího vztahu dostaneme, že

~y (k) = a1λk1~v1 + · · ·+ anλk

n~vn = λk1(a1~v1 + a2(

λ2

λ1)k~v2 · · ·+ an(

λn

λ1)k~vn).

Když označíme ~ε k =∑n

i=1 ai(λi

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

y(k)j = λk

1(a1v1,j + εkj ). (3.1)

Protože λ1 je podle předpokladu největší vlastní číslo, platí pro všechna i ≥ 2,že | λi

λ1| < 1, a tedy lim

k→∞εk

j = 0. Mužeme psát

λ1,j = limk→∞

y(k+1)j

y(k)j

= limk→∞

λk+11 (a1v1j + εk+1

j )

λk1(a1v1j + εk

j ), j = 1, 2, . . . , n.

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

∑nj=1 λ1,j.

Algoritmus (Mocninná metoda)

Vstup: n, A = aij, ~y(0), ε

λ01 = 0

pro k = 1, 2, . . . , dokud |λ(k)1 − λ

(k−1)1 | > ε

pro i = 1, 2, . . . , ny

(k)i =

∑nj=1 aijy

(k−1)j

λ(k)1,i = y

(k)i

y(k−1)i

λ(k)1 = 1

n

∑ni=1 λ1,i

Výstup: λ(k)1

62

Příklad Určete hodnotu dominantního vlastního čísla matice

1 2 00 0 13 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 , λ(1)1 = 5.3434,

~y(3) = (27, 37, 155)T , λ(2)1 = 4.3176,

~y(4) = (101, 155, 613)T , λ(3)1 = 3.9616,

~y(5) = (411, 613, 2459)T , λ(4)1 = 4.0119,

~y(6) = (1637, 2459, 9829)T , λ(5)1 = 3.9972,

~y(7) = (6555, 9829, 39323)T , λ(6)1 = 4.0007,

~y(8) = (26213, 39323, 157285)T , λ(7)1 = 3.9998,

~y(9) = (104859, 157285, 629147)T , λ(8)1 = 4.0000,

~y(10) = (419429, 629147, 2516581)T , λ(9)1 = 4.0000.

Dominantní vlastní číslo λ1 = 4.

Nevýhodou metody je, že předem nelze rozhodnout o tom, zda jsou splněnypožadované předpoklady, tj. zda |λ1| > |λ2| ≥ · · · ≥ |λn|. Také se těžkoodhaduje chyba aproximace.

b) Metoda Rayleighova podílu

Modifikací mocninné metody je metoda Rayleighova podílu, kterou lze apli-kovat 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) =n∑

i=1

ai~vi,

~y (k) = Ak~y (0), k = 1, . . . , n,

pak

λ1 = limk→∞

(~y (k))T~y (k+1)

(~y (k))T~y (k).

Důkaz. Použijeme stejných označení jako v důkazu věty 3.3. Především podle(3.1) položíme y

(k)j = λk

1(a1v1,j + εkj ). Když A je symetrická, pak vlastní

vektory jsou ortonormální a vektory ~y(k) splňují

(~y (k))T~y (k) = λ2k1 (a1~v

T1 + (~ε k)T )(a1~v1 + ~ε k) = λ2k

1 (a12 + (~ε k)T~ε k),

(~y (k))T A~y (k) = λk1(a1~v

T1 +(~ε k)T )λk+1

1 (a1~v1+~ε k+1) = (λk1)2k+1(a1

2+(~ε k)T~ε k+1).

Pro |λ1| > |λ2| ≥ · · · ≥ |λn| je limk→∞(~ε k)T~ε k = limk→∞(~ε k)T~ε k+1 = 0,a tedy

λ1 = limk→∞

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

(~y (k))T~y (k)= lim

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

(~y (k))T~y (k).

Algoritmus (Metoda Rayleighova podílu)

Vstup: n, A = (aij), ~y (0), ε

λ(0)1 = 0

pro k = 1, 2, . . . , dokud |λ(k)1 − λ

(k−1)1 | > ε

λ(k+1)1 = (~y k)T A~y k

(~y k)T ~y k

Výstup: λ(k+1)1

Příklad Najděte dominantní vlastní číslo matice A =

4 1 01 2 10 1 1

.

64

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

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

(~y (0))T ~y (0) = 3, 6667,

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

(~y (1))T ~y (1) = 4, 2667,

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

(~y (2))T ~y (2) = 4, 4086,

~y (4) = A~y (3) = (504, 252, 81)T , λ(4)1 = 4, 4472,

~y (5) = A~y (4) = (2268, 1089, 333)T , λ(5)1 = 4, 4571,

~y (6) = A~y (5) = (10161, 4779, 1422)T , λ(6)1 = 4, 4597,

~y (7) = A~y (6) = (45423, 21141, 6201)T , λ(7)1 = 4, 4603,

~y (8) = A~y (7) = (202833, 93906, 27342)T , λ(8)1 = 4, 4605,

~y (9) = A~y (8) = (905238, 417987, 121248)T , λ(9)1 = 4, 4605,

~y (10) = A~y (9) = (4038939, 1862460, 539235)T , λ(10)1 = 4, 4605.

Úplný problém vlastních čísel

Když chceme určit všechna vlastní čísla dané matice, mužeme zvolit jeden zedvou přístupu. Buď vlastní čísla určíme přímo jako kořeny charakteristickéhopolynomu p(λ) (definice charakteristického polynomu viz strana 15), nebovyuž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 cha-rakteristické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 nerealizova-telný. 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í matice tvaru

A =

a1 c1 0 . . . 0 0b2 a2 c2 . . . 0 00 0 0 . . . bn an

a položíme

f−1(λ) = 0,

f0(λ) = 1,

fk(λ) = (ak − λ)fk−1(λ)− bkck−1fk−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 polynomup(λ).

Algoritmus (Charakteristický polynom třídiagonální matice)

Vstup: n, A,

f−1(λ) = 0, f0(λ) = 1

pro k = 1, 2, . . . , nfk(λ) = (ak − λ)fk−1(λ)− bkck−1fk−2(λ)

Výstup: p(λ) = fn(λ)

Poznámka Tento postup je možné využít i pro obecnou čtvercovou maticipo té, co původní matici převedeme na matici třídiagonální, která má stejnávlastní čísla. (Viz [8] strana 542.)

Příklad Určete charakteristický polynom matice

2 −1 0−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.1podobné 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 = UL je podobná matici A.

Důkaz. Když A = LU, je U = L−1A. Odtud B = UL = L−1AL. 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 = L0U0.Záměnou pořadí násobení obou matic dostaneme novou matici A1 = U0L0.Provedeme LU-rozklad matice A1 : A1 = L1U1, a označíme A2 = U1L1 atd.Dostaneme posloupnost matic Ak, k = 0, 1, . . .Protože Ak = Uk−1Lk−1 a zároveň ze vztahu Ak−1 = Lk−1Uk−1 plyne Uk−1 =L−1

k−1Ak−1, dostáváme

Ak = Uk−1Lk−1 = L−1k−1Ak−1Lk−1 = L−1

k−1Uk−2Lk−2Lk−1 = . . .

= L−1k−1L

−1k−2Ak−2Lk−2Lk−1 = · · · = L−1

k−1 . . . L−10 A0L0 . . . Lk−1.

Označme Mk = L0L1 . . . Lk. Platí následující věta

Věta 3.7 Nechť A je regulární matice. Pokud posloupnost matic Mk kon-verguje 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,že Mk = L0L1 . . . Lk, dostaneme Lk = M−1

k−1Mk, a tedy

limk→∞

Lk = limk→∞

M−1k−1Mk = L−1L = I = lim

k→∞L−1

k .

Zároveň existuje

limk→∞

Uk−1 = limk→∞

AkL−1k−1 = lim

k→∞M−1

k−1A0Mk−1L−1k−1 = M−1AM = U.

Paklimk→∞

Ak = limk→∞

LkUk = limk→∞

Uk = U.

Algoritmus (Určení vlastních čísel pomocí LU-rozkladu)

Vstup: n, A = (ai,j), m

A0 := A, A0 = L0U0

pro k = 1, 2, . . . ,mAk := Uk−1Lk−1

LU-rozklad: Ak = LkUk

Výstup: Am = (amij ), λi ≈ am

ii

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í symet-rická pozitivně definitní matice, použijeme ke stanovení počátečního LU-rozkladu Choleského algoritmus. Označme Lk = (aij)i,j, Lk+1 = (bij)i,j,Ak+1 = Lk+1L

Tk+1. Prvky matice Ak+1 = B spočteme následovně:

Algoritmus (Určení vlastních čísel metodou LU-rozkladu pro symetrickoutřídiagonální matici)

Vstup: n, aii, ai+1,i, i = 1, 2, . . . n

pro i = 1, 2, . . . , nbii = a2

ii + a2i+1,i − b2

i−1,i

b2i+1,i = 1

b2iia2

i+1,ia2i+1,i+1

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

68

c) Určení vlastních čísel metodou ortogonálních trans-formací

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 q-tý sloupec

Qpq =

1...

.... . .

......

1...

.... . . . . . cos α . . . . . . . . . − sin α . . . . . . . . . . . .

... 1...

.... . .

...... 1

.... . . . . . sin α . . . . . . . . . cos α . . . . . . . . . . . .

...... 1

......

. . ....

... 1

je ortogonální matice stejného řádu, pak B = QTpqAQpq je reálná symetrická

matice, která je tvořena prvky

bqq = app sin2 α− apq sin 2α + aqq cos2 α,bpp = app cos2 α + apq sin 2α + aqq sin2 α,bpq = bqp = apq cos 2α + 1

2(aqq − app sin 2α),biq = bqi = −aip sin α + aiq cos α pro i 6= p, q,bip = bpi = aip cos α + aiq cos α pro i 6= p, q,bil = ail, i, l 6= p, q.

69

Poznámka Prvky bpq = bqp = 0, když tg 2α = 2apq

app−aqq, app 6= aqq, tj. když

α = π4 .

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 ma-tici 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í prveks největší absolutní hodnotou. To znamená, že ve stávající matici A vybe-reme prvek v tom řádku p a sloupci q, pro který apq = maxi,j{|aij|}.Hodnoty prvku nově vzniklé matice B spočteme ze vztahů

bpp =app + aqq + r

2,

bqq =appaqq − a2

pq

bpp

,

bpq = bqp = 0,

bip = aipc + aiqs = bpi, i = 1, 2 . . . , n, i 6= p, q,

biq = −aips + aiqc = bqi, i = 1, 2 . . . , n, i 6= p, q

bij = aij jinak.

V předchozích rovnicích jsme použili označení

r2 = (app + aqq)2 + 4apq,

s2 =12− app − aqq

2r,

c2 =12

+app − aqq

2r.

Vypočtené prvky bij matice B přeznačíme na aij. V nové matici vyhledámeten prvek, který má největší absolutní hodnotu a znovu podle uvedenýchvztahů 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(A))2 =n∑

i,j=1,i6=j

|aij|2.

Když apq = maxi,j a2i,j, je (N(A))2 ≤ (n2 − n)a2

pq. Pro největší nediagonálníprvek pak dostaneme

a2pq ≥

(N(A))2

n(n− 1).

Při této volbě je bpq = bqp = 0 a platí(N(B))2 = (N(A))2 − 2a2

pq ≤ q2(N(A))2,kde

q = (1− 2n(n−1))

12 .

Když v k-tém kroku přeznačíme B na {A(k)}, pak pro posloupnost podobnýchmatic {A(k)}, k ∈ N 1, platí

[N(A(k))]2 ≤ qk[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, ε

s =∑n

i=2

∑i−1j=1 a2

ij

dokud s > εurčíme indexy p, q, p 6= q, pro které ve stávající matici apq = max

i,j 6=p,q|ai,j|

podle vzahů na straně 70 spočteme prvky bij matice B.A := Bs :=

∑ni=2

∑i−1j=1 a2

ij

Výstup: λ1 ≈ a1,1, . . . , λn ≈ an,n

Příklad Jacobiovou metodou určete vlastní čísla matice A =

2 −1 0−1 2 −1

0 −1 2

.

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

71

Řešení: Postupně Jacobiovou metodou obdržíme

B =

1 0 0, 7071067810

0 3 0, 7071067810

0, 7071067810 0, 7071067810 2

,

B =

2, 366025404 0, 6279630301 0

0, 6279630301 3 0, 3250575835

0 0, 3250575835 0, 6339745962

,

B =

3, 386446078 0 0, 2768365599

0 1, 979579326 0, 1703641738

0, 2768365599 0, 1703641738 0, 6339745962

,

B =

3, 414013491 0, 01688138882 0

0, 01688138882 1, 979579326 0, 1695257220

0 0, 1695257220 0, 6064071831

,

B =

3, 414013491 0, 01675788864 −0, 002038248435

0, 01675788864 2, 000198603 0

−0, 002038248435 0 0, 5857879069

,

B =

3, 414212095 0 −0, 002038105311

0 1, 999999999 0, 00002415417276

−0, 002038105311 0, 00002415417276 0, 5857879069

,

B =

0, 585786439 −0, 00002415416649 0

−0, 00002415416649 1, 999999999 −0, 00000001740441820

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 budenají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ár-ní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í soustavynelineárních rovnic. Protože pro n > 1 jsou rovnice Pn(x) = 0 rovniceminelineárními, můžeme uvedené metody využít také ke stanovení kořenů po-lynomů. 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é α ∈ 〈a, b〉, pro které f(α) = 0.Abychom si zajistili existenci alespoň jednoho reálného řešení zadané neli-neární rovnice, budeme v celém tomto odstavci předpokládat, že funkce f(x)je spojitá na intervalu 〈a, b〉 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ýmzpů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 grafyfunkcí y = f1(x), y = f2(x), bude x-ová souřadnice průsečíku obou křivekodhadem ř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 hle-daný kořen.

y = x

x0 y = − ln x

Když jsme určili, v jakém intervalu budeme řešení hledat, popř. ze kteréhobodu x0 odstartujeme, budeme pokračovat některou vhodnou numerickoumetodou. Numerické metody, o kterých dále budeme hovořit, rozdělme dodvou skupin, a to na startovací a zpřesňující.Mezi startovací zařadíme ty metody, které konvergují, avšak rychost konver-gence bývá pomalá. Patří sem například metoda bisekce, regula falsi neboprostá 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í, takrychleji než metody startovací. Proto slouží hlavně ke zpřesnění řešení.

74

a) Metoda bisekce

Postupným půlením intervalu 〈a, b〉 sestrojíme posloupnost intervalu, kteréobsahují hledaný kořen. Střed s intervalu 〈a, b〉 bereme za první aproximacikořene. Pokud f(a)f(s) < 0, přeznačíme s na b, jinak přeznačíme s na a.Nový interval 〈a, b〉 opět rozpulí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 b−a2k−1 > ε

s := b+a2

xk := skdyž f(a) f(s) < 0, položíme a := s jinak b := skdyž 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

0

0.2

0.4

0.6

0.8

1

1 1.2 1.4 1.6 1.8 2x

y = arctg(x)

y = 1x

Vidíme, že kořen α leží uvnitř íntervalu 〈1, 2〉. Metodou bisekce budemepočítat tak dlouho, než chyba bude menší než 10−2.

a = 1 b = 2 s1 = 1, 5,a = 1 b = 1, 5 s2 = 1, 25,a = 1 b = 1, 25 s3 = 1, 125,a = 1, 125 b = 1, 25 s4 = 1, 1875,a = 1, 125 b = 1, 1875 s5 = 1, 15625,a = 1, 15625 b = 1, 1875 s6 = 1, 171875,a = 1, 15625 b = 1, 171875 s7 = 1, 1640625

.

Chyba |α− s7| ≤ 127 = 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í formulixk = ϕ(xk−1). (4.1)

Zadáme počáteční hodnotu x0 a přesnost δ, s jakou chceme počítat. Iteru-jeme 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 vubec konver-govat. Uveďme si nyní větu, která nám umožní určit podmínky, za jakýchmetoda prosté iterace konverguje.

Věta 4.1 Nechť ϕ ∈ C1〈a, b〉, supx∈〈a,b〉

|ϕ′(x)| < 1, pak rovnice x = ϕ(x) má

jediné řešení α ∈ 〈a, b〉, ke kterému konverguje posloupnost iterací

xk+1 = ϕ(xk), k = 0, 1, . . . ,

pro libovolné x0 ∈ 〈a, b〉.Pro všechna k ∈ N platí apriorní odhad chyby

|xk − α| ≤ qk

1− q|x1 − x0|

a aposteriorní odhad chyby

|xk − α| ≤ q

1− q|xk − xk−1|.

Důkaz. Když ϕ ∈ C1〈a, b〉, pak podle věty o střední hodnotě (viz str. 112 veskriptech [4]) pro libovolné dva body x, y ∈ 〈a, b〉, x < y, platí

ϕ(x)− ϕ(y) = ϕ′(ξ)(x− y), ξ ∈ (a, b).

Odtud a z předpokladů věty dostaneme

|ϕ(x)− ϕ(y)| ≤ supξ∈(a,b)

|ϕ′(ξ)| |x− y| = q|x− y|.

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í ϕ ∈ C1〈a, b〉, |ϕ′(x)| < 1, pakmetoda 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 〈a, b〉. Pak platí

Věta 4.3 Pokud ϕ ∈ C0〈a, b〉 a existuje q ∈ 〈0, 1) tak, že

|ϕ(x)− ϕ(y)| ≤ q|x− y|, x, y ∈ 〈a, b〉, (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 in-tervalu 〈a, b〉.

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 〈0, 1; 1〉.Jednou z možností, jak vybrat iterační funkci, je položit ϕ(x) = − ln(x).Protože

q = supx∈〈0,1;1〉

|ϕ′(x)| = supx∈〈0,1;1〉

|1x| ≥ 1,

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

supx∈〈0,1;1〉

|ϕ′(x)| = supx∈〈0,1;1〉

e−x = 0.3678794412 = q < 1,

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

x1 = 0, 5664147332, x2 = 0, 5675566373,x3 = 0, 5669089119, x4 = 0, 5672762322,x5 = 0, 5670678984, x6 = 0, 5670678984,x7 = 0, 5671860501, x8 = 0, 5671190401.

Chyba aproximace:

|x− x8| ≤ q

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

V případě iteračních metod není zanedbatelná otázka rychlosti konvergence.Rychlost konvergence iterační metody závisí řádu iterační metody. Metodyvyšší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í

limk→∞

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

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 ∈ 〈1, r), ϕ(r)(α) 6= 0.

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

79

Důkaz. Z Taylorova rozvoje funkce ϕ v bodě α máme

xk − α = ϕ′(α)(xk−1 − α) +ϕ′′(ξ)

2(xk−1 − α)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)

b− a(x− 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. Zkon-struujme novou sečnu. Proces zastavíme, když funkční hodnota v průsečíkubude 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) > δ,

p := a− f(a)f(b)−f(a)(b− 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, žemetoda je řádu 1.

Poznámka Pokud je f ∈ C1〈a, b〉, jsou splněny předpoklady Lagrangeovyvěty (viz skripta [4] str. 77) a existuje tedy ξ ∈ (xk, α) ⊂ 〈a, b〉, že

f ′(ξ) =f(xk)− f(α)

xk − α.

81

Odtud, nezávisle na konkrétním tvaru metody, dostaneme odhad chyby

|xk − α| ≤ f(xk)m1

, kde m1 = minx∈(xk,α)

|f ′(x)|. (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 in-tervalu 〈1, 2〉.

a = 1 b = 2 p1 = 1, 150186819a = 1, 150186819 b = 2 p2 = 1, 161536517a = 1, 161536517 b = 2 p3 = 1, 162287210a = 1, 162287210 b = 2 p4 = 1, 162336388a = 1, 162336388 b = 2 p5 = 1, 162339607a = 1, 162339607 b = 2 p6 = 1, 162339818a = 1, 162339818 b = 2 p7 = 1, 162339832

Chyba |α− p7| ≤ |f(p7)m1

| = |10−9

π4−1 | = 0, 466 · 10−8.

d) Newtonova metoda

Tato metoda byla objevena roku 1669. Newton ji použil při řešení kubickýchrovnic a při aproximaci řešení Keplerových rovnic pro pohyb planet.Pokud f ∈ C1〈a, b〉, můžeme levou stranu rovnice f(x) = 0 nahradit Taylo-rovým polynomem 1. stupně. Obdržíme

f(x0) + f ′(x0)(x− x0) = 0.

Když odtud spočteme x a připojíme indexy, obdržíme Newtonovu iteračníformuli

xk = xk−1 − f(xk−1)f ′(xk−1)

. (4.5)

Iterační proces zastavíme podmínkou |xk − xk−1| ≤ δ pro zvolené δ > 0.

82

Definice 4.2 Nechť f ∈ C1〈a, b〉, f ′(x) 6= 0, x ∈ 〈a, b〉, pak Newtonovametoda pro řešení rovnice f(x) = 0 je dána iteračním schématem

xk+1 = xk − f(xk)f ′(xk)

, k = 0, 1, . . . ,

pro počáteční hodnotu x0 ∈ 〈a, b〉.

Algoritmus: (Newtonova metoda)

Vstup: a, b, f(x), x0, δ > 0

pro k = 1, 2, . . . , dokud |xk − xk−1| > δf ′(xk−1)xk = xk−1 − f(xk−1)

f ′(xk−1)

Výstup: xk

Poznámka Newtonova metoda bývá někdy nazývána metodou tečen. Toproto, ž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í

|xn − α| ≤ 1c|c(x0 − α)|2k

, (4.6)

kde c ≥ 12

f ′′(x1)f ′(x2) pro libovolné x1, x2 ∈ (a, b). Je vidět, že metoda bude určitě

konvergovat, když |c(x0 − α)| < 1, tj. když počáteční aproximace x0 je dostblí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 kon-verguje.

Věta 4.5 Když f ∈ C2(a, b) a kořen α rovnice f(x) = 0 je jednoduchý (tj,f ′(α) 6= 0), pak Newtonova metoda je lokálně konvergentní.

Věta 4.6 Jestliže f ′(x) 6= 0, x ∈ 〈a, b〉 a f ′′(x) v intervalu 〈a, b〉 neměníznaménko, f(a)f(b) < 0, f(a)

f ′(a) < b−a a f(b)f ′(b) < b−a, pak Newtonova metoda

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 ′(x) = 1 + e−x.

x0 = 1 x1 = 0, 5378828427 x2 = 0, 5669869914x3 = 0, 5671432860 x4 = 0, 5671432904.

Když tyto výsledky porovnáme s výsledky příkladu na str. 78, vidíme, žerychlost Newtonovy metody je podstatně vyšší než u metody prosté iterace.

Poznámka Newtonova metoda je metodou řádu 2. Kvadratická konvergenceNewtonovy metody znamená, že počet platných číslic v numerické aproximacise 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ímelineárním polynomem v bodech xk−1, xk. Z nové rovnice

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

(x− xk) + f(xk) = 0

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 ′(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á seuká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 − f(xk) xk−xk−1

f(xk)−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 typu nelineárních algebraickch rovnicPn(x) = 0. Patří sem rovnice

kvadratické a2x2 + a1x + a0 = 0,

kubické a3x3 + a2x

2 + a1x + a0 = 0,

binomické xn − a = 0,

trinomické ax2n + bxn + c = 0,

reciproké a0xn + a1x

n−1 + · · · = 0, kde ak = an−k nebo ak = −an−k.

U zbývajících typu určujeme kořeny numericky. K tomu použijeme buď me-tody, 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) = a0xn + a1x

n−1 + · · ·+ an−1x + 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|}

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

|a0| .

Věta 4.8 (Budan-Fourier) Nechť a, b ∈ R, a < b, koeficient an polynomuPn(x) je kladný a Pn(a) ·Pn(b) 6= 0. Když V (a) je počet znaménkových změnv posloupnosti Pn(a), P ′

n(a), . . . , P(n)n (a) a V (b) je počet znaménkových změn

v posloupnosti Pn(b), P ′n(b), . . . , P

(n)n (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 posloupnostikoeficientů polynomu Pn(x) resp. Pn(−x).

Věta 4.9 (Descartes) Počet kladných reálných kořenů polynomu

Pn(x) = a0xn + a1x

n−1 + · · ·+ an−1x + an

(včetně jejich násobnosti) je roven počtu znaménkových změn v posloupnostikoeficientů a0, a1, . . . , an nebo o sudý počet méně. (Nulové koeficienty neu-važujeme).Počet záporných kořenů je roven počtu znaménkových změn v posloupnosti(−1)na0, (−1)n−1a1, . . . , 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

13

= (1 +max{1, 2}

1)−1 < |α| < 1 +

max{2, 1}1

= 3

x ∈ (−3,−13

) ∪ (13, 3).

Podle věty 4.9 má uvažovaný polynom 2 nebo žádný reálný kladný kořen a1 reálný záporný kořen.Dále P ′(x) = 3x2 − 2, P ′′(x) = 6x, P ′′′(x) = 6,

x 0 1 2 3P (x) + 0 + +P ′(x) − + + +P ′′(x) 0 + + +P ′′′(x) + + + +V (x) 2 0 0 0

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

87

Pokud známe některý kořen polynomu Pn(x) = a0xn + a1x

n−1 + · · ·+ an, jemožné snížit stupeň polynomu prostřednictvím Hornerova schématu. Jestližepříslušný kořen označíme α1, dostaneme

Pn(x) = (x− α1)Pn−1(x), kde Pn−1(x) = b0xn−1 + b1x

n−2 + · · ·+ bn−1.

Pokud α2 je kořenem polynomu Pn−1(x), můžeme znovu snížit stupeň poly-nomu pomocí Hornerova schématu atd. Musíme si však uvědomit, že běhemvý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 po-lynomu 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, . . . , nbk = αbk−1 + ak

Výstup: b1, b2, . . . , bn

Snížení stupně polynomu je možné dosáhnout také tak, že vydělíme kvadra-tickým činitelem d(x) = x2 − px− q. Potom

Pn(x) = d(x)Qn−2(x) + ϕ(x),

kde Qn−2(x) = b0xn−2 + b1x

n−1 + · · ·+ bn−2, ϕ(x) = bn−1x + bn.

Algoritmus: (Dělení kvadratickým polynomem)

Vstup: n, p, q, a0, . . . , an

b−1 = 0, b0 = a0

pro k = 1, 2, . . . , n− 1bk = 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) = a0xn + a1x

n−1 + · · · + an−1x + an máreálné různé kořeny α1 > α2 > · · · > αn a

yk =−1a0

(a1yk−1 + a2yk−2 + · · ·+ anyk−n)

je řešením diferenční rovnice

a0yk + a1yk−1 + a2yk−2 + · · ·+ anyk−n = 0

doplněné podmínkami

y0 = y1 = · · · = yn−2 = 0, yn−1 = 1,

pakα1 = lim

k→∞yk+1

yk

. (4.8)

Důkaz. Algebraická rovnice a0xn + a1x

n−1 + · · ·+ an−1x + an = 0 je charak-teristickou rovnicí příslušnou k diferenční rovnici

a0yk + a1yk−1 + a2yk−2 + · · ·+ anyk−n = 0.

Když α1, . . . , αn jsou reálné různé kořeny charakteristické rovnice, má řešenídiferenční rovnice tvar yk =

∑ni=0 ciα

ki . (Konstanty ci vypočteme z počáteč-

ních podmínek.) Pro c1α1 6= 0 je

yk = c1αk1(1 +

n∑i=2

ci

c1(αi

α1)k).

Za předpokladu, že α1 je největší vlastní číslo a i > 1, je limk→∞( αi

α1)k = 0.

Pak

limk→∞

yk+1

yk

= limk→∞

c1αk+11

c1αk1

= α1.

89

Algoritmus: (Bernoulliova metoda)

Vstup: n, a0, . . . , an, my0 = y1 = · · · = yn−2 = 0, yn−1 = 1pro k = n, n + 1, . . . , m

yk = − 1a0

∑nj=1 aiyk−i

α1 = yk

yk−1

Výstup: α1

Příklad Najděte dominantní kořen rovnice x3 − x2 − x4 + 1

4 = 0.

Řešení: Položíme y0 = y1 = 0, y2 = 1 a pro k = 4, 5, . . . spočítáme iteraceyk = yk−1 + 1

4yk−2 − 14yk−3.

y3 = 1 α11 = 1

y4 = 1, 250000000, α21 = 1, 25,

y5 = 1, 250000000, α31 = 1,

y6 = 1, 312500000, α41 = 1, 05,

y7 = 1, 312500000, α51 = 1,

y8 = 1, 328125000, α61 = 1, 011904762,

y9 = 1, 328125000, α71 = 1.

Dominantní kořen α1 = limn→∞

αn1 = 1.

Poznámka Existuje i modifikace Bernoulliovy metody, kterou lze použítk 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 koefi-cienty kvadratické rovnice. Když rovnice má tvar a0x

2 + a1x + a2 = 0, pakpro její kořeny α1, α2 platí

α1 + α2 = −a1

a0, α1α2 =

a2

a0. (4.9)

90

První z těchto rovnic lze přepsat na tvar

α1(1 +α2

α1) =

−a1

a0.

Když předpokládáme, že |α2α1| << 1, můžeme kořeny zadané kvadratické rov-

nice přibližně spočítat pomocí vztahů

α1 ≈ −a1

a0, α2 ≈ a2

a0(−a0

a1) =

−a2

a1.

Sestavme nyní novou kvadratickou rovnici b0x2 + b1x + b2 = 0, jejíž kořeny

jsou kvadrátem kořenů α1, α2. Platí

β1 + β2 = α21 + α2

2 = (α1 + α2)2 − 2α1α2 =a2

1

a20

− 2a2

a0=

a21 − 2a2a0

a20

= −b1

b0,

β1β2 = α21α

22 =

a22

a20

=b2

b0.

Koeficienty nové kvadratické rovnice proto splňují

b0 = a20,

b1 = −(a21 − 2a0a2),

b2 = a22

a pro kořeny platí

β1 ≈ −b1

b0, β2 ≈ b2

b1.

Kořeny původní rovnice

|α1| ≈√|b1||b0| , |α2| ≈

√|b2||b1| .

Této úvahy lze využít při hledání reálných kořenů polynomu

Pn(z) = a0xn + a1x

n−1 + · · ·+ an,

pokud jsou tyto kořeny uspořádány sestupně podle absolutních hodnot jejichvelikostí. Definujme posloupnost polynomů P 0(x), P 1(x), . . . , P n(x) tak, že

P 0(x) = Pn(x), P k(x) = ak0x

n + ak1x

n−1 + · · ·+ akn−1x + ak

n,

91

kde

a0i = ai, i = 0, 1, . . . , n,

a pro k = 1, 2, . . . , n položme

ak0 = (ak−1

0 )2,

−ak1 = (ak−1

1 )2 − 2ak−10 ak−1

2 ,

ak2 = (ak−1

2 )2 − 2ak−11 ak−1

3 + 2ak−10 ak−1

4 ,

−ak3 = (ak−1

3 )2 − 2ak−12 ak−1

4 + 2ak−11 ak−1

5 − 2ak−10 ak−1

6 ,

. . . . . .

(−1)n−1akn−1 = (ak−1

n−1)2 − 2ak−1n−2a

k−1n ,

(−1)nakn = (ak−1

n )2.

Pro kořeny původní rovnice dostáváme

|αj| ≈ 2k

√| ak

j

akj−1

|, j = 1, 2, . . . , n. (4.10)

Algoritmus: (Graefova metoda)

Vstup: n, a0, . . . , an, m

pro j = 1, 2, . . . , n

a0j = aj

pro k = 1, 2, . . . , m

akj = (−1)j[(ak−1

j )2 + 2∑j−1

s=1(−1)s(ak−1s−1)ak−1

s+1

|αj| ≈ 2k

√| ak

j

akj−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 ak0 ak

1 ak2

0 1 −1, 5 0, 51 1 −1, 25 0, 252 1 −1, 0625 0, 0625

Z hodnot v posledním řádku (viz vztahy (4.10)) dostaneme odhad kořenů

|α1| .= 4

√|a

21

a20

| = 4√

1, 0625 = 1, 015915975,

|α2| .= 4

√|a

22

a21

| = 4

√0, 06251, 0625

= 0, 4924790605.

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é ruzné kořeny. Uvedenámetoda se ale dá modifikovat tak, že ji lze použít i ke stanovení páru kořenukomplexních.

c) Laguerrova metoda

Je to hojně používaná iterační metoda, která slouží k výpočtu kořenů alge-braické rovnice Pn(x) = 0. Metoda je velice univerzální, protože konvergujejak v případě reálných, tak v případě komplexních kořenů. Konvergence jezvlášť rychlá, když kořeny jsou reálné jednoduché.Předpokládejme, že α1 ≤ α2 ≤ · · · ≤ αn jsou kořeny algebraické rovnicePn(x) = 0. Počáteční aproximaci x0 ∈ R zvolme libovolně. Další aproximacekořene určíme ze vztahu

xk+1 = xk − nP (xk)

P ′(xk)±√

H(xk), (4.11)

H(xk) = (n− 1)[(n− 1)(P ′(xk))2 − nP (xk)P ′′(xk)].

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, mpro k = 1, 2, . . . ,m

H(xk) = (n− 1)[(n− 1)(P ′(xk))2 − nP (xk)P ′′(xk)]xk+1 = xk − nP (xk)

P ′(xk)±√

H(xk)

Výstup: xm

Příklad Prostřednictím Laguerrovy metody najděte některý z kořenů po-lynomu x3 − x2 − x

4 + 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

–1.4

–1.2

–1

–0.8

–0.6

–0.4

–0.2

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

y

–1–0.8 –0.4 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2x

p(x) = x3 − x2 − x4 + 1

4

94

Kapitola 5

Aproximace funkcí

Častokrát máme co do činění s funkcemi, které jsme získali jako výsledkylaboratorních měření a u kterých proto máme k dispozici jen jejich hodnotyv 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 funkcif 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 ∈ 〈a, b〉 a n + 1 funkčních hodnotf(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 minimalizu-jeme součet kvadrátů odchylek funkce f od funkce ϕ. Tzn. minimalizujeme

n∑i=0

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

• Čebyševova (stejnoměrná) aproximace: Žádáme, aby

maxxi∈〈a,b〉

|(f(xi)− ϕ(xi))| (5.3)

1Splně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átunormy v prostoru l2〈a, b〉 a normy v prostoru C0〈a, b〉.

Interpolace

Už několik století se provádí interpolace pomocí polynomů. Poprvé bylo in-terplačních polynomů použito v 17. století při tabelování logaritmů. Proč seprávě polynomy těší takové oblibě, vyplývá z jejich vlastností. Polynom jefunkce, 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 jemož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) ∈ C0〈a, b〉, pak ke každému ε > 0 existujen a polynom n-tého stupně p(x), tak, že

maxx∈〈a,b〉

|f(x)− p(x)| < ε.

Nechť tedy funkce, kterou chceme aproximovat, je zadaná funkčními hodno-tami 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čmePn množinu všech polynomů n-tého stupně. Bázi tohoto prostoru tvoří poly-nomy {1, x, x2, . . . , xn}. Naším cílem bude zkonstruovat interpolační polynomp(x) ∈ Pn.

96

a) Lagrangeův interpolační polynom

Nejprve si uveďme postup, který při konstrukci interpolačního polynomuv roce 1794 použil francouzský matematik J.L.Lagrange.

Věta 5.2 (Lagrangeova konstrukce interpolačního polynomu) Když je dánon+1 uzlových bodů x0, x1, . . . , xn ∈ 〈a, b〉 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∑

i=0

li(x)fi, (5.4)

kde

li(x) =n∏

j=0,j 6=i

(x− xj)(xi − xj)

. (5.5)

Důkaz. Existence: Především je vidět, že pro všechna i = 1, . . . , n jsou li(x) ∈Pn. Dále pro i, j = 1, . . . , n platí li(xj) = δij. Pak ale L(xj) =

∑i δijfi = fj.

Tzn. polynom L splňuje interpolační podmínky (5.1).Jednoznačnost: Předpokládejme, že existují dva různé interpolační polynomyL1, L2 ∈ Pn, L1 6= L2. Pak polynom L = L1 − L2 splňuje v uzlových bodechpodmínky

L(xj) = 0, j = 0, . . . , n.

Je tedy L polynom stupně n, který má n + 1 kořenů. To může nastat pouzetehdy, když všechny jeho koeficienty jsou nulové. To znamená, že L = L1 −L2 = 0, a tedy L1 = L2.

Algoritmus: (Lagrangeuv interpolační polynom)Vstup: n, x0, x1, . . . , xn, f0, f1, . . . , fn

pro i = 0, 1, . . . , nli(x) =

∏nj=0,i6=j

(x−xj)(xi−xj)

L(x) =∑n+1

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 3fi 2 1 2 0

Řešení:

l0(x) =x(x− 1)(x− 3)−1(−2)(−4)

=x3 − 4x2 + 3x

−8,

l1(x) =(x + 1)(x− 1)(x− 3)

1(−1)(−3)=

x3 − 3x2 − x + 33

,

l2(x) =(x + 1)x(x− 3)

2.1(−2)=

x3 − 2x2 − 3x

−4,

l3(x) =(x + 1)(x− 1)x

4.3.2=

x3 − x

24.

L(x) = 2l0(x) + l1(x) + 2l2(x) =112

(−5x3 + 12x2 + 5x + 12).

Poznámka Při označení

qn+1(x) =n∏

j=0

(x− xj)

lze psát Lagrangeovy polynomy li ve tvaru

li(x) =qn+1(x)

(x− xi)q′n+1(xi).

Věta 5.3 (Odhad chyby) Když f : 〈a, b〉 → R, f ∈ Cn+1(a, b) a L jeinterpolační polynom zkonstruovaný v bodech x0, . . . , xn ∈ 〈a, b〉, pak provšechna x ∈ 〈a, b〉 je

f(x)− L(x) =f (n+1)(ξ)(n + 1)!

n∏j=0

(x− xj), kde ξ ∈ 〈a, b〉. (5.6)

98

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

qn+1(x) =n∏

j=0

(x− xj).

Pro pevné x, x 6= xi, i = 0, 1, . . . , n, uvažujme funkci g : 〈a, b〉 → R takovou,že

g(y) = f(y)− L(y)− qn+1(y)f(x)− L(x)

qn+1(x), y ∈ 〈a, b〉.

Protože podle předpokladů věty f ∈ Cn+1(a, b), je také g ∈ Cn+1(a, b). Toale znamená, že g musí mít alespoň n + 1 kořenů, g′ má alespoň n kořenů ag(n+1) alespoň jeden kořen ξ ∈ 〈a, b〉. 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, odhadujemechybu interpolačního vzorce vztahem

||f(x)− L(x)||∞ ≤ 1(n + 1)!

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

Z tohoto vztahu je vidět, že chyba interpolace závisí na vlastnostech interpo-lované funkce. Protože vztah (5.7) obsahuje výraz qn+1, je chyba interpolaceovlivněna také rozložením uzlů.

Poznámka Kromě chyby metody může být výpočet ovlivněn ještě zaokrouh-lovacími chybami. Pokud jsou hodnoty v uzlech xi zatíženy chybami, pak jehodnota aproximace

fi = fi + ∆fi, i = 0, . . . , n,

a Lagrangeův interpolační polynom

n∑i=0

li(x)fi =n∑

i=0

li(x)fi +n∑

i=0

li(x)∆fi.

99

Dále

|n∑

i=0

li(x)∆fi| ≤n∑

i=0

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

n∑i=0

|li(x)|,

kde c ≥ |∆fi|, i = 0, 1, . . . , n. Protože∑n

i=0 li(x) = 1, je∑n

i=0 |li(x)| ≥ 1 amůžeme učinit zá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 〈a, b〉rozloženy ekvidistantně (tj. jsou od sebe stejně vzdálené) s krokem h. Prouzlové body pak máme

xi = x0 + ih, i = 0, 1, . . . , n.

Pokud označíme relativní vzdálenost bodu x od bodu x0 symbolem t,

t =x− x0

h,

pak x = x0 + th a pro koeficienty Lagrangeova polynomu (viz vzorec (5.5))dostaneme

li(x0 + th) =n∏

j=0,i6=j

(x− xj)(xi − xj)

=n∏

j=0,i6=j

(x0 + th− x0 − jh)(x0 + ih− x0 − jh)

=n∏

j=0,i 6=j

(t− j)(i− j)

= li(t).

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∑

i=0

li(t)fi = L(t)

je pak také funkcí závisou na t. Podle (5.6) je rozdíl

f(x0 + th)− L(x0 + th) =f (n+1)(ξ)(n + 1)!

n∏j=0

(t− j)h =f (n+1)(ξ)(n + 1)!

hn+1n∏

j=0

(t− j).

(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 zkou-má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ímeznovu 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žímesoustavu 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

a0 = f0,

a1 =f1 − a0

x1 − x0=

f1 − f0

x1 − x0,

a2 =f2 − a0 − a1(x2 − x0)

(x2 − x0)(x2 − x1)=

[f2 − a0 − a1(x1 − x0)] + [a1(x1 − x0)− a1(x2 − x0)](x2 − x0)(x2 − x1)

=

=1

x2 − x0(f2 − f1

x2 − x1+

f1−f0x1−x0

(x1 − x2)

x2 − x1) =

1x2 − x0

(f2 − f1

x2 − x1− f1 − f0

x1 − x0),

. . . . . .

Pokud označíme D0i = fi, i = 0, . . . , n, a

Dki =

Dk−1i −Dk−1

i−1

xi − xi−k

, k = 1, 2, . . . n, i = k, . . . , n, (5.10)

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

N(x) = D00 +

n∑

k=1

Dkk

k−1∏i=0

(x− xi). (5.11)

101

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

x0 D00 = f(x0)

x1 D01 = f(x1) D1

1...

......

......

...

xn−1 D0n−1 = f(xn−1) D1

n−1... Dn−1

n−1

xn D0n = f(xn) D1

n

... Dn−1n Dn

n

Věta 5.4 V Newtonově reprezentaci má interpolační polynom z věty 5.2tvar (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 Lagrange-ova 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, . . . , nD0

i = fi

pro k = 1, 2, . . . , npro i = k, k + 1, . . . , n

Dki =

Dk−1i −Dk−1

i−1

xi−xi−k

Výstup: a0 = D00, a1 = D1

1, . . . , 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řiurčování Newtonova interpolačního polynomu je proto, stejně jako v Horne-rově algoritmu, zapotřebí O(n2) operací.

102

Příklad Určete Newtonuv interpolační polynom, když je funkce zadánatabulkou

xi −1 0 1 3fi 2 1 2 0

Řešení:−1 2

0 1 −11 2 1 13 0 −1 −2

3 − 512

N(x) = 2− (x + 1) + (x− 0)(x + 1)− 512(x− 0)(x + 1)(x− 1) =

= − 512x

3 + x2 + 512x + 1.

Poznámka Jistou analogií uvedené Newtonovy konstrukce je tzv. Nevilluvalgoritmus, který umožňuje stanovit funkční hodnotu interpolačního poly-nomu v bodě α, aniž bychom znali konkrétní tvar interpolačního polynomu.Označme pk

i 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,

p1i (x) =

(x− xi)p0i+1(x)− (x− xi+1)p0

i (x)xi+1 − xi

=

=[(x− xi+1) + (xi+1 − xi)]p0

i+1(x)− (x− xi+1)p0i (x)

xi+1 − xi

=

= p0i + (x− xi)

p0i (x)− p0

i−1(x)xi − xi−1

, i = 0, . . . , n− 1,

p2i (x) =

(x− xi)p1i+1(x)− (x− xi+2)p1

i (x)xi+2 − xi

=

= p1i (x) + (x− xi)

p1i (x)− p1

i−1(x)xi − xi−1

, i = 0, . . . , n− 2,

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

pki (x) = pk−1

i (x) + (x− xi)pk−1

i (x)− pk−1i−1 (x)

xi − xi−1, i = k, . . . , n.

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, . . . , np0

i = fi

pro k = 1, 2, . . . , npro i = k, k + 1, . . . , n

pki = pk−1

i + (α− xi)pk−1

i −pk−1i−1

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 3fi 2 1 2 0

Řešení: Výpočet provedeme v tabulce

xi α− xi fi p1i p2

i p3i

−1 3 20 2 1 −11 1 2 3 5

3 −1 0 1 5/3 5/2

Zde

p11 = 1 + 2

1− 20− (−1)

= −1,

p12 = 2 + 1

2− 11− (0)

= 3,

p13 = 0− 1

0− 23− 1

= 1,

p22 = 3 + 1

3 + 11− (−1)

= 5, atd.

104

Hodnota interpolačního polynomu p(2) = 52 .

Při ekvidistantních uzlech můžeme Newtonův interpolační polynom zapsatve tvaru, v němž vystupují diference vpřed a vzad.

V uzlech xi = x0 + ih, i = 0, . . . , n, definujeme diference vpřed následovně

∆0fi = fi,

∆1fi = fi+1 − fi,

∆2fi = ∆(fi+1 − fi) = fi+2 − 2fi+1 + fi,

. . . . . . . . .

∆kfi = ∆(∆k−1fi) =k∑

j=0

(kj

)fi+k−j.

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

Dki =

∆kfi

k!hk,

můžeme psát, že Newtonův interpolační polynom

N+(x0 + th) = f0 +

(t

1

)∆f0 + · · ·+

(t

n

)∆nf0 =

n∑i=0

(t

i

)∆if0. (5.12)

V případě, že vycházíme při konstrukci interpolačního polynomu z uzlů napravém konci interpolačního intervalu, můžeme využít diferencí vzad. Tytodiference definujeme vztahy

∇0fi = fi,

∇1fi = fi − fi−1,

∇2fi = ∇(fi − fi−1) = fi − 2fi−1 + fi−2,

. . . . . . . . .

∇kfi = ∇(∇k−1fi) =k∑

j+0

(kj

)fi−k+j.

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í:

∆kfi−k = ∇kfi,

můžeme Newtonův interpolační polynom psát ve tvaru:

N−(x) = fn +

(−t

1

)∇fn + · · ·+

(−t

n

)∇nf0 =

n∑i=0

(−1)i

(−t

i

)∇if0.

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 pouzev těch případech, kdy nejsme příliš vzdáleni od koncových bodu zadanéhointervalu.

Příklad Když znáte hodnoty funkce sin x pro x = 20◦ a 21◦, určete (inter-polací) sin 20◦18′ a (extrapolací) sin 21◦18′.

xi 20 21sin xi 0, 34202 0, 35837

Řešení: Při interpolaci máme

f(x) ∼ L1(x) =x− x1

x0 − x1f0 +

x− x0

x1 − x0f1,

sin 20◦18′ ∼ 20, 3− 2120− 21

0, 34202 +20, 3− 2021− 20

0, 35837 = 0, 34693.

Chyba

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

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

106

V našem případě je n = 1, f ′′(x) = − sin x a

M2 = maxx∈〈20,21〉

| − sin x| = 0, 35837.

To znamená, že

E(x) ≤ 0, 358372

(0, 3π

1800, 7π

180) = 1, 1 · 10−5.

Extrapolací dostaneme

sin 21◦18′ ∼ 21, 3− 2120− 21

0, 34202 +21, 3− 2021− 20

0, 35837 = 0, 36328,

E(x) ≤ 0, 358372

(0, 3π

1801, 7π

180) = 3 · 10−5.

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žíttzv. Richardsonovy extrapolace.Označme veličinu, kterou počítáme, F a T (h) výsledek numerického výpočtunějakou metodou. Předpokládejme, že při zjemňování kroku h dostávámepřesnější řešení,

limh→0+

T (h) = F,

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 zkon-struovaným v bodech Hp, (qH)p,

L1(h) =h− (qH)p

Hp − (qH)pT (H)− h− (H)p

(qH)p −HpT (qH).

Richardsonovou extrapolací budeme rozumět hodnotu tohoto polynomu v boděh = 0. 2

L1(0) =−(qH)p

Hp − (qH)pT (H) +

(H)p

(qH)p −HpT (qH) =

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

=

= T (H) +T (H)− T (qH)

qp − 1.

2O 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 chybyaproximace je v tomto případě roven hr, r > p.Zobecněním tohoto postupu je tzv. opakovaná Richardsonova extrapolace. Jemožné ji realizovat, pokud numerická metoda T (h) umožňuje vyjádřit chybuaproximace ve formě polynomu s-tého stupně v mocninách hp, tj. ve tvaru

T (h)− F ≈ c1hp + c2h

2p + · · ·+ cshsp.

Pak můžeme ke konstrukci nové formule použít interpolační polynom s-téhostupně. Jeho hodnota v bodě h = 0 bude představovat přesnější aproximaciveličiny F. K výpočtu hodnoty interpolačního polynomu Ls(0) je možné po-uží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 kte-rou interpolační proces nekonverguje stejnoměrně. Speciálně problematickáje interpolace spojitých funkcí při použití ekvidistantních uzlů. (Např. profunkci y = |x| konverguje interpolační proces s ekvidistantními uzly pouzev bodech x = −1, 0, 1 intervalu 〈−1, 1〉.)Navíc polynomy vyšších stupňů se na koncích více „vlníÿ a to nemusí býtv ř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 aproxi-maci po částech polynomiální funkci - tzv. splajn. Splajny byly studoványuž začátkem 20. století. Označení „splineÿ použil poprvé v roce 1946 Shoen-berg. Původně byla tímto názvem označována speciální pružná dřevěná nebokovová 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. Funkcis(x) : 〈a, b〉 → R, pro kterou platí

s ∈ Cm−1〈a, b〉,s |<xi−1,xi>∈ Pn, n ≤ m, i = 1, 2, . . . , n,

budeme nazývat splajnem m-tého stupně.Množinu všech těchto splajnů označíme Sn

m.

Nyní si ukážeme způsob, jak je možné zkonstruovat kubický splajn. V uzlo-vých bodech x0, x1, . . . , xn ∈ 〈a, b〉 jsou dány funkční hodnoty f0, f1, . . . , fn

interpolované funkce. Požadujeme, aby kubický splajn s(x), který chcemenajít, splňovals ∈ C2〈a, b〉,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

s′′i (x) = 2ci + 6di(x− xi). (5.14)

Předpokládejme na chvíli, že známe s′′(xi) = Mi. (Mi jsou tzv. momenty.)Funkci s′′i (x) pak můžeme na intervalu < xi, xi+1 > nahradit lineárním poly-nomem

s′′i (x) = Mixi+1 − x

hi

+ Mi+1x− xi

hi

. (5.15)

Když tento vztah integrujme, dostaneme

s′i(x) = −Mi(xi+1 − x)2

2hi

+ Mi+1(x− xi)2

2hi

+ Ai, (5.16)

si(x) = Mi(xi+1 − x)3

6hi

+ Mi+1(x− xi)3

6hi

+ Ai(x− xi) + Bi. (5.17)

109

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 = fi − Mih2i

6. (5.18)

Pro x = xi+1 z (5.17) je

Ai =fi+1

hi

−Mi+1h3

i

6h2i

− fi

hi

+ Mih2

i

6hi

= (5.19)

=fi+1 − fi

hi

− (Mi+1 −Mi)hi

6.

To znamená, že pro x ∈ 〈xi, xi+1〉, i = 0, . . . , n− 1, je s(x) = si(x) a

si(x) = Mi(xi+1 − x)3

6hi

+ Mi+1(x− xi)3

6hi

+ Ai(x− xi) + Bi.

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 = f ′′0 , Mn = f ′′n .

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 〈a, b〉,spojitou první derivaci tj.

s′i−1(x+i ) = s′i(x

−i ), i = 1, . . . , n− 1.

Pak podle (5.16), (5.19)

fi − fi−1

hi−1+

hi−1

3Mi +

hi−1

6Mi−1 =

fi+1 − fi

hi

− hi

3Mi − hi

6Mi+1,

hi−1

hi−1 + hi

Mi−1 + 2Mi +hi

hi−1 + hi

Mi+1 =6

hi−1 + hi

(fi−1 − fi

hi

− fi − fi−1

hi−1).

Označíme-li

αi =hi−1

hi−1 + hi

,

110

βi =hi

hi−1 + hi

,

γi =6

hi−1 + hi

(fi−1 − fi

hi

− fi − fi−1

hi−1),

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

2 β1 0 0 . . . 0 0α2 2 β2 0 . . . 0 0. . . . . . . . . . . . . . . . . . . . .

0 0 0 0 αn−2 2 βn−2

0 0 0 0 0 αn−1 2

M1

M2

. . .Mn−2

Mn−1

=

γ1 − α1f′′0

γ2

. . .γn−2

γn−1 − βn−1f′′n

.

Algoritmus (Konstrukce kubického splajnu)

Vstup: n, x0, . . . , xn, f0, . . . , fn, f ′′0 , f ′′npro i = 0, . . . , n− 1

hi = xi+1 − xi

pro i = 1, . . . , n− 1αi = hi−1

hi−1+hi

βi = hi

hi−1+hi= 1− αi

γi = 6hi−1+hi

(fi−1−fi

hi− fi−fi−1

hi−1)

M0 = f ′′0 , Mn = f ′′nM1, . . . , Mn−1 — řešení soustavy (5.20)

pro i = 1, . . . , n− 1 spočteme

Bi = fi − Mih2i

6

Ai = fi+1−fi

hi− (Mi+1−Mi)hi

6

Výstup: si(x) = Mi(xi+1−x)3

6hi+ Mi+1

(x−xi)3

6hi+ Ai(x− xi) + Bi

111

Příklad Funkce je zadána tabulkou

xi −1 0 1 3fi 3 3 1 2

Najděte přirozený kubický splajn (tj. takový, pro který jsou momenty M0 =Mn = 0.)

i 0 1 2 3xi −1 0 1 3fi 3 3 1 2hi 1 1 2 −αi − 1

213 −

βi − 12

23 −

γi − 0 6 −

M0 = M3 = 0

2M1 +12M2 = 0

13M1 + 2M2 = 6

Odtud:

M1 = −1823

, M2 =7223

.

Dále

A0 =323

, A1 = −6123

, A2 =7146

,

B0 = 3, B1 =7223

, B2 = −2523

.

Dostáváme tedy

s0(x) =7223− 3

23(x + 1)3 +

323

x pro x ∈ 〈−1, 0〉

s1(x) = − 323

(1− x)3 +1223

x3 − 6123

x +7223

pro x ∈ 〈0, 1〉

112

s2(x) =623

(3− x)3 − 12146

− 7146

x pro x ∈ 〈1, 3〉.

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

–0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4

–1 1 2 3x

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

uk(x) = (x− x0)k, k = 0, . . . m, (5.21)

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

Na jednotlivých subintervalech intervalu 〈a, b〉 potom máme

s(x) =m∑

k=0

αkuk(x) +i−1∑

k=1

βkvk(x), x ∈ 〈x0, xi〉, i = 1, . . . , n. (5.22)

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

3Zde značíme

(x− xk)+ =

{x− xk, když x− xk ≥ 0,0 jinak.

113

α0, . . . , αm, β0, . . . , βn−1. Přidáme proto ještě dalších m − 1 podmínek proderivace 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é koeficientypak určíme ze soustavy

m∑

k=0

αkuk(xi) +n−1∑

k=1

βkvk(xi) = fi, i = 0, . . . , n, (5.23)

m∑

k=0

αku(i)k (a) +

n−1∑

k=1

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

m∑

k=0

αku(i)k (b) +

n−1∑

k=1

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

Tato soustava je sice jednoznačně řešitelná, ale kvůli zvoleným bázovým funk-cím (5.21), které mají globální nosič, je špatně podmíněná. Z tohoto důvodubývá zvykem za bázové funkce volit B-splajny

B0(x) =

{1 pro |x| ≤ 0, 5,0 pro |x| > 0, 5.

Bm+1(x) =∫ x+ 1

2

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

B2(x) =

2− (|x| − 12)2 − (|x|+ 1

2)3 pro |x| ≤ 12 ,

(|x| − 32)2 pro 1

2 < |x| < 32 ,

0 pro |x| ≥ 32 .

Když m = 3, dostaneme

B3(x) =16

(2− |x|)2 − 4(1− |x|)3 pro |x| ≤ 1(2− |x|)3 pro 1 < |x| < 20 pro |x| ≥ 2

114

Tyto funkce už mají lokální nosič. Dá se ukázat, že když uzly jsou v intervalu〈a, b〉 rozloženy ekvidistantně, h = b−a

n, n ≥ 2, m = 2l − 1, pak za bázi

prostoru Snm můžeme vzít B-splajny

Bm(x− xk

h), x ∈ 〈a, b〉, k = −l + 1, . . . , n + l − 1.

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 uve-dené interpolační polynomy právě nejvhodnější. Místo prostoru polynomůPn budeme v takovém případě raději pracovat s prostorem trigonometric-kých polynomů. Poprvé tento typ interpolace použili Clairout v roce 1759 aLagrange v roce 1782.

Definice 5.2 Trigonometrickým polynomem rozumíme funkci

g(x) =n∑

k=0

ak cos kx +n∑

k=1

bk sin kx,

kde n ∈ N, an, bn ∈ C, |an| + |bn| > 0. Množinu všech takových trigonomet-rický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 za-dána tabulkou

xi 0 2π3

4π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,

a0 − 12a1 +

√3

2b1 = 0,

a0 − 12a1 −

√3

2b1 = −1.

Odtud a0 = 0, a1 = 1, b1 =√

33 , g(x) = cos x +

√3

3 sin x.

Předpokládejme, že funkce, kterou aproximujeme, je periodická s periodou2π, a že interpolujeme na intervalu 〈0, 2π〉. Bez újmy na obecnosti budemenadále předpokládat, že uzlové body jsou v intervalu 〈0, 2π〉 rozloženy ekvi-distantně.

Věta 5.5 Když jsou dány ekvidistantní uzlové body x0, . . . , x2n ∈ 〈0, 2π) av nich jsou zadány funkční hodnoty f0, . . . , f2n ∈ R, pak existuje právě jedeninterpolační polynom g ∈ Tn tak, že

g(x) =a0

2+

n∑

k=1

[ak cos kx + bk sin kx], (5.24)

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

g(2πj

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

Jeho koeficienty jsou určeny vztahy

ak =2

2n + 1

2n∑j=0

fj cos2πj

2n + 1k, k = 0, . . . , n, (5.25)

bk =2

2n + 1

2n∑j=0

fj sin2πj

2n + 1k, k = 1, . . . , n. (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 figurujefunkce 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 ∈〈0, 2π) a hodnoty f0, . . . , f2n−1 ∈ R, pak existuje právě jeden interpolačnípolynom g ∈ Tn,

g(x) =a0

2+

n−1∑

k=1

[ak cos kx + bk sin kx] +an

2cos nx, (5.27)

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

g(πj

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

Jeho koeficienty jsou určeny vztahy

ak =1n

2n−1∑j=0

fj cosπj

nk, k = 0, . . . , n, (5.28)

bk =1n

2n−1∑j=0

fj sinπj

nk, k = 1, . . . , n− 1. (5.29)

Příklad Určete trigonometrický interpolační polynom, když je funkce za-dána tabulkou

xi 0 2π3

4π3

fi 1 0 −1

Řešení: Počet ekvidistantních uzlů je lichý. Podle věty 5.5 obdržíme trigono-metrický polynom g(x) = a0

2 + a1 cos x + b1 sin x, kde

a0 =23

(1 · cos 0 + 0 · cos 0− 1 · cos 0) = 0,

a1 =23

(1 · cos 0 + 0 · cos2π

3− 1 · cos

4π

3) = 1,

b1 =23

(1 · sin 0 + 0 · sin 2π

3− 1 · sin 4π

3) =

√3

3.

117

Odtud g(x) = cos x +√

33 sin x.

Poznámka Uvedená aproximace periodických funkcí trigonometrickými po-lynomy 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í trigonomet-rického polynomu jako funkce komplexní proměnné. Díky Eulerovým vzor-cům je

g(x) =a0

2+

n∑

k=1

[ak cos kx + bk sin kx] =a0

2+

n∑

k=1

[ak12

(eikx + e−ikx)+

+ bk12i

(eikx − e−ikx)] =a0

2+

12

n∑

k=1

[(ak + ibk)e−ikx + (ak − ibk)eikx].

Když označíme c0 = a02 , ck = 1

2(ak − ibk) a c−k = 12(ak + ibk), k = 1, . . . n,

lze trigonometrický polynom psát ve tvaru

g(x) =n∑

k=−n

ckeikx. (5.30)

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

ck =1

2n + 1

2n∑j=0

fje−ikxj , k = −n, . . . , n.

Poznámka Když položíme w = e−ixj , pak ck = P (wk). Vidíme, že ko-eficienty ck jsou rovny hodnotám polynomů (2n + 1)-ního stupně v bodechwk, 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 aTukey. Předpokládejme, že k dispozici máme N = 2n ekvidistantních uzlů.Označme

w = e−i 2πN , Fj =

fj

N.

Pak

ck =N∑

j=0

Fjwkj.

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

ck =

N2 −1∑j1=0

F2j1(w2)kj1 + wk

N2 −1∑j1=0

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

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

(w2)kj1 = (w2)r N2 j1(w2)k1j1 = (wN)rj1(w2)k1j1 = (w2)k1j1 .

Vidíme, že

ck =

N2 −1∑j1=0

F2j1(w2)k1j1 + wk

N2 −1∑j1=0

F2j1+1(w2)k1j1 = Pk1(wk) + wkPk1(wk).

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

Hermitova interpolace

Pokud kromě funkčních hodnot známe i hodnoty derivací funkce v uzlovýchbodech, je možné konstruovat interpolační polynom stupně 2n + 1.

Věta 5.7 (Hermitova interplace) Když jsou dány body x0, . . . , xn ∈ 〈a, b〉a hodnoty f0, . . . , fn, f ′0, . . . , f

′n ∈ R, pak existuje právě jeden interpolační

119

polynom p ∈ P2n+1, který vyhovuje interpolačním podmínkám p(xi) = fi,p′(xi) = f ′i , i = 0, . . . , n. Tento tzv. Hermitův polynom je definován vztahem

H(x) =n∑

i=0

[fiH0i + f ′iH

1i ],

kde H0i (x) = [1− 2l′i(xi)(x− xi)][li(x)]2, H1

i (x) = (x− xi)[li(x)]2 a li(x) jsouLagrangeovy polynomy (5.5).

Důkaz. Protože H0i , H1

i ∈ P2n+1, i = 0, . . . , n, je Hermitův polynom H ∈P2n+1. Protože

H0i (xj) = H1

i′(xj) = δi,j,

H0i′(xj) = H1

i (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í polynomyH1, H2 ∈ P2n+1, H1 6= H2. Pak polynom H = H1 − H2 splňuje v uzlovýchbodech vztahy

H(xj) = H ′(xj) = 0, j = 0, . . . , n.

H je polynom stupně 2n + 1 a má n + 1 dvojnásobných kořenů. To aleznamená, ž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 au-tomobilové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é

4Bernsteinovy polynomy hrály důležitou roli v Bernsteinově důkazu Weierstrassovyvě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é reprezentacikř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 je-jich 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í odhadnouttvar křivky.V dalším výkladu se pro jednoduchost omezíme pouze na interval 〈0, 1〉. Tovšak nebude újmou na obecnosti, protože k obecnému intervalu 〈a, b〉 lzedospět prostřednictvím transformace t = x−a

b−a.

Definice 5.3 Když pro t ∈ 〈0, 1〉 označíme

Bni (t) =

(n

i

)ti(1− t)n−i, i = 0, . . . , n, (5.31)

Bernsteinův polynom n-tého stupně, pak

p(t) =n∑

i=0

biBni (t) (5.32)

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. Je-jich konvexní obal

conv {b1, b2, . . . , bn} =n∑

i=0

αibi, αi ≥ 0,

n∑i=1

αi = 1,

je tzv. Bézierův polygon.

Příklad Nakreslete Bézierovu křivku, když jsou dány řídicí bodya) 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 ∈ 〈0, 1〉. Parame-trické rovnice křivky:

x(t) = (1− t)2 + 8t(1− t) + 5t2,

y(t) = 3(1− t)2 + 8t(1− t) + 2t2, t ∈ 〈0, 1〉.

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 ∈ 〈0, 1〉.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 ∈ 〈0, 1〉.

Poznámka Důležité vlastnosti Bézierových polynomů:1. V krajních bodech intervalu 〈a, b〉 je hodnota Bézierova polynomu rovnahodnotě odpovídajícího kontrolního bodu.2. Tečny Bézierova polynomu a hrany Bézierova polygonu v krajních bodechsplývají. Viz následující obrázek (a).3. Jestliže Bézierův polygon je konvexní (konkávní), je také konvexní (kon-ká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í seve 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 dvaBézierovy polynomy stýkají, pak splynou. Viz následující obrázek (b).

Obrázek: Vlastnosti Bézierových polynomu 2.stupně.

(a)b1

b0 b2

(b)b′1

b0 b2 = b′0 b′2

b1

122

Určení hodnot Bézierova polynomu p(t) pro t ∈ 〈0, 1〉 se provádí Casteljau-ový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∑

i=0

biB2i (t) = (1− t2)2b0 + 2t(1− t)b1 + t2b2 =

= (1− t) [(1− t)b0 + tb1]︸ ︷︷ ︸b10

+ [(1− t)b1 + tb2]︸ ︷︷ ︸b11

.

To znamená, že když zvolíme t0 ∈ 〈0, 1〉 pevně, můžeme určit bod p(t0), kterýleží na Bézierově křivce tak, že položíme

b10 = (1− t0)b0 + t0b1 a b1

1 = (1− t0)b1 + t0b2

a spočteme hodnotup(t0) = (1− t0)b1

0 + t0b11.

Obrázek: Konstrukce bodu na Bézierově křivce n = 2, t = 12 .

b1

b10 b2

0 b11

b0 b2

Obecně, když je dáno n+1 řídících bodů, položíme b00 = b0, b0

1 = b1, . . . , b0n =

bn a postupně pro k = 0, 1, . . . , n a j = 0, 1, . . . , n−k počítáme subpolynomy

bkj (t) = (1− t)bk−1

j (t) + tbk−1j+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−10

bn b1n−1

... bn−11 bn

0

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

i + tbk−1i+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ázelazadanými body, je pokusit se najít funkci, která je těmto hodnotám v jis-tém smyslu co nejblíž. Pokud velikost odchylky aproximující a dané funkcebudeme 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∑

j=0

cjϕj(x) (5.33)

a snažíme se najít konstanty c1, . . . , cm tak, aby číslo

S(c1, c2, . . . , cm) =n∑

i=0

[f(xi)− ϕ(xi)]2 =

n∑i=0

[f(xi)−m∑

j=0

cjϕj(xi)]2

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

∂S

∂ck

= 2n∑

i=0

[f(xi)−m∑

j=0

cjϕj(xi)]ϕk(xi) = 0, k = 1, . . . , m.

Dostáváme tak soustavu tzv. normálních rovnic

m∑j=0

cj

n∑i=0

ϕj(xi)ϕk(xi) =n∑

i=0

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

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∑j=0

(ϕj, ϕk)cj = (f, ϕk), k = 0, . . . , n. (5.35)

Poznámka Pokud zvolíme ϕj(xi) = xji , j = 1, . . . ,m, znamená to, že funkci

f aproximujeme polynomem m-tého stupně. Pak funkce ϕ(x) =∑m

j=0 cjxji a

soustava normálních rovnic (5.34) přejde na tvar

m∑j=0

(n∑

i=0

xj+ki )cj =

n∑i=0

f(xi)xki , k = 0, . . . , n.

5Za základní funkce lze zvolit např. {1, x, x2, . . . , xm} nebo{1, cos x, sin x, cos 2x, sin 2x, . . . , sin mx}

125

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

ϕj,k =∑n

i=0 ϕj(xi)ϕk(xi)

ϕk,j = ϕj,k

fk =∑n

i=0 fiϕj(xi).

c0, c1, . . . , cm — řešení soustavy

ϕ0,0 ϕ0,1 . . . ϕ0,m

ϕ1,0 ϕ1,1 . . . ϕ1,m

. . . . . . . . . . . .ϕm,0 ϕm,1 . . . ϕm,m

c0

c1

. . .cm

=

f0

f1

. . .fm

.

Výstup: ϕ(x) =∑

j cjϕj(x).

Příklad Funkce je zadána tabulkou

xi 1 2 3 5fi 3 3 1 2

Pomocí metody nejmenších čtverců aproximujte tuto funkci kvadratickýmpolynomem.

Řešení: Funkci aproximujeme polynomem ϕ(x) = c0 + c1x + c2x2, jehož ko-

eficienty 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∑

i=0

ϕ0(xi)ϕ0(xi) =3∑

i=0

1 = 4,

(ϕ0, ϕ1) = (ϕ1, ϕ0) =3∑

i=0

ϕ1(xi)ϕ0(xi) =3∑

i=0

xi = 11,

(ϕ0, ϕ2) = (ϕ2, ϕ0) =3∑

i=0

ϕ2(xi)ϕ0(xi) =3∑

i=0

x2i = 39,

(ϕ1, ϕ1) =3∑

i=0

ϕ1(xi)ϕ1(xi) =3∑

i=0

x2i = 39,

(ϕ1, ϕ2) = (ϕ2, ϕ1) =3∑

i=0

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

i=0

x3i = 161,

(ϕ2, ϕ2) =3∑

i=0

ϕ2(xi)ϕ2(xi) =3∑

i=0

x4i = 723,

(f, ϕ0) =3∑

i=0

f(xi)ϕ0(xi) = 9,

(f, ϕ1) =3∑

i=0

f(xi)ϕ1(xi) = 22,

(f, ϕ1) =3∑

i=0

f(xi)ϕ1(xi) = 74.

Řešíme tedy soustavu4c0 + 11c1 + 39c2 = 9,

11c0 + 39c1 + 161c2 = 22,

39c0 + 161c1 + 723c2 = 74.

Odtud c2 = 14 , c1 = −37

20 a c0 = 4910 . Metodou nejmenších čtverců jsme získali

aproximaci zadané funkce ve tvaru

ϕ(x) =14x2 − 37

20x +

4910

.

127

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

0

1

2

3

4

5

6

y

1 2 3 4 5 6x

ϕ(x) = 14x

2 − 3720x + 49

10 .

128

Kapitola 6

Numerická kvadraturaa derivace

Numerickou metodu pro integrování neboli numerickou kvadraturu použi-jeme v případě, když integrál nelze určit analyticky prostřednictvím elemen-tárních funkcí nebo když integrovaná funkce je zadaná tabulkou. Principnumerické kvadratury spočívá v tom, že funkci f(x) nahradíme jinou (jed-nodušší) funkcí ϕ a tu pak analyticky integrujeme. Pokud je ϕ dobrou apro-ximací funkce f na nějakém intervalu 〈a, b〉, je také

∫ b

aϕ(x) dx dobrou apro-

ximací∫ b

af(x) dx.

Když x0, x1, . . . , xn ∈ 〈a, b〉 a známe hodnoty integrované funkce v uzlovýchbodech, můžeme ji nahradit třeba Lagrangeovým polynomem

f(x) ≈ Ln(x) =∑n

i=0 f(xi)li(x).

Integrací obdržíme∫ b

af(x) dx ≈ ∫ b

aLn(x) dx =

∫ b

a

∑ni=0 f(xi)li(x) dx =

∑ni=0 f(xi)(

∫ b

ali(x) dx).

Když označíme wi =∫ b

ali(x) dx, obdržíme kvadraturní formuli∫ b

a

f(x) dx ≈n∑

i=0

wif(xi). (6.1)

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

129

chybu aproximace platit

|∫ b

a

f(x) dx−∫ b

a

Ln(x) dx| ≤∫ b

a

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

Vidíme, že chyba numerické kvadratury závisí na tom, jak velká je chybaaproximace funkce.

Poznámka V dalším textu budeme používat pro zadaný integrál označení

I =∫ b

a

f(x) dx (6.3)

a pro kvadraturní formuli

Qn =n∑

i=0

wifi. (6.4)

Poznamenejme ještě, že integrál i kvadratura představují lineární operátory

I, Qn : C0〈a, b〉 → R.

Definice 6.1 Pokud uzlové body volíme tak, že a = x0 a b = xn, hovořímeo 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ě roz-ložených uzlech. Označme tyto uzly

x0, x0 + h, . . . , x0 + nh ∈ 〈a, b〉, h =b− a

n.

130

Příklad Pro n = 1 zkonstruujte Newtonovu-Cotesovu kvadraturní formulia odhadněte chybu kvadraturního vzorce.

Řešení: Položíme h = b−a2 , x0 = a, x1 = b. Integrovanou funkci f aproximu-

jeme lineárním polynomem L1(x) = x−ba−b

f(a) + x−ab−a

f(b). Pak

∫ b

a

f(x) dx ≈∫ b

a

L1(x) dx =b− a

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

Nyní odhadneme chybu. Lze psát∫ b

a

f(x) dx− b− a

2(f(a) + f(b)) =

∫ b

a

(x− a)(x− b)f(x)− L1(x)(x− a)(x− b)

dx.

Pokud je funkce f ∈ C2(a, b), pak z věty o střední hodnotě máme: existujetakové y ∈ 〈a, b〉, že

∫ b

a

f(x) dx− b− a

2(f(a) + f(b)) =

f(y)− L1(y)(y − a)(y − b)

∫ b

a

(x− a)(x− b) dx.

Odtud a z odhadu chyby (5.8)∫ b

a

f(x) dx− b− a

2(f(a) + f(b)) =

f ′′(ξ)2

(−b− a

6). (6.5)

Poznámka Kvadraturní formule

Q1 =b− a

2(f(a) + f(b)), (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 naintervalu 〈a, b〉 kvadraturní formuli prvního a druhého řádu.

Řešení: Podle vztahu (6.4) je kvadraturní formule 1. řádu Q1 = w0f0 + w1f1.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))

w0 + w1 =∫ b

a

1 dx = b− a,

−w0a + w1b =∫ b

a

x dx =b2 − a2

2.

Odtud w0 = w1 = b−a2 a

∫ b

af(x) dx ≈ b−a

2 (f0 + f1) = Q1(x).Ukázali jsme, že lichoběžníkové pravidlo představuje kvadraturní formuli 1.řádu.

Podle (6.4) kvadraturní formule 2. řádu Q2 = w0f0 + w1f1 + w2f2. Zvolmex0 = a, x1 = a+b

2 , x2 = b. Z požadavku, aby kvadraturní formule byla přesnápro polynomy p0(x) = 1, p1(x) = x, p2(x) = x2, získáme soustavu

w0 + w1 + w2 =∫ b

a

1 dx = b− a,

w0a + w1a + b

2+ w2b =

∫ b

a

x dx =b2 − a2

2,

w0a2 + w1

(a + b)2

4+ w2b

2 =∫ b

a

x2 dx =b3 − a3

3.

Odtud a0 = a2 = b−a3 , a1 = b−a

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

Q2(x) =b− a

3(f(a) + 4f(

a + b

2) + f(b)) = Q2(x), (6.7)

tzv. Simpsonovo pravidlo.

Poznámka Obdobně lze odvodit hodnoty wi i pro ostatní Newtonovy-Cotesovy vzorce. Když interval, přes který integrujeme, je 〈−1, 1〉, pak prokoeficienty v Newtonových-Cotesových vzorcích dostaneme

n wi

1 Lichoběžníkové pravidlo 12

12

2 Simpsonovo pravidlo 13

43

13

3 Newtonovo 3/8 pravidlo 38

98

98

38

4 Milnovo pravidlo 1445

6445

2445

6445

1445

132

Složené kvadraturní vzorce

Položme si nyní otázku, zda lze nějakým způsobem dosáhnout zvýšení přes-nosti kvadraturního vzorce. Přidáním uzlových bodu 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 kvadratur-ních vzorcu. Postupujeme tak, že interval 〈a, b〉 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

∫ b

a

f(x) dx =m−1∑j=0

(∫ aj+1

aj

f(x) dx) ≈m−1∑j=0

nj∑i=0

wijfi, (6.8)

kde

wij =∫ aj+1

aj

lij(x) dx.

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 〈xj, xj+1〉 nahradíme funkcif(x) hodnotou f(xj + h

2 ), dospějeme ke kvadraturní formuli

R(f, h) =n−1∑j=0

∫ xj+1

xj

f(x) dx ≈ h

n−1∑j=0

f(xj +h

2). (6.9)

Chyba aproximace je

|∫ b

a

f(x) dx−∫ b

a

R(f, h) dx| ≤ (b− a)h2

24f ′′(ξ), ξ ∈ 〈a, b〉. (6.10)

Algoritmus: (Obdélníkové pravidlo)

Vstup: a, b, f(x), nh = b−a

n

pro j = 0, 1, . . . , nxj := a + jh

R(f, h) = h∑n−1

j=0 f(xj + h2 )

Výstup: R(f, h)

133

b) Lichoběžníkové pravidlo. Pokud na intervalu 〈xj, xj+1〉 aproximujeme funkcif(x) aritmetickým pruměrem 1

2(f(xj) + f(xj+1)), pak získáme kvadraturníformuli

T (f, h) =n−1∑j=1

∫ xj+1

xj

f(x) dx ≈ h

n−1∑j=0

12

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

=h

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

Chyba aproximace je zde

|∫ b

a

f(x) dx−∫ b

a

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

12f ′′(ξ), ξ ∈ 〈a, b〉. (6.12)

Algoritmus: (Lichoběžníkové pravidlo)

Vstup: a, b, f(x), nh = b−a

n

pro j = 0, 1, . . . , nxj := a + jh

T (f, h) = h2

∑n−1j=0 (f(xj) + f(xj+1))

Výstup: T (f, h)

c) Simpsonovo pravidlo. Další kvadraturní vzorec obdržíme, když na inter-valu 〈xj, xj+2〉 funkci f(x) nahradíme kvadratickým polynomem

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

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

(x− xj)(x− 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∑j=0

∫ xj+1

xj−1

f(x) dx ≈ h

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

(6.13)

134

Chyba aproximace

|∫ b

a

f(x) dx−∫ b

a

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

180f (4)(ξ), ξ ∈ 〈a, b〉. (6.14)

Algoritmus: (Simpsonovo pravidlo)

Vstup: a, b, f(x), n

h = b−an

pro j = 0, 1, . . . , nxj := a + jh

S(f, h) = h3

n2−1∑j=0

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

Výstup: S(f, h)

Příklad Simpsonovým pravidlem spočítejte∫ 2

1ex

xdx.

Řešení: Pro n = 6 obdržíme∫ 2

1

ex

xdx

.=

118

e+421

e76 +

112

e43 +

427

e32 +

115

e53 +

433

e116 +

136

e2 = 3, 059142296,

|∫ 2

1

ex

xdx− S(f, h)| ≤ 1

18016

15 = 0, 0052083.

Ve vzorcích (6.10), (6.12) a (6.14) se vyskytují hodnoty 2. a 4. derivace funkcef. Když však je integrovaná funkce zadaná tabulkou, nemůžeme hodnotyjejích derivací stanovit přesně. Můžeme jen získat jejich hrubý odhad. Toznamená, že také dostáváme jen hrubý odhad chyby. Z tohoto důvodu sek 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ů pro-vedeme odhad prostřednicvím vztahů

|∫ b

a

f(x) dx−∫ b

a

R2n(f, h) dx| ≤ 13|Rn(f, h)−R2n(f, h)|,

135

|∫ b

a

f(x) dx−∫ b

a

T2n(f, h) dx| ≤ 13|Tn(f, h)− T2n(f, h)|,

|∫ b

a

f(x) dx−∫ b

a

S2n(f, h) dx| ≤ 115|Sn(f, h)− S2n(f, h)|.

Příklad Lichoběžníkovým pravidlem spočítejte∫ 2

1ex

xdx. (Volte n = 6.) Od-

had chyby kvadratury proveďte pomocí vztahu (6.12) i Rungeovou metodoudvojnásobného kroku.

Řešení: Pro n = 6 je∫ 2

1

ex

xdx ≈ T6 =

112

e+124

e2+17

e76 +

18

e43 +

19

e32 +

110

e53 +

111

e116 = 3, 063385879,

|∫ 2

1

ex

xdx− T6(f, h)| ≤ h3

12M2 = (

16

)3 112

2, 718 = 0, 10487 · 10−2.

Pro n = 3 dostaneme∫ 2

1

ex

xdx ≈ T3 =

16

e +14

e43 +

15

e53 +

112

e2 = 3, 076116630,

|∫ 2

1

ex

xdx− T6(f, h)| ≤ 1

3|T3(f, h)− T6(f, h)| = 0, 4243583667 · 10−2.

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ásle-dující skutečnost: Pokud funkce f má derivace až do (m + 1)-ního řádua B2 = 1

6 , B4 = − 130 , B6 = 1

42 , . . . , Bm jsou Bernoulliova čísla, lze chybulichoběžníkového pravidla psát ve tvaru

T (f, h)−∫ b

a

f(x) dx =m∑

s=0

Bs

2s!h2s[f (2s−1)(b)− f (2s−1)(a)]

(tzv. Eulerův-Maclaurinův vzorec).Chybu složené lichoběžníkové formule je tedy možné rozvinout do řady sudýchmocnin integračního kroku h

T (f, h)−∫ b

a

f(x) dx = a1h2 + a2h

4 + · · ·+ amh2m + O(h2m+2). (6.15)

136

To znamená (viz poznámku na str. 107), že je možné aplikovat Richardsonovuextrapolaci.Odvození nových formulí: Když funkce f ∈ C2〈a, b〉, h = b−a

n, lichoběžníkové

pravidlo pro krok h označíme T 00 = T (h, f) a pro krok h

2 označíme T 01 =

T (h2 , f), pak podle vztahu (6.15) platí

∫ b

a

f(x) dx = T 00 + a1h

2 + O(h4),

∫ b

a

f(x) dx = T 01 + a1

h2

4+ O(h4).

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

∫ b

a

f(x) dx =13

(4T 01 − T 0

0 ) + O(h4).

Dále tuto lineární kombinaci složených lichoběžníkových pravidel budemeznačit

T 11 =

13

(4T 01 − T 0

0 ).

Obdobně pro f ∈ C4〈a, b〉 obdržíme

∫ b

a

f(x) dx = T 11 + a2h

4 + O(h6),

∫ b

a

f(x) dx = T 12 + b2

h4

42+ O(h6).

Opět eliminujeme a2 a dostaneme novou kvadraturní formuli

T 22 =

115

(16T 12 − T 1

1 ),

která je ale už přesnější než předchozí formule T 11 , protože aproximuje hod-

notu 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

Tm+1k =

4mTmk+1 − Tm−1

k

4m − 1.

137

Pokud f ∈ C2m〈a, b〉, platí pro odhad chyby

|∫ b

a

f(x) dx− Tmk | ≤ c(m)‖f (2m)‖∞(

h

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

Protože Rombergova kvadratura představuje v podstatě aplikaci Richard-sonovy extrapolace na lichoběžníkové pravidlo, lze výpočet hodnoty Tm

k přiručním počítání provést v tabulce

k 2k m = 0 m = 1 m = 20 1 T 0

0

1 2 T 01 T 1

1

2 4 T 02 T 1

2 T 22

Hodnoty T 11 , T 2

2 , . . . na diagonále tabulky představují aproximace hodnotyintegrálu s chybami řádu O(h4), O(h6), . . . Podle nich hned zjistíme, kdy jsmedosá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ší řádkytabulky. 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

T 0i = T (f, h

2i )

pro j = i, i + 1, . . . , m

T ji = T j−1

i +T j−1

i −T j−1i−1

4j−1

Výstup: I ≈ Tmm

Příklad Spočítejte∫ 1

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

138

Řešení:

i T 0i T 1

i T 2i T 3

i T 4i T 5

i

0 0.1732871 0.248829 0.2740102 0.256458 0.272334 0.2722223 0.270769 0.272206 0.272197 0.2721974 0.271841 0.272199 0.272198 0.272198 0.2721985 0, 272109 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ž inte-grovaná funkce je dost hladká. Není to však metoda univerzální. Pro integrályz 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

Uvažujme nyní integrál v obecnějším tvaru∫ b

aω(x)f(x) dx, kde ω je tzv. vá-

hová funkce 1. Opět budeme konstruovat kvadraturní formuli ve formě součtu∑i wif(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 va-hou ω na intervalu 〈a, b〉. 2 Výhodou Gaussových kvadratur je, že při stejnémpočtu uzlů dávají přesnější výsledky než Newtonovy-Cotesovy vzorce.

1Tato funkce je na uvažovaném intervalu spojitá a kladná. Např. ω(x) = 1, ω(x) =√1− x2 nebo ω(x) = 1√

1−x2 .2Viz definice (1.13). Příslušný skalární součin polynomů pi(x), pj(x) má tvar

(pi, pj) =∫ b

a

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

139

Pro ω(x) = 1 a 〈a, b〉 = 〈−1, 1〉 dostáváme tzv. Gaussův-Legendrův kvadra-turní vzorec ∫ 1

−1f(x) dx =

n∑i=0

wif(xi). (6.16)

Uzlové body xi zde volíme jako kořeny Legendrových polynomů

Pn(x) =1

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

dxn

a koeficienty wi jsou definovány vztahy

wi =−2

(n + 2)Pn+2(xi)P ′n+1(xi)

.

Hodnoty uzlových bodů a koeficientů v Gaussově-Legendrově vzorci nemu-sí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

12 0 8

9

±√

35

59

3 ±0, 861136311 0, 347854845±0, 339981043 0, 652145154

Poznámka Gaussovy-Legendrovy kvadraturní vzorce pro

n = 0 :∫ 1−1 f(x) dx = f(0) + 1

3f′′(ξ), ξ ∈ (−1, 1).

n = 1 :∫ 1−1 f(x) dx = f(− 1√

3) + f( 1√

3) + 1

135f(4)(ξ), ξ ∈ (−1, 1).

n = 2 :∫ 1−1 f(x) dx = 5

9f(−√

35) + 8

9f(0) + 59f(

√35) + 1

15750f(6)(ξ),

ξ ∈ (−1, 1).

Poznámka Když pro x ∈ 〈−1, 1〉 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

|∫ 1

−1f(x) dx−

k∑i=0

wif(xi)| ≤ Mdn+1.

Poznámka Jestliže chceme integrovat přes interval 〈a, b〉 a stanovit hodnotuintegrálu

∫ b

af(x) dx, převedeme jej substitucí x = 1

2(t(b− a) + a + b) (dx =12(b−a) dt, t = 2x−a−b

b−a), na integrál b−a

2

∫ 1−1 f(1

2((b−a)t+a+ b)) dt. Pro jehovýpočet lze použít buď jednoduchý Gaussův-Legendrův vzorec

∫ b

a

f(x) dx =b− a

2

∫ 1

−1f(

12

(b− a)t + a + b)) dt ≈

≈ b− a

2

n∑i=0

wif(12

(xi(b− a) + a + b)),

kde chyba

∫ b

a

f(x) dx−∑

wif(xi) =f (2n)(ξ)

(2n + 2)!

∫ 1

−1qn+1(x)2 dx.

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

∫ b

a

f(x) dx =m−1∑j=0

∫ aj+1

aj

f(x) dx ≈ h

2

m−1∑j=0

(n∑

i=0

wif(12xih + a + h(j +

12

)),

ve kterém h = b−am

, xi = a + hi. Chyba tohoto kvadraturního vzorce je řáduO(h2m).

Algoritmus: (Složený Gaussův-Legendrův kvadraturní vzorec)

Vstup: a, b, f(x), m, n, w0, . . . , wn

h = b−am

pro j = 0, 1, . . . , nxj = a + hj.

G(f, h) = h2

∑mj=0

∑ni=0 wif(1

2xih + a + h(j + 12)).

Výstup: I ≈ G(f, h)

141

Příklad Pomocí Gaussova-Legendrova kvadraturního vzorce stanovte hod-notu integrálu

∫ 21

ex

xdx. (Volte m = 2, n = 1.)

Řešení: Integrační interval rozdělíme na dva subintervaly a na každém z nichintegrál nahradíme Gaussovým-Legendrovým kvadraturním vzorcem 1. řádu.

∫ 2

1

ex

xdx

.=

14

(f(5 + 1√

3

4)+f(

5− 1√3

4)+f(

7 + 1√3

4)+f(

7− 1√3

4)) = 3, 059035425.

Numerická derivace

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

f(x + h) = f(x) + hf ′(x) +h2

2!f ′′(x) +

h3

3!f ′′′(x) +

h4

4!f (4)(x) + O(h5),

f(x− h) = f(x)− hf ′(x) +h2

2!f ′′(x)− h3

3!f ′′′(x) +

h4

4!f (4)(x) + O(h5).

Odtudf(x + h)− f(x− h)

2h= f ′(x) +

h2

3!f ′′′(x) + O(h4).

Derivaci funkce f můžeme tedy přibližně spočítat pomocí vzorce

f ′(x) ≈ f(x + h)− f(x− h)2h

.

Obecný tvar vzorců pro numerickou derivaci můžeme odvodit také tak, žefunkci aproximujeme interpolačním polynomem a pak derivujeme. Nevýho-dou 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∑

i=0

li(t)fi +hn+1 qn+1(t)

(n + 1)!f (n+1)(ξ).

Zde t = x−xi

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

f (k)(x) =1hk

m∑i=0

l(k)i (t)fi +

hn+1−k qn+1−k(t)(n + 1)!

f (n+1)(ξ), (6.17)

Zatímco u numerické interpolace je chyba metody úměrná kroku h a zao-krouhlovací 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 me-tody 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, abyzaokrouhlovací 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íchpolynomů 1. nebo 2. stupně. Když známe hodnoty funkce f v ekvidistantníchuzlech x− h, x, x + h, pak je možné spočítat

f ′(x) =f(x + h)− f(x)

h− h

2f ′′(ξ),

f ′(x) =f(x)− f(x− h)

h+

h

2f ′′(ξ),

f ′(x) =f(x + h)− f(x− h)

2h− h2

6f ′′′(ξ),

f ′′(x) =f(x + h)− 2f(x) + f(x− h)

h2− h2f ′′′(ξ).

143

Příklad Stanovte přibližně f ′(0, 5), když f(x) = ex

xa h = 0, 1.

Řešení:xi 0.4 0.5 0.6fi 3, 729561745 3, 297442542 3, 036864667

Nejprve spočteme derivaci z hodnoty funkce v bodě x a jednom jeho soused-ním bodě

f ′(x) =f(x + h)− f(x)

h,

f ′(0, 5).=

f(0, 6)− f(0, 5)0.1

= −2, 60577875.

Když derivaci určíme z hodnot funkce ve dvou bodech, které sousedí s bodemx, máme

f ′(x) =f(x + h)− f(x− h)

2h,

f ′(0.5).=

f(0.6)− f(0.4)0.2

= −3, 46348539.

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

144

Literatura

[1] Buchanan J.I.,Turner P.R.: Numerical Methods and Analysis, New York1992

[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, 8aproximace funkcí, 95

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

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

diference vpřed, 105diference vzad, 105dobře podmíněná úloha, 12

extrapolace, 106

Gaussova eliminační metoda, 28Gaussova eliminační metoda s vý-

běrem hlavního prvku, 31Gaussova-Seidelova iterační metoda,

46Gaussovy kvadraturní vzorce, 139Gaussův-Legendruv kvadraturní vzo-

rec, 141Graefova metoda, 90

Hilbertův prostor, 19

charakteristický polynom, 15

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

počítači, 9

interpolace na ekvidistantních uz-lech, 105

iterační metody, 41

Jacobiova iterační metoda, 42Jacobiova metoda, 44

korektní úloha, 12

Lagrangeův interpolační polynom,97

Laguerrova metoda, 93Legenderovy polynomy, 21lichoběžníkové pravidlo, 131, 134lineární operátor, 22LU-rozklad matice, 35

metoda bisekce, 75metoda LU-rozkladu, 33metoda nejmenších čtverců, 124metoda prosté iterace, 76metoda Rayleighova podílu, 63metoda regula falsi, 80metoda sečen, 85mocninná metoda, 61

nejlepší aproximace vzhledem k mno-žině, 26

146

Nevillův algoritmus, 103Newtonova metoda, 82Newtonovy-Cotesovy kvadraturní vzorce,

130Newtonův interpolační polynom, 101norma, 18norma matice, 18norma operátoru, 22normovaný prostor, 18numerická derivace, 142numerická úloha, 8

obdélníkové pravidlo, 133odhad polohy vlastních čísel, 59ortogonální matice, 14ortogonální prvky, 20ortonormální prvky, 20ostře diagonálně dominantní matice,

14

platné číslice, 11podmíněnost soustav lineárních rov-

nic, 55podobné matice, 17pozitivně definitní matice, 14pásová matice, 15přesnost řešení, 42přímé metody řešení soustav line-

árních rovnic, 28přímý výpočet vlastních čísel, 65

Rayleighův podíl, 17Richardsonova extrapolace, 107Rombergova kvadratura, 136Rungeova metoda dvojnásobného kroku,

135

řešení nelineární rovnice, 73řešení rovnic Pn(x) = 0, 86

řešení soustav nelineárních rovnic,85

Simpsonovo pravidlo, 132, 134skalární součin, 19složené kvadraturní vzorce, 133splajny, 108spojitý lineární operátor, 22stabilita algoritmu, 13stanovení hodnoty determinantu, 39stanovení inverzní matice, 40superrelaxační metoda, 54, 55symetrická matice, 14

určení vlastních čísel metodou LU-rozkladu, 67

určení vlastních čísel metodou or-togonálních transformací, 69

úplný normovaný prostor, 21úplný problém vlastních čísel, 59,

65

vlastní číslo matice, 15

Weierstrassova věta, 96

147