2

´ Uvod do numerick´e matematiky. Pˇredn´aˇska pro posluchaˇce informatiky. Zimn´ı resp. Letn´ı semestr 2/2 Ivo Marek, Petr Mayer a Bohuslav Sekerka

1

´ Uvodn´ ı pozn´ amky.

Vymezen´ı problematiky vystihuje n´asleduj´ıc´ı charakteristika. Numerick´ a (v´ypoˇctov´ a) matematika je realizace matematick´ych model˚ u na v´ypoˇcetn´ı technice. Numerick´a matematika je tedy souˇc´ ast´ı matematiky. Je to vˇsak souˇc´ ast maj´ıc´ı oproti teoretick´e matematice jeden v´ yznamn´ y nedostatek. Ten spoˇc´ıv´ a v tom, ˇze ˇcistˇe teoretick´e matematick´e modely podstatn´ ym zp˚ usobem operuj´ı s realn´ ymi ˇc´ısly zat´ımco numerick´ a matematika je z´avisl´ a na v´ ypoˇcetn´ı technice a tud´ıˇz se mus´ı bez tohoto prostˇredku obej´ıt. Typick´ ym projevem t´eto skuteˇcnosti jsou zaokrouhlovac´ı chyby, coˇz vede ve sv´ ych d˚ usledc´ıch k numerick´e nestabilitˇe. Zdroj numerick´e nestability vˇsak nemus´ı b´ yt v´ yhradnˇe chyby ze zaokrouhlen´ı. Numerick´a nestabilita m˚ uˇze b´ yt zapˇriˇcinˇena t´eˇz vlastnostmi matematick´ ych model˚ u sam´ ych, na pˇr. realizace matematick´eho ˇspatnˇe podm´ınˇen´eho modelu na poˇc´ıtaˇci je pˇrirozen´ ym zdrojem numerick´e nestability. Zde pouˇzit´ y pojem ˇspatnˇe podm´ınˇen´y model je charakterizov´ an nespojitou z´avislost´ı v´ ystupn´ı informace na informac´ıch vstupn´ıch. ˇ Pˇ r´ıklad 1.1 Bud A=

”

a11 a21

a12 a22

•

,

kde ajk jsou re´aln´a ˇc´ısla. Symbolem h(A) oznaˇcme hodnost matice A. Snadno si uvˇedom´ıme, ˇze h nen´ı spojitou funkc´ı prvk˚ u ajk . Abychom si tuto okolnost uvˇedomili, staˇc´ı poloˇziti • ” a11 0 6 0, , a11 = A0 = 0 0 1

takˇze na jedn´e stranˇe h(A0 ) = 1, avˇsak, na stranˇe druh´e, pro libovoln´e re´aln´e b 6= 0 h(A(b)) = 2, kde A(b) =

”

a11 0

0 b

•

.

Pˇredchoz´ı pˇr´ıklad ukazuje jednak nespojitost hodnosti matice jakoˇzto funkce jej´ıch prvk˚ u jednak dalˇs´ı pozoruhodnou vlastnost hodnosti. Ta je obsahem n´asleduj´ıc´ı vˇety. Vˇ eta 1.1 Hodnost matice je zdola polospojitou funkc´ı jej´ıch prvk˚ u. Definice 1.1 Skal´arn´ı re´aln´ yv´ a zdola polospojitou v bodˇe s0 , a funkce f se naz´ jestliˇze mnoˇzina {s : f (s) > f (s0 )} je otevˇren´a, t. j. existuje δ > 0 takov´e, ˇze

f (s) − f (s0 ) ≥  pro s ∈ (s0 − δ, s0 + δ). Dalˇs´ım pˇr´ıkladem uvedeme numerickou nestabilitu poplatnou chyb´ am ze zaokrouhlen´ı. al˚ u Pˇ r´ıklad 1.2 Urˇceme hodnoty integr´ Z 1 1 n t In = t e dt (1.1) e 0 pro n = 1, 2, . . . Uˇzit´ım metody ”per partes” zjist´ıme, ˇze (1.2)

In = 1 − nIn−1

pˇri ˇcemˇz (1.3) Zˇrejmˇe plat´ı (1.4) jakoˇz i (1.5)

I0 =

e−1 . e

0 ≤ In ≤ In−1 ≤ 1 lim In = 0.

n→∞

2

V´ ypoˇcet na poˇc´ıtaˇci se realizuje nikoliv v term´ınech zaveden´ ych veliˇcin In avˇsak v term´ınech veliˇcin poˇc´ıtan´ ych I˜n . Vzhledem k zaokrouhlovac´ım chyb´ am zjist´ıme, ˇze existuj´ı indexy n1 a n2 tak, ˇze (1.6) I˜n1 < 0 a (1.7)

I˜n2 > 1.

Vid´ıme tedy, ˇze vztahy (1.4) neplat´ı pro veliˇciny I˜n . D˚ uvodem pro platnost vztah˚ u (1.6) a (1.7) je skuteˇcnost vyj´adˇren´ a rovnostmi I˜n = In + δn , n = 1, 2, ... z nichˇz plyne, ˇze δn = nδn−1 , odkud je patrn´e, ˇze chyba se na kaˇzd´em kroku v´ ypoˇctu n´asob´ı indexem kroku. Toto zes´ılen´ı chyby je velikosti faktori´alu, coˇz znaˇc´ı, ˇze ”destrukce” se explicitnˇe projev´ı po p kroc´ıch, kde p znaˇc´ı poˇcet cifer pouˇz´ıvan´ ych v pr˚ ubˇehu v´ ypoˇctu.

2

Diferenˇ cn´ı rovnice

Pˇripomeˇ nme si pojem obyˇcejn´e diferenci´ aln´ı rovnice. ˇ f = f (s, y, z) skal´arn´ı re´aln´ a funkce, pˇri ˇcemˇz Bud (2.1)

s ∈ I = (a, b), y ∈ S1 , z ∈ S2 ,

kde S1 ⊂ R, S2 ⊂ R, pˇri standartn´ım znaˇcen´ı mnoˇziny re´aln´ ych ˇc´ısel symbolem R. ´ Ulohu nal´ezti funkci x = x(s), s ∈ I takovou,ˇze jednak (2.2) a jednak (2.3)

x(s) ∈ S1 , x0 (s) ∈ S2 f (s, x(s), x0 (s)) = 0,

se naz´ yv´a obyˇcejnou diferenci´ aln´ı rovnic´ı 1. ˇr´ adu. Pˇ r´ıklad 2.1 Poloˇzme f (s, y, z) = z − αy. Snadno se pˇresvˇedˇc´ıme, ˇze x(s) = Cexp{αs}

3

je pro s ∈ [0, T ), kde T > 0, ˇreˇsen´ım diferenci´aln´ı rovnice x0 = αx splˇ nuj´ıc´ı poˇc´ateˇcn´ı podm´ınku x(0) = C. Pˇ r´ıklad 2.2 Nechˇt I = (−∞, +∞), S0 = S1 = [−1, 1] a f (s, y, z) = z 2 − (1 − y 2 ). Pak m´ame co do ˇcinˇen´ı s diferenci´aln´ı rovnic´ı x02 = 1 − x2 , jej´ıˇz ˇreˇsen´ı jsou d´ana formul´ı x(s) = sin(s − a), kde a ∈ R je libovoln´e. D´ale si pˇripomeˇ nme pojem diferenci´aln´ı rovnice N-t´eho ˇr´ adu. ˇ Bud f = f (s, y0 , y1 , ..., yN ) skal´arn´ı re´aln´a funkce definovan´ a pro s ∈ I ⊂ R a y0 , ..., yN leˇz´ıc´ı v podmnoˇzin´ ach re´aln´ ych ˇc´ısel S0 , ...SN . Probl´em nal´ezt funkci x = x(s) definovanou pro s ∈ I, maj´ıc´ı N derivac´ı na I a splˇ nuj´ıc´ı relace (2.4) a (2.5)

x(k) (s) ∈ Sk k = 0, 1, ..., N 0

f (s, x(s), x (s), ..., x(N ) (s)) = 0

pro vˇsechna s ∈ I, se naz´ aln´ı rovnic´ı ˇr´ adu N ; symbolicky yv´a obyˇcejnou diferenci´ p´ıˇseme f (s, x(s), x0 (s), ..., x(N ) ) = 0. Pˇ r´ıklad 2.3 Nechˇt N = 2 a I = S0 = S2 = R. Nechˇt d´ale f (s, y0 , y1 , y2 ) = y0 + y2 . Potom kaˇzd´a funkce x tvaru x(s) = Acoss + Bsins, kde A a B jsou konstanty, je ˇreˇsen´ım odpov´ıdaj´ıc´ı diferenci´aln´ı rovnice x00 + x = 0. 4

Nyn´ı pˇristoup´ıme k u ´loh´ am podobn´ ym diferenci´aln´ım rovnic´ım. Analogie bude zˇrejm´a z definice. Budeme se zab´ yvati diferenˇcn´ımi rovnicemi. ˇ Uvedme napˇred pojem diferenˇcn´ı rovnice 1. ˇr´ adu. Budiˇz d´ana soustava f = {fn (y, z)},

kde fn = fn (y, z), n = 1, 2, ..., jsou skal´ arn´ı re´aln´e funkce definovan´e na mnoˇzinˇe ´ ych ˇc´ısel a y, z ∈ S ⊂ R. Uloha nal´ezt I ⊂ Z∞ , pˇri ˇcemˇz Z∞ znaˇc´ı mnoˇzinu cel´ posloupnost {xn } ⊂ R, n ∈ I, splˇ nuj´ıc´ı (2.6) a (2.7)

xn ∈ S, xn−1 ∈ S fn (xn , xn−1 ) = 0

se naz´ yv´a diferenˇcn´ı rovnice 1. ˇr´ adu. Kaˇzd´ a posloupnost {xn } splˇ nuj´ıc´ı (2.6) a (2.7) se naz´ yv´a ˇreˇsen´ım diferenˇcn´ı rovnice (2.7). ˇ I mnoˇzina cel´ ych ˇc´ısel. Potom Pˇ r´ıklad 2.4 Bud fn (y, z) = y − z − 1 vede k diferenˇcn´ı rovnici Jej´ı ˇreˇsen´ı maj´ı tvar

xn − xn−1 = 1, xn = n + c,

kde c ∈ R je libovoln´e. Podobnˇe Pˇ r´ıklad 2.5 Nechˇt I je mnoˇzina nez´aporn´ ych cel´ ych ˇc´ısel a fn (y, z) = y − z − n, takˇze naˇse diferenˇcn´ı rovnice m´a tvar xn − xn−1 = n. Jej´ım ˇreˇsen´ım je posloupnost {xn }, kde xn =

n(n − 1) . 2

ˇ I mnoˇzina vˇsech cel´ ych ˇc´ısel a Pˇ r´ıklad 2.6 Bud fn (y, z) = y − qz, Nen´ı obt´ıˇzn´e uk´azat, ˇze ˇreˇsen´ı splˇ nuj´ıc´ı x0 = 1 odpov´ıdaj´ıc´ı rovnice xn = qxn−1 m´a tvar xn = q n . 5

Nyn´ı pˇristoup´ıme k vyˇsetˇrov´ an´ı diferenˇcn´ıch rovnic ˇr´ adu N ≥1. ˇ Bud f = {fn (y0 , y1 , ..., yN )}

posloupnost funkc´ı definovan´ ych na mnoˇzinˇe I ⊂ Z∞ , pˇr ˇcemˇz yj ∈ Sj ⊂ R, j = 0, 1, ..., N . ´ Ulohu nal´ezt posloupnost {xn }, n ∈ I a splˇ nuj´ıc´ı n´asleduj´ı poˇzadavky (2.8)

xn ∈ Sn , xn−1 ∈ Sn−1 , ...xn−N ∈ Sn−N

a (2.9)

fn (xn , xn−1 , ..., Xn−N ) = 0.

nuj´ıc´ı (2.8) a (2.9) se naz´ yv´ a ˇreˇsen´ı diferenˇcn´ı rovnice Posloupnost {xn } splˇ (2.9). Pˇ r´ıklad 2.7 Nechˇt I = Z∞ a fn (y0 , y1 , y2 ) = y0 − 2cosφy1 + y2 , ˇ sen´ı takto vznikl´e diferenˇcn´ı rovnice kde φ ∈ R. Reˇ xn − 2cosφ xn−1 + xn−2 = 0 je d´ano formul´ı xn = cosnφ. Podobnˇe jako v teorii diferenci´aln´ıch rovnic se pro nˇekter´e tˇr´ıdy u ´loh jednoznaˇcn´e ˇreˇsitelnosti dociluje vhodnou volbou poˇc´ ateˇcn´ıch podm´ınek. arn´ı. V D˚ uleˇzit´ ym pˇr´ıpadem diferenˇcn´ıch rovnic jsou diferenˇcn´ı rovnice line´ tom pˇıpadˇe (2.10)

fn (y0 , y1 , ..., yN ) = a0n y0 + a1n y1 + ... + AN n yN + bn ,

kde ajn , bn ∈ R, j = 0, 1, ..., N, n = 0, 1, .... Pˇ r´ıklad 2.8 Vˇsechny diferenˇcn´ı rovnice uveden´e v pˇr´ıkladech 2.3 − 2.7 jsou line´arn´ı. Podobnˇe diferenˇcn´ı rovnice xn + 5nxn−1 + n2 xn−2 = 2 je line´arn´ı, zat´ımco je neline´arn´ı.

xn−1 − 2x2n−1 = 0

6

Analogicky jako v teorii diferenci´aln´ich rovnic, obecn´e ˇreˇsen´ı diferenˇcn´ıch rovnic um´ıme plnˇe charakterizovat pro pˇr´ıpad diferenˇcn´ıch rovnic line´arn´ıch. Zab´ yvejme se u ´lohou (analytick´eho) sestrojen´ı ˇreˇsen´ı diferenˇcn´ıch rovnic 1. ˇr´adu. Tedy, ˇreˇsme rovnici (2.11)

+ xn = an xn−1 + bn , n ∈ I ⊂ Z∞ , an 6= 0,

+ znaˇc´ı mnoˇzinu vˇsech nez´aporn´ ych cel´ ych ˇc´ısel. kde Z∞

Vyˇsetˇrujme nejprve homogenn´ı diferenˇcn´ı rovnici (2.12)

xn = an xn−1 .

Snadno zsjist´ıme, ˇze ˇreˇsen´ı splˇ nuj´ıc´ı podm´ınku (2.13)

x0 = c

je d´ano v´ yrazem (2.14)

xn = cπn ,

kde (2.15)

π0 = 1, πn =

n Y

ak , n = 1, 2, ...

k=1

K urˇcen´ı ˇreˇsen´ı nehomogenn´ı rovnice (2.12) splˇ nuj´ıc´ı x0 = c, pouˇzijeme metody zn´am´e z teorie line´arn´ıch diferenci´aln´ıch rovnic pod n´azvem metoda variace konstanty. Poloˇzme (2.16)

xn = cn πn ,

kde {cn } podl´eh´a urˇcen´ı. Z rovnost´ı x0 = c0 π0 = c0 , vid´ıme, ˇze c0 = c. Po dosazen´ı (2.16) do (2.11) zjist´ıme, ˇze cn πn = an cn−1 πn−1 + bn = cn−1 πn + bn . Z pˇredpokladu an 6= 0 plyne, ˇze t´eˇz πn = 6 0. Pˇredchoz´ı vztahy implikuj´ı rovnost cn = cn−1 + 7

bn πn

a d´ale pak rovnosti cn = c0 +

n X

k=1

(ck − ck−1 ) = c +

n X bk . πk

k=1

V´ ysledek shrneme ve tvaru vˇety. ˇ sen´ı diferenˇcn´ı rovnice Vˇ eta 2.1 Nechˇt an 6= 0, n = 1, 2, ... Reˇ xn = an xn−1 + bn , splˇ nuj´ıc´ı x0 = c , je d´ ano v´yrazem   (2.17)

x n = πn

n X bk c+ πk k=1

pˇri ˇcemˇz π0 = 1, πn =

!

, n = 0.1, ...

n Y

ak .

k=1

ˇ I mnoˇzina nez´aporn´ Pˇ r´ıklad 2.9 Bud ych ˇc´ısel, fn (y0 , y1 , y2 ) = y0 − y1 − y2 . Odpov´ıdaj´ıc´ı diferenˇcn´ı rovnice (2.18)

xn − xn−1 − xn−2 = 0

definuje Fibonacciova ˇc´ısla (x0 = 0, x1 = 1). Podobnˇe jako v problematice diferenci´aln´ıch rovnic je jedn´ım z moˇzn´ ych zp˚ usob˚ u ˇreˇsen´ı line´arn´ıch rovnic ˇr´ adu N jejich pˇrevod na soustavy line´arn´ıch rovnic 1. ˇr´adu. Uk´aˇzeme si zm´ınˇen´ y postup na pˇr´ıkladˇe rovnice (2.18). Nechˇt Xn =

’

xn xn−1

“

, X0 =

’

1 0

a

A=

’

1 1 1 0 8

“

“

, n = 1, 2, ...

Snadno nahl´edneme, ˇze rovnici (2.18) lze vyj´adˇriti ve tvaru ’ “ 1 Xn = AXn−1 , X0 = (2.19) . 0 Podobnˇe jako v pˇr´ıpadˇe line´arn´ı diferenˇcn´ı rovnice 1. ˇr´ adu, ˇreˇsen´ı rovnice (2.19) m´a tvar Xn = An X0 . Odtud plyne platnost zaj´ımav´eho vztahu √ 1+ 5 xn = . n→∞ xn−1 2 lim

Vid´ıme tedy, ˇze limitn´ı pomˇer veliˇcin xn a xn−1 je d´an ˇc´ıslem charakterizuj´ıc´ım t. zv. zlat´y ˇrez. Protoˇze rovnice (2.18) m´a vskutku zaj´ımavou biologickou interpretaci, m˚ uˇzeme ve v´ yˇse uveden´em v´ ysledku spatˇrovati u ´kaz estetiky projevuj´ıc´ı se v nˇekter´ ych oblastech ˇzivota na naˇs´ı planetˇe. Model popisuj´ıc´ı jist´ y zp˚ usob rozmnoˇzov´an´ı kr´al´ık˚ u poch´ az´ı ze 13. stolet´ı; jeho autorem je pr´avˇe Fibonacci (r. 1228). Zcela jinou aplikac´ı line´arn´ıch diferenˇcn´ıch rovnic je n´asleduj´ıc´ı u ´loha. ˇ Bud n X pn (x) = bk xn−k , b0 6= 0, n = 0, 1, ..., N k=0

a poloˇzme si u ´lohu stanovit hodnotu

pn (z), kde z je dan´ y bod na re´aln´e ose. Naˇs´ım c´ılem je pokud moˇzno minimalitovat pˇri tom poˇcet aritmetick´ ych operac´ı ve snaze sn´ıˇzit riziko numerick´e nestability. Algoritmus 2.1 Poˇc´ıtejme veliˇciny x0 , x1 , ..., xN rekurentnˇe pomoc´ı relac´ı (2.20)

x0 = b0 , xn = zxn−1 + bn , n = 1, ..., N.

Posloupnost {xn } dan´a v (2.20) je ˇreˇsen´ı diferenˇcn´ı rovnice (2.11) splˇ nuj´ıc´ı podm´ınku x0 = b0 , kde an = z pro n = 1, 2, ..., . V tomto pˇr´ıpadˇe πn = z n

9

a tud´ıˇz, na z´akladˇe vˇety 2.1, xn =

n X

bk z k

k=0

a, speci´alnˇe pro n = N , xN =

N X

bk z k = pN (z).

k=0

Plat´ı tedy Vˇ eta 2.2 Veliˇciny xn tvoˇren´e pomoc´ı algoritmu 2.1 jsou hodnoty polynom˚ u pn definovan´ych pomoc´ı pn (x) =

n X

bk xk , n = 0, 1, ..., N

k=0

v bodˇe x = z. Algoritmus 2.1 se naz´ yv´a Hornerov´ym sch´ematem. Nyn´ı si uk´aˇzeme jak lze uveden´ y algoritmus zobecnit pro v´ ypoˇcet hodnot derivac´ı polynomu pn v bodˇe x = z. Opˇet je nutn´e zd˚ uraznit, ˇze pˇr´ım´ a aplikace standardn´ıch pravidel derivov´ an´ı ˇclen po ˇclenu nen´ı vhodn´a a algoritmy typu Hornerova jsou daleko stabilnˇejˇs´ı a jednoduˇsˇs´ı z hlediska teorie sloˇzitosti. Vyjdeme ze zn´am´ ych vztah˚ u p(x) = p(z + h) =

N X

ck hN −k ,

k=0

kde (2.21)

cN −k =

1 (k) p (z), k = 0, 1, ..., N. k! (0)

Algoritmus 2.2 Oznaˇcme prvky tvoˇren´e pomoc´ı algoritmu 2.1 symboly {xn }. (k) Pro k = 1, 2, ..., N tvoˇrme posloupnosti {xn } rekurzivnˇe pomoc´ı diferenˇcn´ıch sch´emat (2.22)

(k)

(k−1)

x0 = x0

(k)

(k−1) , x(k) , n = 0, 1, ..., N − k. n = zxn−1 + xn

Vˇ eta 2.3 Oznaˇcuje-li pn polynom pn (x) =

n X

bk xn−k , n = 0, 1, ..., N,

k=0

10

(0)

x0

(1)

x0

(2) x0 (3) x0 (4) x0

(0)

x1

(1)

x1

(2) x1 (3) x1

(0)

x2 ↓ (1) x2

(2) x2

(0)

x3

(0)

x4

(1)

x3

Figure 1: Obr´azek ˇc. 2.1 pak veliˇciny urˇcovan´e algoritmem 2.2 splˇ nuj´ı (2.23)

x(k) n =

1 (k) p (z), n = 1, , ..., N, k = 0, ..., N − n. k! n+k

Tud´ıˇz koeficienty cN −k poˇzadovan´e ve formuli (2.21) jsou pr´avˇe ty, jeˇz se nal´ezaj´ı v doln´ı diagon´ale na obr´azku 2.1 pro N = 4 v posloupnosti veliˇcin tvoˇren´ ych algoritmem 2.2.

3 3.1

Iteraˇ cn´ı metody ˇ reˇ sen´ı neline´ arn´ıch rovnic a jejich soustav Postupn´ e aproximace

Zaˇcnˇeme pˇr´ıkladem u ´lohy naj´ıt koˇren funkce g, t. j. nal´ezt takov´ y bod x ˆ∈I= [a, b], −∞ < a ≤ b < +∞, tak, aby g(ˆ x) = 0 za pˇredpokladu, ˇze g je dan´a spojit´a funkce zobrazuj´ıc´ı I = I1 do I2 ⊂ R. Speci´alnˇe, g m˚ uˇze b´ yt polynom. A jiˇz v tomto speci´aln´ım pˇr´ıpadˇe si m˚ uˇzeme uvˇedomit nˇekolik d˚ uleˇzit´ ych skuteˇcnost´ı. Tak pˇredevˇs´ım, obecnˇe nemus´ı v I existovat koˇren funkce g. Ale i kdyˇz takov´ y koˇren existuje, nemus´ı existovat analytick´e vyj´adˇren´ı tohoto koˇrene pomoc´ı parametr˚ u charakterizuj´ıc´ıch g (na pˇr. koeficienty polynomu); takov´e formule jsou zn´amy pro pˇr´ıpad polynom˚ u stupˇ n˚ u nepˇrevyˇsuj´ıc´ıch 4. A jak v´ıme, cena takov´ ych formul´ı je sporn´a, protoˇze takov´e formule jsou mnohdy prakticky nepouˇziteln´e (t. zv. casus irreducibilis). D˚ usledky tˇechto skuteˇcnost´ı jsou zˇrejm´e: Nadˇeji, ˇze se podaˇr´ı nal´ezt koˇren g maj´ı sp´ıˇse numerick´e neˇz analytick´e metody. Mezi numerick´ ymi metodami se z d˚ uvod˚ u snadn´e algoritmick´e realizovatelnosti prosadily metody iteraˇcn´ı a to navzdory jejich neunivers´ alnosti a pohˇr´ıchu pomˇernˇe pomal´e konvergenci. Pr´avˇe uveden´ y nedostatek dal podnˇet k rozvoji metod urychlov´an´ı konvergence obecn´ ych posloupnost´ı, jmenovitˇe pak takov´ ych

11

generovan´ ych iteraˇcn´ımi procesy. Z´akladn´ım iteraˇcn´ım postupem jsou postupn´e aproximace. Jsou urˇceny k sestrojov´an´ı aproximac´ı ˇreˇsen´ı rovnic typu s = f (s),

(3.1)

kde f je dan´a funkce, tedy, g(s) = s − f (s). Algoritmus 3.1 Zvolme x0 libovolnˇe a tvoˇrme posloupnost {xn } rekurz´ıvnˇe pomoc´ı relace (3.2) xn = f (xn−1 ). Vˇs´ımnˇeme si, ˇze k tomu, abychom mohli prov´ adˇeti jednotliv´e kroky dle (3.2) mus´ı f kromˇe spojitosti m´ıti jeˇstˇe dalˇs´ı vhodn´e vlastnosti. Budeme pˇredpokl´adat, ˇze obor hodnot funkce f patˇr´ı do I1 ; k tomu staˇc´ı, aby I2 ⊂ I1 . V takov´em pˇr´ıpadˇe spojitost f nejen zaruˇcuje neomezenou proveditelnost procesu (3.2) ale t´eˇz existenci hledan´eho koˇrene v I1 . Abychom to ozˇrejmili, poloˇzme g(s) = s − f (s), a ≤ s ≤ b, takˇze g(a) ≤ 0, a g(b) ≥ 0.

Ze spojitosti f plyne spojitost g a odtud pak skuteˇcnost, ˇze g nab´ yv´ a vˇsech hodnot z intervalu [g(a), g(b)] = I3 ⊂ I1 . Existuje tedy alespoˇ n jeden bod sˆ takov´ y, g(ˆ s) = 0. To znaˇc´ı, ˇze sˆ = f (ˆ s) a sˆ je n´ami hledan´ y koˇren. Obecnˇe, mnoˇzina M vˇsech samodruˇzn´ ych bod˚ u funkce f , t. j. M = {s ∈ I1 : s = f (s)} m´a mohutnost moh M ≥ 1. Jednoznaˇcnosti se dociluje pomoc´ı dalˇs´ıch pˇredpoklad˚ u. Definice 3.1 Pˇredpokl´adejme, ˇze f zobrazuje I ⊂ R do R a ˇze existuje konstanta L takov´ a, ˇze nerovnost (3.3)

|f (s) − f (t)| ≤ L|s − t|

plat´ı pro libovoln´a s, t ∈ I . V takov´em pˇr´ıpadˇe se f naz´ yv´ a lipschitzovskou a konstanta L Lipschitzovou konstantou funkce f. 12

Vˇs´ımnˇeme si, ˇze funkce lipschitzovsk´ a na I je nutnˇe spojit´a na I. Je - li f diferencovateln´a na I a plat´ı sup {|f 0 (s)| : s ∈ I} ≤ L, kde L > 0 je konstanta, pak je f na I lipschitzovsk´ a a L je Lipschitzovou konstantou funkce f . Pˇredpokl´adejme, ˇze Lipschitzova konstanta L splˇ nuje nerovnost 0 < L < 1.

(3.4)

Snadno nahl´edneme, ˇze pro s1 , s2 ∈ M plat´ı |s1 − s2 | = |f (s1 ) − f (s2 )| ≤ L|s1 − s2 | a tud´ıˇz s1 = s2 . M˚ uˇzeme pˇristoupit k formulaci z´akladn´ıho v´ ysledku o konvergenci postupn´ ych aproximac´ı. Vˇ eta 3.1 Budˇ I = [a, b] uzavˇren´y koneˇcn´y interval a nechˇt funkce f splˇ nuje n´ asleduj´ıc´ı podm´ınky (i) − (iii). (i) f (s) ∈ I pro s ∈ I. a s Lipschitzovou konstantou L < 1 . (ii) f je lipschitzovsk´ Potom pro libovoln´e nulov´e pˇribl´ıˇzen´ı x0 ∈ I posloupnost {xn } definovan´ a pomoc´ı algoritmu 3.1 konverguje k (jedin´emu) ˇreˇsen´ı rovnice s = f (s). D˚ ukaz. Z pˇredchoz´ıch u ´vah jiˇz v´ıme, ˇze existuje alespoˇ n jedno ˇreˇsen´ı rovnice s = f (s) v I a ˇze, d´ıky pˇredpokladu o L, toto ˇreˇsen´ı je urˇceno jednoznaˇcnˇe. Staˇc´ı tedy dok´azati konvergenci posloupnosti {xn }. Nechˇt x ˆ = f (ˆ x). Zˇrejmˇe xn − x ˆ = f (xn−1 ) − f (ˆ x), takˇze a tud´ıˇz (3.5) Protoˇze

ˆ| ≤ L|xn−1 − x |xn − x ˆ| ˆ|. ˆ| ≤ Ln |x0 − x |xn − x lim Ln = 0,

n→∞

obdrˇz´ıme ˇz´adan´ y v´ ysledek ˆ| = 0. lim |xn − x

n→∞

||| 13

Pˇ r´ıklad 3.1 Pomoc´ı postupn´ ych aproximac´ı hledejme ˇreˇsen´ı rovnice s = e−s . Protoˇze e−s > 0 pro s ∈ (−∞, +∞) a s > 1 ≥ e−s pro s > 1, staˇc´ı omezit se na interval [0, 1]. Snadno si uvˇedom´ıme, ˇze obor hodnot funkce f (s) = e−s pro s ∈ [0, 1] leˇz´ı v intervalu [e−1 , 1] ⊂ [0, 1] a ˇze pro s1 s2 ∈ [0, 1] plat´ı |f (s1 ) − f (s2 )| = |f 0 (˜ s)||s1 − s2 |, pˇri ˇcemˇz f 0 (˜ s) = −e−˜s , s˜ ∈ [s1 , s2 ]. Protoˇze max {|f 0 (s)| : s ∈ [0, 1]} = 1, Lipschitzova konstanta pro f na [0, 1] je rovna 1, coˇz k platnosti vˇety 3.1 nestaˇc´ı. Avˇsak zvol´ıme-li za I interval [ 12 , log 2], snadno zjist´ıme, ˇze pro 12 ≤ s ≤ log 2, 1 1 = e−log2 ≤ e−s ≤ e− 2 < log 2 2 a tud´ıˇz [ 12 , log 2] = I ⊂ [0, 1]. Nav´ıc prvn´ı derivace funkce f , kde f (s) = e−s , je klesaj´ıc´ı, takˇze 1 max {|f 0 (s)| : s ∈ I} = f 0 ( ) ≈ 0, 606531 < 1. 2 Iteraˇcn´ı proces (3.2) pro funkci f (s) = e−s je tedy konvergentn´ı v I. Odhad (3.5) nem´a bezprostˇredn´ı praktick´ y v´ yznam, z´avis´ı totiˇz na znalosti hledan´eho ˇreˇsen´ı. V praxi jsme totiˇz vˇzdy nuceni omezit se na urˇcit´ y koneˇcn´ y poˇcet iterac´ı xn v (3.2). R´adi bychom proto znali nˇejak´ y realistick´ y odhad chyby po n - t´em kroku. Vˇs´ımnˇeme si, ˇze pro libovoln´e k ≥ 1 plat´ı vztahy |xk+1 − xk |

= |f (xk ) − f (xk−1 )| ≤ L|xk − xk−1 | ≤ Lk |x1 − x0 |.

ˇ n zafixov´ano pevnˇe a nechˇt m = n + p > n. Potom Bud xm − xn =

p−1 X

(xk+1 − xk ).

k=n

14

Po aplikaci troj´ uheln´ıkov´e nerovnosti obdrˇz´ıme vztahy Pn+p−1 |xm − xn | ≤ k=n Lk |x1 − x0 | = Ln

≤ Ln

Pp−1

k=0

P∞

k=0

Lk |x1 − x0 | Lk |x1 − x0 |

= Ln (1 − L)−1 |x1 − x0 |. Protoˇze pro m → ∞ plyne, ˇze xm → x ˆ, tedy t´eˇz p → ∞, a plat´ı (3.6)

|xn − x ˆ| ≤

Ln |x1 − x0 |. 1−L

u vˇety 3.1 je odhad chyby po n kroc´ıch D˚ usledek 3.1 Pˇri splnˇen´ı pˇredpoklad˚ postupn´ych aproximac´ı definovan´ych pomoc´ı algoritmu 3.1 d´ an v´yrazem (3.6). Zkoumejme bl´ıˇze chov´an´ı posloupnosti {xn } za pˇredpoklad˚ u vˇety 3.1 doplnˇen´ ych o pˇredpoklad n´asleduj´ıc´ı. Nechˇt f je spojitˇe diferencovateln´ a a nechˇt v cel´em intervalu I je tato derivace nenulov´ a. ˇ kles´a nebo roste v I. To znamen´a, ˇze f bud ˆ = f (ˆ x), pak f (xn ) 6= x Je-li x0 6= x ˆ pro libovoln´e n = 1, .2, ... Abychom to nahl´edli, staˇc´ı si uvˇedomit, ˇze kdyby xn = f (xn−1 ) = f (xn ), pˇri ˇcemˇz xn 6= xn−1 , pak 0 = f (xn−1 ) − f (xn ) = f 0 (ξ) (xn−1 − xn ),

kde ξ leˇz´ı mezi xn−1 a xn , t.j.

(xn−1 − ξ) (ξ − xn ) < 0.

Protoˇze xn−1 − xn = 6 0, mus´ı b´ yt f 0 (ξ) = 0, ale to je ve sporu s pˇredpokladem. Odtud vypl´ yv´ a, ˇze chyba dn = xn − x ˆ nen´ı pro ˇz´adn´e n ≥ 0 nulov´ a. Existuje limita

lim

n→∞

dn+1 dn

a kdyˇz ano, pak ˇcemu je rovna? 15

Snadno zjist´ıme, ˇze dn+1 = xn+1 − x ˆ = f (xn ) − f (ˆ x) x + θn dn ) dn , = f (ˆ x + dn ) − f (ˆ x) = f 0 (ˆ

kde

Definujme εn pomoc´ı relace

0 < |θn | < 1.

x + θn dn ) = f 0 (ˆ x ) + εn . f 0 (ˆ Potom (3.7)

x) + εn )dn dn+1 = (f 0 (ˆ

a tud´ıˇz, jeˇzto εn → 0 pro n → ∞ plyne na z´akladˇe spojitosti f 0 rovnost lim

(3.8)

n→∞

dn+1 = f 0 (ˆ x). dn

x) Lze t´eto skuteˇcnosti vyuˇz´ıt pro praktick´e u ´ˇcely? Veliˇcinu x ˆ a tud´ıˇz i f 0 (ˆ nezn´ame. Uk´aˇzeme, ˇze explicitn´ı znalost uveden´ ych veliˇcin nen´ı nutn´ a k tomu, abychom rovnosti (3.8) efektivnˇe vyuˇzili k sestrojen´ı algoritmu, kter´ y poskytne aproximace konverguj´ıc´ı k hledan´emu ˇreˇsen´ı x ˆ rychleji neˇz posloupnost {xn }. Pˇredpokl´adejme na okamˇzik, ˇze ve formuli (3.7) εn = 0 pro nˇejak´e pevn´e n ≥ 1. Nechˇt f 0 (ˆ x) = A, takˇze xn+1 − x ˆ = A (xn − x ˆ),

Snadno zjist´ıme, ˇze

ˆ = A (xn+1 − x ˆ). xn+2 − x A=

takˇze x ˆ= = xn +

xn+2 − xn+1 , xn+1 − xn

1 (xn+1 − Axn ) 1−A

1 (xn+1 − xn )2 (xn+1 − xn ) = xn − 1−A xn+2 − 2xn+1 + xn

Vid´ıme tedy, ˇze za pˇredpokladu, ˇze εn = 0, pro nˇejak´e n ≥ 1, lze pˇresn´e ˇreˇsen´ı ˆ obdrˇzeti pouˇzit´ım tˇr´ı po sobˇe jdouc´ıch iterac´ı. To vˇsak je dost akademick´ x a situace. Aˇckoliv veliˇciny εn 6= 0, lze oˇcek´ avat, ˇze ve srovn´ an´ı s modulem f 0 (ˆ x) jsou mal´e. 16

Pro velk´a n tak posloupnost (3.9)

yn = xn −

(xn+1 − xn )2 xn+2 − 2xn+1 + xn

sk´ yt´a pro x ˆ lepˇs´ı aproximaci neˇz xn . M´a tedy smysl vyˇsetˇrovat Algoritmus 3.2 Budˇ {xn } libovoln´ a posloupnost a nechˇt yn je tvoˇreno pomoc´ı (3.9). ˇ Zavedme oznaˇcen´ı ∆xn = xn+1 − xn , n = 0, 1, ... a Na pˇr

∆k+1 xn = ∆(∆k xn ), k ≥ 1. ∆2 xn = ∆(xn+1 − xn ) = xn+2 − 2xn+1 + xn

Formuli (3.9) lze pak ps´at ve tvaru

yn = xn −

(3.10)

(∆xn )2 . ∆ 2 xn

Posloupnost {yn } tvoˇren´ a pomoc´ı algoritmu 3.2 definuje t. zv. Aitken˚ uv urychlovac´ı ∆2 - proces. Jest na m´ıstˇe zd˚ uraznit, ˇze Aitken˚ uv proces je definov´ an pro posloupnosti ne nutnˇe generovan´e pomoc´ı algoritmu 3.1. Pro takto obecn´ y proces plat´ı pomˇernˇe siln´e tvrzen´ı. a posloupnost takov´ a, ˇze Vˇ eta 3.2 Budˇ {xn } dan´ x ˆ = lim xn . n→∞

Nechˇt (3.11) a (3.12)

dn = xn − x ˆ 6= 0 n = 0, 1, ... dn+1 = (A + εn ) dn ,

kde A je konstanta takov´ a, ˇze a

|A| < 1 εn → 0 pro n → ∞. 17

Potom posloupnost {yn } tvoˇren´ a pomoc´ı (3.9) v algoritmu 3.2 je dobˇre definovan´ a pro dostateˇcnˇe velk´ a n a nav´ıc ˆ yn − x = 0, n→∞ xn − x ˆ lim

(3.13)

ˆ rychleji neˇz posloupnost {xn }. coˇz znaˇc´ı, ˇze posloupnost {yn } konverguje k x D˚ ukaz. Pouˇzit´ım formule (3.12) obdrˇz´ime

dn+2 = (A + εn+1 ) (A + εn ) dn , takˇze

∆ 2 xn

kde

= xn+2 − 2xn+1 + xn = ‚dn+2 − 2dn+1 + ƒ dn = (A − 1)2 + εn0 dn ,

ε0n = A[εn + εn+1 ] − 2εn + εn εn+1 .

Z podm´ınky εn → 0 pro n → ∞ plyne, ˇze

ε0n → 0 pro n → ∞.

(3.14) Odtud dost´av´ame,

(A − 1)2 + ε0n 6= 0

pro dostateˇcnˇe velk´a n, ˇreknˇeme pro n > n0 . Tud´ıˇz ∆2 xn 6= 0 pro n > n0 . D´ale pak a tud´ıˇz

∆xn = ∆dn = (A − 1 + εn ) dn yn − x ˆ = dn − = dn −

(∆xn )2 ∆2 xn

(A − 1 + εn )2 dn . (A − 1)2 + ε0n

Z εn → 0 pro n → ∞ a (3.14) plyne, ˇze

ˆ yn − x ε0 − 2(A − 1) εn − ε2n = n → 0, dn (A − 1)2 + ε0n coˇz jsme mˇeli dok´azati. ||| D˚ usledek 3.2 Pˇredpokl´ adejme, ˇze f splˇ nuje podm´ınky vˇety 3.1 a nav´ıc, ˇze f m´ a spojitou 1. derivaci na I a f 0 (s) 6= 0 pro s ∈ I. Pro x0 6= x ˆ posloupnost {xn } vyhovuje pˇredpoklad˚ um vˇety 3.2, takˇze algoritmus 3.2 poskytuje urychlen´ı konvergence. 18

D˚ ukaz. Zˇrejmˇe staˇc´ı ovˇeˇrit platnost vztahu (3.12). To vˇsak je vztah (3.7) odvozen´ y za obecnˇejˇs´ıch pˇredpoklad˚ u. Zb´ yv´ a tud´ıˇz ovˇeˇrit poˇzadavek |A| < 1. Platnost tohoto vztahu vˇsak plyne odtud, ˇze proces postupn´ ych aproximac´ı je konvergentn´ı. T´ımto konstatov´ an´ım m˚ uˇzeme d˚ ukaz ukonˇcit. ||| Znovu pˇripomeˇ nme ˇsirˇs´ı aplikabilitu urychlovac´ıho Aitkenova ∆2 - procesu i mimo oblast metody postupn´ ych aproximac´ı.

3.2

Kvadratick´ a konvergence a Newtonova metoda

6 0 pro s ∈ I a obdrˇzeli jsme V pˇredchoz´ım odstavci jsme pˇredpokl´adali, ˇze f 0 (s) = tak charakteristiku rychlosti konvergence danou pomoc´ı (3.12), t. zv. line´ arn´ı konvergenci. D´a se ˇr´ıci, ˇze v takov´em pˇr´ıpadˇe poˇcet platn´ ych cifer je line´arn´ı funkc´ı indexu iterac´ı. Uˇciˇ nme nyn´ı pˇredpoklad (3.15) f 0 (ˆ x) = 0. Tento pˇredpoklad je pomˇernˇe velmi siln´ y a zabezpeˇcuje splnˇen´ı nˇekter´ ych pˇredpoklad˚ u vˇety 3.1. Vˇ eta 3.3 Budˇ I ⊂ (−∞, +∞), (nikoliv nutnˇe ohraniˇcen´y interval) a nechˇt f je definovan´ a na I a pˇritom plat´ı (i) f a f 0 jsou spojit´e na I. (ii) Rovnice s = f (s) m´ a ˇreˇsen´ı x ˆ ∈ I takov´e, ˇze plat´ı (3.15). Potom existuje δ > 0 takov´e, ˇze algoritmus 3.1 poskytuje posloupnost {xn } ˆ| ≤ δ. konverguj´ıc´ı k x ˆ pro vˇsechna |x0 − x D˚ ukaz. Oznaˇcme x − δ, x ˆ + δ]. Iδ = [ˆ

ˇ 0 < L < 1 . Potom existuje δ0 > 0 a Pro dostateˇcnˇe mal´a δ1 > 0, Iδ1 ⊂ I. Bud L > 0 takov´e, ˇze |f 0 (s)| ≤ L pro s ∈ Iδ0 . Nechˇt δ = min{δ0 , δ1 }. Podle vˇety o stˇredn´ı hodnotˇe,

|f (s) − x ˆ| = |f (s) − f (ˆ x)| ≤ |f 0 (ξ)||s − x ˆ| ≤ L|s − x ˆ| < δ. Tud´ıˇz obor hodnot funkce f je obsaˇzen v Iδ a algoritmus 3.1 vede ke konvergenci podle vˇety 3.1. ||| Zesilme naˇse poˇzadavky na hladkost funkce f pˇredpokladem existence spojit´e druh´e derivace f 00 na Iδ a jej´ı nenulovosti. ˆ implikuje platnost Podobnˇe jako v pˇredchoz´ıch u ´vah´ ach pˇredpoklad x0 6= x vztah˚ u xn 6= x ˆ, n = 1, 2, ... Postupn´e aproximace nemohou tedy poskytnout 19

pˇresn´e ˇreˇsen´ı v koneˇcn´em poˇctu krok˚ u. Uˇzit´ım Taylorova rozvoje se zbytkem, obdrˇz´ıme vztahy ˆ dn+1 = xn+1 − x = f (xn ) − f (ˆ x) 1 = f 0 (ˆ x + θn dn )d2n , x)dn + f 00 (ˆ 2 kde 0 < |θn | < 1. Z naˇsich pˇredpoklad˚ u plyne, ˇze (3.16)

dn+1 =

1 00 f (ˆ x + θn dn )d2n . 2

Protoˇze dn 6= 0, n = 1, 2, ... a dn → 0 pro n → ∞, odvod´ıme, ˇze (3.17)

1 dn+1 x). = f 00 (ˆ 2 n→∞ dn 2 lim

To je pozoruhodn´ y vztah, podle nˇehoˇz chyba po (1+n)- t´em kroku je u ´mˇern´ a a ˇctverci chyby po n - t´em kroku. Tato skuteˇcnost se oznaˇcuje jako kvadratick´ konvergenvce. Poˇcet platn´ ych cifer se zdvojuje po kaˇzd´em iteraˇcn´ım kroku. ˇ a > 0 a nechˇt f (s) = (1/2)[s + (a/s)] pro s > 0. Rovnice Pˇ r´ıklad 3.2 Bud √ √ √ s = f (s) m´a ˇreˇsen´ı a. Zˇrejmˇe f 0 ( a) = 0 a f 00 (s) > 0 pro s > 0 a f 00 ( a) = √1 . Tud´ ıˇz a) ’ “ 1 a (3.18) xn + xn+1 = 2 xn √ √ a k a. Pozdˇeji uk´aˇzeme, konverguje k a kvadraticky pro x0 dostateˇcnˇe bl´ızk´ ˇze tento proces konverguje pro libovoln´e v´ ychoz´ı pˇribl´ıˇzen´ı x0 > 0. Nyn´ı pˇristoup´ıme k snad nejpouˇz´ıvanˇejˇs´ı metodˇe pˇribliˇzn´eho ˇreˇsen´ı neline´arn´ıch rovnic a lze ˇr´ıci, ˇze nejen rovnic ale neline´arn´ıch u ´loh - k metodˇe Newtonovˇe. V pˇredchoz´ıch u ´vah´ach jsme dospˇeli ke kvadratick´e konvergenci za v´ yraznˇe akademick´eho pˇredpokladu nulovosti prvn´ı derivace v pˇresn´em ˇreˇsen´ı. Ani se nechce vˇeˇrit, ˇze tento akademismus m´a tak n´adhern´e uplatnˇen´ı, k jehoˇz v´ ykladu ˇ pˇristupujeme. ted ˇ F definovan´a a dvakr´ Bud at spojitˇe diferencovateln´ a na intervalu I = [a, b] ⊂ R a nechˇt F 0 (s) 6= 0 pro s ∈ I. D´ale nechˇt rovnice (3.19)

F (s) = 0

m´a ˇreˇsen´ı x ˆ ∈ (a, b) (nikoliv nutnˇe jedin´e). D´a se uk´azat, ˇze toto ˇreˇsen´ı lze nal´ezt pomoc´ı postupn´ ych aproximac´ı s funkc´ı f (s) = s + M F (s), 20

kde M je konstanta splˇ nuj´ıc´ı urˇcit´e podm´ınky, na pˇr. M = (α + β)/2, kde I = [0, 1], F (0) < 0 < F (1) a 0 < α ≤ F 0 (s) ≤ β. V´ıme t´eˇz, ˇze genericky obdrˇz´ıme line´arn´ı konvergenci. Jak to tedy zaˇr´ıdit, abychom sestrojili algoritmus poskytuj´ıc´ı kvadratickou konvergenci? Coˇz pˇripustit m´ısto konstantn´ıho faktoru M vhodnou funkci, tedy snaˇzit se sestrojit f ve tvaru f (s) = s + h(s)F (s), kde hodnoty funkce h se ”snadno” poˇc´ıtaj´ı. Snaˇzme se tedy splnit tyto poˇzadavky: x x) a f 0 (ˆ x) = 0, ˆ = f (ˆ To n´as vede k relac´ım f 0 (s) = 1 + h0 (s)F (s) + h(s)F 0 (s) a tud´ıˇz a jeˇzto F (ˆ x) = 0,

x) = −1 h0 (ˆ x) + h(ˆ x)F (ˆ x)F 0 (ˆ h(s) = −

1 F 0 (s)

.

Dospˇeli jsme tedy k n´asleduj´ıc´ımu algoritmu. Algoritmus 3.3 Volme x0 a urˇceme posloupnost {xn } pomoc´ı vztah˚ u (3.20)

xn+1 = xn −

1 F (xn ), n = 0, 1, ... F 0 (xn )

Algoritmus 3.3 nese n´azev Newtonova metoda, ˇci metoda Newtonova - Raphsonova. Z pˇredchoz´ıho v´ ykladu jiˇz v´ıme, ˇze, jeˇzto f (s) = (1/F 0 (s))F (s), f 0 (ˆ x) = 0, takˇze pˇr´ıpadn´a konvergence posloupnosti {xn } je kvadratick´ a. Jednotliv´e kroky Newtonovy metody maj´ı velice n´azornou interpretaci. V bodˇe xn aproximuje graf funkce F pomoc´ı teˇcny ke grafu F v tomto bodˇe. Pr˚ useˇc´ık t´eto teˇcny s osou x urˇcuje dalˇs´ı aproximaci xn+1 . Um´ıme si snadno pˇredstavit grafy funkce F , pro nˇeˇz Newton˚ uv proces diverguje. (To je velice hezk´e cviˇcen´ı.) Podobnˇe jako pro vˇsechny iteraˇcn´ı metody typu postupn´ ych aproximac´ı, slabina Newtonovy metody spoˇc´ıv´ a mimo jin´e t´eˇz v siln´e z´avislosti na volbˇe poˇc´ateˇcn´ıho pˇribl´ıˇzen´ı. Uk´aˇzeme si jak lze zabezpeˇcit glob´aln´ı konvergenci; poˇzadujeme vyˇsˇs´ı hladkost F - totiˇz znam´enkovou st´alost druh´e derivace F 00 . 21

Vˇ eta 3.4 Pˇredpokl´ adejme, ˇze funkce F je definovan´ a a dvakr´ at spojitˇe diferencovateln´ a v I = [a, b] a nechˇt splˇ nuje tyto poˇzadavky (i) F (a)F (b) < 0; (ii) F 0 (s) 6= 0 pro s ∈ I; (iii) Budˇ F 00 (s) ≥ 0 s ∈ I nebo F 00 (s) ≤ 0, s ∈ I; (iv) Nechˇt (α) c = a jestliˇze |F 0 (a)| ≤ |F 0 (b)| , (β) c = b, jestliˇze |F 0 (b)| ≤ |F 0 (a)| . D´ ale nechˇt

Œ Œ Œ Œ Œ F (c) Œ ≤ b − a. Œ F 0 (c) Œ

Potom Newtonova metoda konverguje pro libovoln´e v´ychoz´ı pˇribl´ıˇzen´ı x0 k jedin´emu ˇreˇsen´ı rovnice F (s) = 0. Pozn´ amky. Podm´ınka (i) zˇrejmˇe zaruˇcuje existenci koˇrene d´ıky spojitosti ˇ F . Z (ii) plyne jednoznaˇcnost koˇrene v I. Podm´ınka (iii) ˇr´ık´ a, ˇze F je bud ˇ konk´ avn´ı nebo konvexn´ı v I. Posl´eze (iv) zajiˇstuje, aby teˇcna ke kˇrivce y = F (x) v koncov´em bodˇe, kde |F 0 (x)| je menˇs´ı, prot´ınala osu x v intervalu I. D˚ ukaz vˇ ety 3.4. Vˇeta 3.4 zahrnuje tyto ˇctyˇri odliˇsn´e situace (a) F (a) < 0, F (b) > 0, F 00 (s) ≤ 0 (c = b), (b) F (a) > 0, F (b) < 0, F 00 (s) ≥ 0 (c = b), (c) F (a) < 0, F (b) > 0, F 00 (s) ≥ 0 (c = a), (d) F (a) > 0, F (b) < 0, F 00 (s) ≤ 0 (c = a). Pˇr´ıpady (b) a (d) se snadno redukuj´ı na (a) a (c) volbou −F na m´ıstˇe F . (Tato u ´prava nemˇen´ı dokonce ani prvky posloupnosti {xn }.) Pˇr´ıpad (c) se pˇrevede na (a) volbou −x → x. (V tomto pˇr´ıpadˇe obdrˇz´ıme m´ısto {xn } → {−xn }. Staˇc´ı se tedy zab´ yvati pˇr´ıpadem (a). M´ame pak co do ˇcinˇen´ı s grafem rostouc´ı konk´avn´ı funkce F takov´e,ˇze F (a) < 0 < F (b). ˇx Bud ˆ jedin´ y koˇren rovnice F (s) = 0. Napˇred pˇredpokl´adejme, ˇze a ≤ x0 ≤ x ˆ. Potom, na z´akladˇe pˇredpokladu, ˇze F 0 (x0 ) ≥ 0, x1 = x0 −

F (x0 ) ≥ x0 . F 0 (x0 )

Indukc´ı lze uk´azati, ˇze (3.21)

ˆ. x0 ≤ x1 ≤ ... ≤ xn ≤ ... ≤ x

22

Indukˇcn´ı krok je realizov´ an takto. Z vˇety o stˇredn´ı hodnotˇe, −F (xn ) = F (ˆ x) − F (xn ) = F 0 (ξ)(ˆ x − xn ) a pˇritom xn ≤ ξn ≤ x ˆ. Z podm´ınky F 00 (s) ≤ 0 plyne, ˇze F 0 je nerostouc´ı a tud´iˇz 0 0 F (ξn ) ≤ F (xn ), takˇze −F (xn ) ≤ F 0 (xn )(ˆ x − xn ) a tedy xn+1 = xn −

F (xn ) ≤ xn + (ˆ x − xn ) = x ˆ. F 0 (xn )

Aplikac´ı F k obˇema stran´am nerovnosti obdrˇz´ıme (F je neklesaj´ıc´ı), ˇze F (xn+1 ) ≤ 0 a tud´ıˇz xn+2 = xn+1 −

F (xn+1 ) ≥ xn+1 . F 0 (xn+1 )

Tedy plat´ı (3.21). Na tomto m´ıstˇe aplikujeme jednu z fundament´ aln´ıch vˇet matematick´e anal´ yzy: Kaˇzd´ a ohraniˇcen´ a monotonn´ı posloupnost re´ aln´ych ˇc´ısel je konvergentn´ı. Oznaˇcme tedy x ˜ = lim xn . n→∞

Z (3.21) plyne, ˇze x ˜≤x ˆ. Na druh´e stranˇe z definice posloupnosti {xn } obdrˇz´ıme pˇri ∞ rovnost F (˜ x) x ˜=x ˜− 0 F (˜ x) x) = 0. Z jednoznaˇcnosti pak x ˜=x ˆ. T´ım je dok´az´ a tedy F (˜ ano tvrzen´ı vˇety 3.4 pro pˇr´ıpad x0 ≤ x ˆ. D´ale vyˇsetˇrujme p´r´ıpad x0 ≥ x ˆ. Opˇet pouˇzit´ım vˇety o stˇredn´ı hodnotˇe F (x0 ) = F (x0 ) − F (ˆ x) = F 0 (ξ0 )(x0 − x ˆ), kde x ˆ ≤ ξ0 ≤ x0 a d´ıky monotonii F 0 , F (x0 ) ≥ F 0 (x0 )(x0 − x ˆ). Snadno zjist´ıme, ˇze F (x0 ) ≤ x0 − (x0 − x ˆ) = x ˆ x1 = x0 − 0 F (x0 ) a tud´ıˇz F (x0 ) ≤ F (b) − (b − x0 )F 0 (b). Na z´akladˇe podm´ınky (iv) obdrˇz´ıme, ˇze plat´ı x1 = x0 −

F (x0 ) F (x0 ) ≥ x0 − 0 F 0 (x0 ) F (b) 23

≥ x0 −

F (b) + b − x0 F 0 (b)

= x0 − (b − a) + b − x0 = a. Tud´ıˇz a ≤ x1 ≤ x ˆ. T´ım jsme se dostali do situace, jiˇz jsme vyˇsetˇrovali v pˇredchoz´ı ˇc´asti d˚ ukazu s t´ım rozd´ılem, ˇze na m´ıstˇe x0 vystupuje x1 . Podle jiˇz dok´azan´eho, {xn } konverguje k x ˆ. T´ım je d˚ ukaz vˇety 3.4 proveden. ||| Vraˇtme se jeˇstˇe k pˇr´ıkladu 3.2. O procesu v nˇem vyˇsetˇrovan´em jiˇz v´ıme, ˇze je konvergentn´ ı pro poˇc´ateˇcn´ı pˇribl´ıˇzen´ı dostateˇcnˇe bl´ızk´ a k pˇresn´emu ˇreˇsen´ı, t. √ j. k c, kde c > 0 je dan´e kladn´e ˇc´ıslo, jehoˇz odmocninu hled´ame. Poloˇzme tedy F (s) = s2 − c a f (s) = s −

F (s) . F 0 (s)

Snadno zjist´ıme, ˇze F 0 (s) = 2s > 0, a F 00 (s) = 2 > 0 a tud´ıˇz nast´av´a v tomto pˇr´ıpadˇe situace (c) popsan´a v d˚ ukazu vˇety 3.4. V √ kaˇzd´em intervalu (a, b), kde 0 < a < c < b, plat´ı, ˇze Œ Œ Œ F (a) Œ c − a2 Œ Œ Œ F 0 (a) Œ = 2a ≤ b − a, pro kaˇzd´e b splˇ nuj´ıc´ı nerovnost

b≥

1 (a + c/a). 2

Z vˇety 3.4 plyne konvergence posloupnosti ’ “ √ c 1 xn+1 = xn + → c 2 xn pro libovoln´e x0 > 0. yhodn´e je pouˇzit´ı metody Newtonovy k hled´an´ı nulov´ ych bod˚ u polyZvl´aˇsˇt v´ nom˚ u. V tom pˇr´ıpadˇe vyuˇzijeme t´e skuteˇcnosti, ˇze um´ıme efektivnˇe stanoviti jak hodnoty p(z) tak p0 (z), kde p je dan´ y mnohoˇclen a z je dan´a veliˇcina.

24

ˇ tedy Bud p(z) =

N X

ak z k ,

k=0

pak pomoc´ı vztah˚ u (3.22) a (3.23)

b0 = a0 , bn = zbn−1 + an , n = 1, 2, ..., N ; c0 = b0 , cn = zcn−1 + bn , n = 1, 2, ..., N − 1,

obdrˇz´ıme rovnosti

bN = p(z) a cn−1 = p0 (z). Sch´ematicky postupujeme jak zn´azornˇeno n´ıˇze a2 a0 a1 ... ↓ ↓ ↓ ... b0 ⇒ b1 ⇒ b2 ⇒ ... ↓ c0 → c1 → c2 → ...

aN −1 aN ↓ ↓ bN −1 ⇒ bN cN −1

pˇri ˇcemˇz → znamen´a sˇc´ıt´an´ı a ⇒ n´asoben´ı faktorem z a sˇc´ıt´ an´ı. ˇ nyn´ı z koˇren polynomu p. Pak Bud p(x) = (x − z)q(x), kde q je polynom stupnˇe N − 1, q(x) =

N −1 X

bk xN −k−1 , bN = 0.

k=0

Snadno zjist´ıme, ˇze

a0 = b0 a Tud´ıˇz

an−1 = −bn z + bn , n = 1, ..., N − 1. N X

k=0

pˇri ˇcemˇz neboli, (3.24)

ak xN −k =

N −1 X k=0

bk xN −k − z

−1 N X

bk xN −k−1 ,

k=0

an = bn − zbn−1 , bn = an + zbn−1 .

Vid´ıme, ˇze jsme obdrˇzeli vztahy (2.20) algoritmu 2.1. Plat´ı tedy 25

Vˇ eta 3.5 Je-li z koˇren polynomu p, kde p(x) =

N X

k=0

ak xN −k , a0 6= 0,

a koeficienty b0 , ..., bN −1 se poˇc´ıtaj´ı podle formule (3.24), pak q(x) = (x − z)−1 p(x) =

N −1 X

bk xN −k , b0 = a0 .

k=0

Mohlo by se zd´at, ˇze je v´ yhodn´e hledat dalˇs´ı koˇreny polynomu p jakoˇzto a procedura nulov´e body polynomu niˇzˇs´ıho stupnˇe q. Ukazuje se vˇsak, ˇze takov´ je numericky nestabiln´ı a je tud´ıˇz vhodnˇejˇs´ı poˇc´ıtat vˇsechny koˇreny polynomu pomoc´ı p˚ uvodn´ıho polynomu p a nikoliv modifikovan´eho polynomu q. usobem obdobn´ym anal´yze Newtonovy metody metodu Cviˇ cen´ı 3.1 Analyzujte zp˚ n´ asleduj´ıc´ı, metodu regula falsi, definovanou pˇredpisem xn+1 = xn −

(xn − xn−1 )F (xn ) , F (xn ) − F (xn−1 )

kde F je funkce, jej´ıˇz koˇreny hled´ ame, t. j. F (x) = 0.

Cviˇ cen´ı 3.2 Jak´e vlastnosti m´ a Newtonova metoda v pˇr´ıpadˇe, kdy F 0 (x∗ ) = 0, kde F (x∗ ) = 0. Cviˇ cen´ı 3.3 Analyzujte modifikovanou Newtonovu sonovu) metodu danou pˇredpisem

(Whittakerovu - Robin-

1 F (xn ), m kde m je pevnˇe zvolen´ a konstanta, na pˇr. m = F 0 (x0 ). xn+1 = xn −

4 4.1

Iteraˇ cn´ı metody ˇ reˇ sen´ı soustav neline´ arn´ıch rovnic Vˇ eta o kontrakci

Pˇredmˇetem naˇseho studia bude vyˇsetˇrov´ an´ı n´asleduj´ıc´ı u ´lohy. Nal´ezt prvek x? ∈ RN takov´ y, ˇze   x1 = f1 (x1 , ..., xN ),     x2 = f2 (x1 , ..., xN ),   . (4.1) .     .    xN = fN (x1 , ..., xN ), 26

neboli, ps´ano vektorovˇe, (4.2) T

kde x = (x1 , ..., xN )

x = f (x),

a f = (f1 , ..., fN )T .

Pˇ r´ıklad 4.1 Nechˇt N = 2 a f1 (x1 , x2 ) = x21 + x22 , f2 (x1 , x2 ) = x21 − x22 .Potom m´a (4.1) tvar x1 = x21 + x22 , x2 = x21 − x22 .

V´ yraz x21 + x22 − x1 = 0 je rovnic´ı kruˇznice se stˇredem v bodˇe (1/2, 0) a − x22 − x2 = 0 je rovnic´ı hyperboly se stˇredem v bodˇe (0, −1/2). Obˇe tyto kˇrivky proch´azej´ı poˇc´atkem, coˇz znaˇc´ı, ˇze jedn´ım z nulov´ ych bod˚ u je x1 = 0, x2 = 0. Pomoc´ı graf˚ u pro (x1 − 1/2)2 + x22 = 1/4 a (x2 + 1/2)2 − x21 = 1/4 ˜2 = 0, 4 leˇz´ı dalˇs´ı (netrivi´aln´ı) snadno nahl´edneme, ˇze v okol´ı bodu x ˜1 = 0, 8, x nulov´ y bod soustavy (4.1). x21

Pˇri anal´ yze d˚ ukaz˚ u konvergence metod studovan´ ych v pˇredchoz´ı kapitole zjist´ıme, ˇze v´ yznamn´ ym d˚ ukazov´ ym prostˇredkem byly vlastnosti absolutn´ı hodnoty re´aln´ ych ˇc´ısel: |x| = 0 ⇔ x = 0; a pro c, x a y ∈ R.

|cx| = |c||x|;

|x + y| ≤ |x| + |y|

V souvislosti s prostory RN dimenze N ≥ 2, to vede k nutnosti zobecnit pojem absolutn´ı hodnoty. Matematika na pˇrelomu 19. a 20. stolet´ı dospˇela pˇri tomto zobecˇ nov´an´ı k pojmu normy a normovan´eho prostoru. arn´ı prostor nad tˇelesem re´ aln´ych ˇc´ısel. Funkce ν : Definice 4.1 Budˇ E line´ 1 E → R+ se naz´yv´ a normou na E, jestliˇze plat´ı n´ asleduj´ıc´ı relace: (i) ν(x) = 0 ⇔ x = 0; (ii) ν(cx) = |c|ν(x) x ∈ E, c ∈ R1 ; (iii) ν(x + y) ≤ ν(x) + ν(y), x, y ∈ E. Existuje-li na E norma, tak se prostor E se naz´yv´ a normovan´ ym. Pozn´ amka. Obvykle se norma na E oznaˇcuje symbolem || . || . Je-li zapotˇreb´ı zv´ yraznit, ˇze norma je br´ ana na prostoru E, pak se znaˇc´ı n´asledovnˇe || . ||E . a konvergentn´ı k x ∈ E a x Definice 4.2 Posloupnost {xn }, xn ∈ E se naz´yv´ jej´ı limitou, jestliˇze plat´ı lim ||xn − x|| = 0. n→∞

27

a cauchyovskou, jestliˇze pro Definice 4.3 Posloupnost {xn }, xn ∈ E, se naz´yv´ kaˇzd´e ε > 0 existuje N0 takov´e, ˇze ||xn+p − xn || < ε pro n ≥ N0 a p ≥ 1. Normovan´y prostor E se naz´yv´ a Banachov´ ym prostorem je-li v nˇem kaˇzd´ a cauchyovsk´ a posloupnost konvergentn´ı. Pˇ r´ıklad 4.2 Nechˇt E = RN a (a)   ||x||p =

N X

k=1

(b)

|xk |

!1/p

, 1 ≤ p < +∞,

||x||∞ = max {|xk | : k = 1, ..., N } .

Vˇs´ımnˇeme si, ˇze pro N = 1

||x||∞ = |x|, takˇze absolutn´ı hodnota v R1 je norma. D´ale nechˇt C = (cjk ) je regul´arn´ı matice typu N xN . Potom v´ yraz ||x||C = ||Cx||RN je norma na RN . Speci´alnˇe nechˇt E = (RN , ||.||2 ) a C = (cjk ), detC 6= 0. D´ale nechˇt N X (x, y) = xk yk . k=1

T

Poloˇzme B = C C a

||x||22,B = (Bx, x).

Snadno se uk´aˇze, ˇze ||.||2,B je norma na RN .

Vraˇtme se vˇsak k vyˇsetˇrov´ an´ı soustav neline´arn´ıch rovnic v RN . Metoda postupn´ ych aproximac´ı je definov´ ana podobnˇe jako pro pˇr´ıpad jedn´e rovnice. u {xn } Algoritmus 4.1 Zvolme prvek x ∈ E a poˇc´ıtejme posloupnost vektor˚ rekurz´ıvnˇe pomoc´ı formule (4.3)

xn+1 = f (xn ), n = 1, ...,

Opˇet si klademe ot´azku, zda je proces (4.3) konvergentn´ı a jak´e vlastnosti m´a pˇr´ıpadn´a limita; zˇrejmˇe se nab´ız´ı ˇreˇsen´ı rovnice (4.2). Odpovˇedi na poloˇzen´e ot´azky jsou obsahem n´asleduj´ıc´ıho tvrzen´ı. 28

QN Vˇ eta 4.1 Nechˇt I = j=1 Ij , kde Ij = [aj , bj ], j = 1, ..., N , a nechˇt funkce f1 , ... fN vyhovuj´ı n´ asleduj´ıc´ım podm´ınk´ am: (i) f1 , ..., fN jsou definov´ any a jsou spojit´e na I; (ii) Pro kaˇzd´y prvek x ∈ I t´eˇz vektor f (x) ∈ I; (iii) Existuje konstanta L, 0 < L < 1, takov´ a, ˇze pro libovoln´ a x1 , x2 ∈ I plat´ı (4.4) ||f (x1 ) − f (x2 )|| ≤ L||x1 − x2 ||. Potom plat´ı n´ asleduj´ıc´ı v´yroky (a) Rovnice (4.2) m´ a pr´ avˇe jedno ˇrevsen´ı x∗ ∈ I. a v (b) Pro libovoln´y v´ychoz´ı prvek x0 ∈ I posloupnost {xn } urˇcen´ algoritmu 4.1 je definov´ ana pro kaˇzd´e n a konverguje k x∗ . (c) Pro libovoln´e n ≥ 1 plat´ı odhad (4.5)

||xn − x∗ || ≤

Ln ||x1 − x0 ||. 1−L

Definice 4.4 Konstanta L, poˇzadovan´ a v (4.4) se naz´yv´ a Lipschitzovou konstantou (vzhledem k normˇe ||.||) vektorov´e funkce f. Vztah (4.4) se naz´yv´ a Lipschitzovou podm´ınkou pro f vzhledem k normˇe ||.||. Pozn´ amky. Podm´ınky kladen´e ve vˇetˇe 4.1 na Lipschitzovu konstntu L napov´ıdaj´ı interpretaci metody. Vzd´alenost obraz˚ u je pˇri takov´em zobrazen´ı u, coˇz d´av´ a takov´emu zobrazen´ı jm´eno : f ostˇre majorizov´ana vzd´alenost´ı vzor˚ kontraktivn´ı zobrazen´ı, ˇcili kontrakce. Je na m´ıstˇe podotknout, ˇze tvrzen´ı (a) je v pˇr´ıpadˇe N ≥ 2 daleko hlubˇs´ım v´ yrokem neˇz je tomu pro N = 1, kdy existenci ˇreˇsen´ı lze nahl´ednout na pˇr. vyˇsetˇren´ım grafu zkouman´e funkce. Tvrzen´ı (b) vypov´ıd´ a o konvergenci algoritmu, zat´ımco v (c) je d´an horn´ı odhad chyby. ana D˚ ukaz vˇ ety 4.1. Z poˇzadavku (ii) plyne, ˇze posloupnost {xn } je definov´ pro kaˇzd´e n ≥ 1 a x0 ∈ I. Na z´akladˇe (iii) odvod´ıme snadno, ˇze (4.6)

||xn+1 − xn || ≤ Ln ||x1 − x0 ||, n = 0, 1, ...

Fixujme n pevnˇe a vezmˇeme m > n, m = n + p, p ≥ 1 . Vid´ıme, ˇze xm − xn =

p X

k=1

a tedy, ||xm − xn || ≤

(xn+k − xn+k−1 ) ,

p X

k=1

||xn+k − xn+k−1 | |.

29

Pouˇzit´ım (4.6) obdrˇz´ıme vztahy ||xm − xn || ≤

(4.7)

Pp

k=1

≤

Ln+k−1 ||x1 − x0 | |

Ln 1−L ||x1

− x0 ||.

Protoˇze E = RN je u ´pln´ y (Banach˚ uv) prostor, plyne odtud konvergence posloupnosti {xn }. Nechˇt ˆ = lim xn , x n→∞

coˇz je tot´eˇz co ˆ || = 0. lim ||xn − x

n→∞

Vzhledem ke spojitosti funkc´ı f1 , ..., fN obdrˇz´ıme z (4.7) lim f (xn ) = f (ˆ x)

n→∞

a tud´ıˇz ˆ = lim xn+1 = lim f (xn ) = f (ˆ x x), n→∞

n→∞

ˆ = x∗ je ˇreˇsen´ım soustavy rovnic (4.2). a tak x Jednoznaˇcnost se dokazuje zcela analogicky jako v pˇr´ıpadˇe N = 1. Nechˇt x1∗ a x∗2 jsou dvˇe ˇreˇsen´ı soustavy (4.2). Potom, na z´akladˇe (4.4) ||x∗1 − x2∗ || = ||f (x∗1 ) − f (x2 )|| ≤ L||x1∗ − x∗2 ||, coˇz je moˇzn´e pouze kdyˇz x∗1 = x2∗ . Platnost (4.5) plyne z (4.7) pro p → ∞ : ||xn − x∗ || = lim ||xn − xn+p || ≤ n→∞

T´ım je d˚ ukaz vˇety 4.1 proveden. |||

4.2

Ln ||x1 − x0 ||. 1−L

Kvadratick´ a konvergence a Newtonova metoda pro soustavy

Zab´ yvejme se ot´azkou explicitn´ıho vyj´adˇren´ı vektoru chyby dn = xn − x∗ pˇri n → ∞. nuje pˇredpoklady vˇety 4.1 a nav´ıc, ˇze funkce f1 , ..., fN Pˇredpokl´adejme, ˇze f splˇ maj´ı v I spojit´e parc´aln´ı derivace aˇz do ˇr´ adu 2. vˇcetnˇe. Aplikac´ı Taylorovy vˇety pro funkce N promˇenn´ ych zjist´ıme, ˇze dn+1 = xn+1 −x∗ = f (xn )−f (x∗ ) = f (x∗ +dn )−f (xn ) = J (x∗ )dn +O(||dn ||2 ), (4.8) 30

kde (4.9)



 (f1 )x1 (x1 , ..., xN ) ... (f1 )xN (x1 , ..., xN )  ... ... ... J (x) =  (fN )x1 (x1 , ..., xN ) ... (fN )xN (x1 , ..., xN )

a O(||d2n ||) znaˇc´ı v´ yrazy maj´ıc´ı majorantu tvaru c ||dn ||2 , pˇri ˇcemˇz c je konstanta nez´avisl´a jak na n tak na x ∈ I. Definice 4.5 Matice definovan´ a v (4.9) se naz´yv´ a Jacobiovou matic´ı soustavy (4.1) resp. (4.2).

aln´ım analogem vztahu Vˇs´ımnˇeme si toho, ˇze vztah (4.8) je N -dimenzion´ (3.12). Je-li J (x∗ ) nenulov´a matice, pak (4.8) ˇr´ık´ a, ˇze vektor chyby se v kaˇzd´em iteraˇcn´ım kroku n´asob´ı Jacobiovou matic´ı (4.9), coˇz odpov´ıd´ a line´ arn´ı konvergenci. Je-li J (x∗ ) = 0, je z (4.8) patrn´e, ˇze chyba po 1+n- t´em kroku je u ´mˇern´ a ˇctverci chyby po n-t´em kroku - v tom pˇr´ıpadˇe m´ame co do ˇcinˇen´ı s kvadratickou konvergenc´ı. Opˇet m˚ uˇzeme formulovat analog vˇety 3.3 pro N > 1. Vˇ eta 4.2 Nechˇt funkce f1 , ..., fN definovan´e na I splˇ nuj´ı na I tyto poˇzadavky: (i) Existuj´ı spojit´e derivace ∂fj , 1 ≤ j , l ≤ N; ∂xl (ii) Soustava (4.1) m´ a uvnitˇr I ˇreˇsen´ı x∗ takov´e, ˇze J (x∗ ) = 0 . Potom existuje ˇc´ıslo δ > 0 takov´e, ˇze posloupnost definovan´ a pomoc´ı algoritmu 4.1 konverguje k x∗ pro libovoln´y startovac´ı vektor x0 , pro nˇejˇz ||x0 − x∗ || < δ. ychoz´ı pˇribl´ıˇzen´ı dostateˇcnˇe bl´ızk´e pˇresn´emu Tedy opˇet, jako pro N = 1, je-li v´ ˇreˇsen´ı, doch´az´ı ke konvergenci. Glob´aln´ı konvergence je pro v´ıcerozmˇern´ y pˇr´ıpad probl´em principi´alnˇe komplikovanˇejˇs´ı neˇz pro N = 1 (konvexita ev. konkavita ˇreˇsen´ı se jev´ı pro N = 1 dost akademicky). Podm´ınky zaruˇcuj´ıc´ı glob´aln´ı konvergenci jsou dost sloˇzit´e a jejich anal´ yza patˇr´ı mezi speci´aln´ı partie numerick´ ych metod. Hledejme ˇreˇsen´ı soustavy (4.10) F(x) = 0, kde F = (F1 , ..., FN )T .

31

Algoritmus 4.2 (Newtonova metoda pro soustavy) Volme x0 ∈ I a sestrojme posloupnost {xn } podle schematu xn+1 = xn − hn ,

(4.11)

kde hn = (h01 , ..., h0N )T je ˇreˇsen´ım soustavy (4.12)

J (xn )hn = −F(xn )

a



J =

∂F1 ∂x1

... ... ...

.

∂FN ∂xN

∂F1 ∂xN

.

∂FN ∂xN

.

 

Jak je patrn´e z (4.11) pˇri prov´ adˇen´ı algoritmu 4.2 mus´ıme pˇri kaˇzd´em iteraˇcn´ım kroku ˇreˇsiti soustavu line´arn´ıch algebraick´ ych rovnic (4.12). Jest proto pˇrirozen´e poˇzadovati (4.13) detJ (xn ) 6= 0, n = 0, 1, ... Cviˇ cen´ı 4.1 Naleznˇete matici H = H(x) tak, aby uˇzit´ı postupn´ych aproximac´ı na rovnici x = x − H(x)F (x) poskytlo kvadratickou konvergenci. Rozhodnˇete, zda se vˇzdy obdrˇz´ı Newtonova metoda.

5 5.1

ˇ sen´ı soustav line´ Reˇ arn´ıch algebraick´ ych rovnic Obecn´ e soustavy a ˇ reˇ sen´ı ve smyslu nejmenˇ s´ıch ˇ ctverc˚ u

Budeme se zab´ yvat ˇreˇsen´ım soustav typu (5.1)

Ax = b,

kde x ∈ RN , b ∈ RM , A = (ajk ), 1 ≤ j ≤ M, 1 ≤ k ≤ N.

Pozn´ amka. Je zˇrejm´e, ˇze lze vˇzdy doc´ılit situace, kdy M = N a to rozˇs´ıˇren´ım matice A o vhodn´ y poˇcet nulov´ ych sloupc˚ u pro pˇr´ıpad N < M ˇci ˇr´ adk˚ u pro pˇr´ıpad N > M . ˇ sitelnost soustavy v klasick´em smyslu t. j. kdy vztah (5.1) je ch´ Reˇ ap´ an jako soustava rovnost´ı N X ajk xk = bj , j = 1, ..., M, k=1

32

nelze vyˇzadovat. Jako pˇr´ıklad m˚ uˇze slouˇzit ˇcast´ y pˇr´ıpad, kdy soustava (5.1) representuje nˇejak´ y model v nˇemˇz nˇekter´e z prvk˚ u ajk a bj se obdrˇz´ı mˇeˇren´ım. Obecnˇe nelze vylouˇcit, aby odpov´ıdaj´ıc´ı soustava pro pro dvˇe r˚ uzn´ a mˇeˇren´ı byla v jednom p´r´ıpadˇe ˇreˇsiteln´a a v druh´em nikoliv. V praktick´ ych aplikac´ıch se ukazuje potˇreba n´asleduj´ıc´ıho tvrzen´ı. Kaˇzd´ a soustava typu (5.1) m´ a pr´ avˇe jedno ”ˇreˇsen´ı”. Takov´e tvrzen´ı nem˚ uˇze platit pro klasick´e ˇreˇsen´ı (podm´inka o hodnostech). Aby poˇzadovan´e tvrzen´ı platilo, mus´ıme vhodn´ ym zp˚ usobem zobecnit pojem ˇreˇsen´ı. Proto zavedeme zobecnˇen´e ˇreˇsen´ı. ˇ (. , .)M a (. , .)N skal´ Budte arn´ı souˇciny na RM a RN . Definujme (5.2) a hledejme x ∈ R (5.3)

N

F (x, x) = (Ax − b, Ax − b)M

tak, aby

F (x) = min{F (x) : x ∈ RN }.

n jeden vektor x ∈ RN takov´y, ˇze plat´ı (5.3). Vˇ eta 5.1 Existuje alespoˇ D˚ ukaz. Snadno ovˇeˇr´ıme, ˇze plat´ı relace (5.4)

RM = rangeA ⊕ kerA∗

(5.5)

RN = rangeA∗ ⊕ kerA,

kde A∗ je definov´ana pomoc´ı vztahu A∗ v = y, pˇri ˇcemˇz (Ax, v)M = (x, y)N pro libovoln´a x ∈ RN a v ∈ RM . Existence i jednoznaˇcnost y je d˚ usledkem Rieszovy vˇety o reprezentaci spojit´eho line´arn´ıho funkcion´alu na Hilbertovˇe prostoru (RN , (. , .)N ). Na z´akladˇe (5.4) a (5.5) obdrˇz´ıme, ˇze plat´ı x = A∗ u + v, Av = 0, a b = Ac + d, A∗ d = 0. Odtud pak F (x) = (Au − Ac, Au − Ac)M + (d, d)M . 33

Vzhledem k tomu, ˇze d je pevn´ y vektor, staˇc´ı tedy poloˇziti ¯ = u = c, x takˇze (5.6)

F (¯ x) = (d, d)M = min{F (x) : x ∈ RN }.

T´ım je vˇeta 5.1 dok´az´ana.

Pozn´ amka. Protoˇze (5.6) je hodnota glob´aln´ıho minima kvadratick´e funkce F , je veliˇcina (5.6) t´aˇz pro vˇsechna eventu´ aln´ı zobecnˇen´ a ˇreˇsen´ı x soustavy Ax = yv´a cenou uveden´e u ´lohy naj´ıt zobecnˇen´e ˇreˇsen´ı, neboli b. Tato hodnota se naz´ cenou zobecnˇen´eho ˇreˇsen´ı. Vˇs´ımnˇeme si t´e skuteˇcnosti, ˇze cena zobecnˇen´eho ˇreˇsden´ı x je nulov´a pr´avˇe kdyˇz x je ˇreˇsen´ı klasick´e. Nechˇt M = {x ∈ RN : F (x) = (d, d)M }.

Z vˇety 5.1 plyne, ˇze M 6= ∅ a ze spojitosti F pak, ˇze M je uzavˇren´ a. D´ıky subaditivitˇe normy snadno nahl´edneme, ˇze M je konvexn´ı. Existuje tedy jedin´ y prvek x∗ ∈ M takov´ y, ˇze ||x∗ ||N = min{||x||N : x ∈ M}. Zobecnˇen´e ˇreˇsen´ı x∗ se naz´ yv´ a norm´ aln´ım ˇreˇsen´ım soustavy (5.1). Vˇ eta 5.2 Existuje pr´ avˇe jedno norm´ aln´ı ˇreˇsen´ı libovoln´e soustavy Ax = b. Vid´ıme tedy, ˇze n´ami poˇzadovan´e tvrzen´ı uveden´e na poˇc´ atku t´eto kapitoly, plat´ı pro n´ami zaveden´e zobecnˇen´e ˇreˇsen´ı. Pozn´ amka. Zobecnˇen´ a ˇreˇsen´ı a tedy t´eˇz norm´aln´ı ˇreˇsen´ı jsou podstatnˇe z´avisl´a na skal´arn´ıch souˇcinech (. , .)M a (. , .)N . Vol´ıme-li speci´alnˇe (5.7)

(u, v)M =

M X

(u)j (v)j ,

j=1

kde (u)j znaˇc´ı j-tou komponentu vektoru u ∈ RM , obdrˇz´ıme funkci F ve tvaru F (x) =

"N M X X j=1

k=1

ajk (x)k − (b)j

#2

.

Odtud z´ıskalo zobecnˇen´e ˇreˇsen´ı takto zaveden´e n´azev ˇreˇsen´ı ve smyslu nejˇ sen´ı ve smyslu nejmenˇs´ıch ˇctverc˚ menˇs´ıch ˇctverc˚ u. Reˇ u patˇr´ı mezi nejrozˇs´ıˇrenˇejˇs´ı. Toto zobecnˇen´e ˇreˇsen´ı bylo zavedeno a velmi zevrubnˇe aplikov´ ano jiˇz Gaussem. Zobecnˇ en´ e inverzn´ı oper´ atory 34

ˇ A matice typu M × N . Bud Vyˇsetˇrujme soustavu relac´ı (i)AXA = A,

(ii)XAX = X,

(iii) (AX)∗ = AX,

(iv) (XA)∗ = XA.

Pˇ r´ıklad. Nechˇt M = N a detA 6= 0. Poloˇzme X = A−1 . Splnˇen´ı vztah˚ u (i) - (iv) s X = A−1 je oˇcividn´e. Napˇred si uk´aˇzeme, ˇze matice X typu N × N splˇ nuj´ıc´ı vztahy (i) - (iv) existuje nejv´ yˇse jedna. ˇ tedy X1 a X2 dvˇe takov´e matice. Jest potom Budte X1 = X1 AX1 = X1 (AX1 )∗ = X1 (AX2 AX1 )∗ = X1 (AX1 )∗ (AX2 )∗ = X1 AX1 AX2 = X1 AX2 = (X1 A)∗ X2 AX2 = (X1 A)∗ (X2 A)∗ X2 . Na druh´e stranˇe, (X2 AX1 A)∗ X2 = X1 AX2 AX2 = (X1 A)∗ (X2 A)∗ X2 , takˇze X1 = (X1 AX2 A)∗ X2 = X2 AX2 = X2 . Jiˇz v´ıme, ˇze v pˇr´ıpadˇe M = N a detA 6= 0 jedin´ ym ˇreˇsen´ım soustavy vztah˚ u (i) - (iv) je inverzn´ı matice A−1 . Existuje ˇreˇsen´ı v obecn´em pˇr´ıpadˇe? Ano! Poloˇzme T = A∗ A. Snadno nahl´edneme, ˇze matice T je typu N × N , je symetrick´ a a t´eˇz je pozitivnˇe semidefinitn´ı, t. j. (T x, x)N ≥ 0. Existuj´ı tedy projekce Pj , j = 1, ...s, s ≤ N takov´e, ˇze plat´ı vztahy (5.8)

T =

s X

λk Pk , T P0 = 0,

k=1

35

pˇri ˇcemˇz s X

(5.9) Pj Pk = Pk Pj = δjk Pk , Pj∗ = Pj ,

k=1

a d´ale pak (5.10) Poloˇzme

Pk = I − P0 , j, k = 0, 1, ..., N

(T − λj )Pj = 0, λj > 0, j = j, ..., N, λ0 = 0.

(5.11)

Z=

N X

∗ λ−1 j Pj A .

j=1

Vˇ eta 5.3 Ke kaˇzd´e matici A typu M × N existuje pr´ avˇe jedna matice typu N × M takov´ a, ˇze pro ni plat´ı vztahy (i) - (iv). a Definice 5.1 Matice, jej´ıˇz existenci a jednoznaˇcnost zaruˇcuje vˇeta 5.3, se naz´yv´ (Mooreovou - Penroseovou) pseudoinverzn´ı matic´ı a oznaˇcuje se symbolem A+ . ukaz jednoznaˇcnosti je element´ arn´ı a lze jej pˇrenechat ˇcten´ aˇri. D˚ ukaz. D˚ Na z´akladˇe (5.11) a (5.8) snadno zjist´ıme, ˇze AZA = A

N X

∗ λ−1 j Pj A A = A

N X N X

λj−1 λk Pj Pk .

j=1 k=1

j=1

Z tohoto vyj´adˇren´ı pomoc´ı (5.9) odvod´ıme, ˇze AZA = A(I − P0 ).

ˇ x ∈ RN . Podle pˇredpokladu (5.8) plat´ı, ˇze Bud 0 = (T P0 x, P0 x)N = (AP0 x, AP0 x)N , takˇze AP0 x = 0. Protoˇze x je libovoln´ y, plyne odtud, ˇze AP0 = 0. nuje vztah (i). Dok´azali jsme tak, ˇze , d´ıky (5.12), matice Z splˇ Platnost zb´ yvaj´ıc´ıch vztah˚ u (ii) - (iv) lze dok´azat analogicky, coˇz pˇrenech´ ame ˇcten´aˇri a pokl´ad´ame t´ım vˇetu 5.3 za dok´azanou s t´ım, ˇze A+ = Z, pˇri ˇcemˇz Z je definov´ana v (5.11). ||| Poloˇzme x+ = A+ b. (5.12) Dok´aˇzeme, ˇze plat´ı 36

ˇ A matice typu M ×N a b ∈ RM , b = b1 +b2 , b1 ∈ range(A), b2 ∈ Vˇ eta 5.4 Budte ker(A). Potom plat´ı rovnost x+ = x∗ . (5.13) D˚ ukaz. Zˇrejmˇe F (x∗ ) = (Ax+ − b, Ax+ − b)M = (AA+ b − b, AA+ b − b)M = ((I − P0 )b − b, (I − P0 )b − b)M = (b2 , b2 )M . Staˇc´ı tedy uk´azati, ˇze x+ ∈ R(A∗ ). To vˇsak plyne odtud, ˇze P0 x+ = 0 a tud´ıˇz, x+ = (I − P0 )x+ , ˇc´ımˇz je tvrzen´ı dok´az´ ano. ||| Vyˇsetˇrujme opˇet soustavu Ax = b, kde A je M × N matice a b ∈ RN a b = b1 + b2 , pˇri ˇcemˇz b1 ∈ range(A) t.j. b1 = Ac1 , d´ale pak a A∗ b2 = 0. Snadno zjist´ıme, ˇze plat´ı A∗ Ax = A∗ b = A∗ b1 , kter´aˇzto soustava se naz´ yv´a norm´ aln´ı soustavou. Tato soustava m´a vˇzdy klasick´e ˇreˇsen´ı. To vypl´ yv´a odtud, ˇze (A∗ A)∗ = A∗ A a pro kaˇzd´e ˇreˇsen´ı y homogenn´ı soustavy A∗ Ay = 0 plat´ı vztahy (y, A∗ b)N = (y, A∗ Ac1 )N ) = (A∗ Ay, c1 )N = 0. Hledat ˇreˇsen´ı norm´aln´ı soustavy rovnic je numericky velmi n´aroˇcn´e, protoˇze ˇc´ıslo podm´ınˇenosti smax κ(A∗ A) = , smin kde 1/2 smax = kA∗ Ak2

a smin je odmocnina nejmenˇs´ı kladn´e vlastn´ı hodnoty matice A∗ A, je zpravidla velmi velik´e; v pˇr´ıpadˇe ˇctvercov´e matice A je rovno ˇctverci ˇc´ısla podm´ınˇenosti nme, ˇze ˇc´ıslo podm´ınˇenosti ˇctvercoiv´e matice A je definov´ matice A. Pˇripomeˇ ano jakoˇzto v´ yraz kAkkA−1 k.

Norm´aln´ı ˇreˇsen´ı je tedy (klasick´e) ˇreˇsen´ı norm´aln´ı soustavy leˇz´ıc´ı v range(A+ ). D´ale plat´ı

Lemma 5.1 Kaˇzd´e zobecnˇen´e ˇreˇsen´ı x ∈ M m´ a tvar (5.14)

x = x+ + y,

kde Ay = 0. 37

ˇ x dalˇs´ı ˇreˇsen´ı soustavy Ax = b a nechˇt Proof. Bud x = a1 + a2 , a1 = A∗ c1 , Aa2 = 0. Potom plat´ı, ˇze

x − x = x+ − a1 + y − a2 .

Protoˇze Aa1 = b1 , odvod´ıme, ˇze plat´ı vztahy

A(x+ − a1 ) = AA+ b1 − b1 = 0. Jest tud´ıˇz kA+ b1 − A∗ c1 k2 = (A∗˜b1 − A+ c1 )N = (˜b1 − c1 , A[(x+ − a1 )M = 0, kde A+ b1 = A∗˜b1 a tedy

x = x = x+ .

ano. ||| Protoˇze x − x ∈ ker(A), tvrzen´ı lemmatu je dok´az´

5.2

Pomocn´ e prostˇ redky ke konstrukci speci´ aln´ıch reprezentac´ı matic

Nejprve si zopakujme nˇekter´e pojmy, jeˇz se uk´aˇz´ı jako vhodn´e pro konstrukci a vyˇsetˇrov´an´ı urˇcit´ ych algoritm˚ u. Pˇredpokl´ad´ame, ˇze vˇsechny matice v tomto ˇcl´ anku jsou obecnˇe komplexn´ı typu N × N . Matice A je hermiteovsky zdruˇzen´ a, jestliˇze AH = A, AH = (a∗jk ) pˇri ˇcemˇz a∗jk = a ¯kj , j, k = 1, ..., N. Je-li matice A re´aln´a, t. j. je-li a ¯jk = ajk naz´ yv´ame matici hermiteovskou symetrickou. Jest totiˇz AT = A, AT = (aTjk ), aTjk = akj , j, k = 1, ..., N. Permutaˇcn´ı matice P = (pjk ) splˇ nuje vztahy ( N N X X 0 , pjk = pjk = pkj = 1, 1 k=1

a

detP 6= 0. 38

k=1

Zˇrejmˇe P H = P T = P −1 . yv´ a matice splˇ nuj´ıc´ı Projekˇcn´ı matic´i ˇci projekc´ı se naz´ P 2 = P. Speci´alnˇe tedy, I 2 = I a 02 = 0. uz Permutaˇcn´ı matice, kter´a realizuje z´amˇenu l - t´e a k - t´e sloˇzky vektor˚ RN , se naz´ yv´a element´ arn´ı permutaˇcn´ı matic´ı neboli transpoziˇcn´ı matic´ı. Takov´a element´arn´ı permutaˇcn´ı matice m´a tvar k l ↓ ↓   1 . . . 0 . . . 0 . . . 0  . . .     . . .     . . .    k →  0 . . . 1 0   0   . . . . . .   .   0 l →  1 . . . 0 0    . . .     .  . .    . . .  0 . . . 0 . . . 0 . . . 1 a pˇri jej´ı aplikaci plat´ı vztahy y = P x, pˇri ˇcemˇz yj = Xj , j 6= k, l, zat´ımco yk = xl a yl = xk . Unit´ arn´ı matice U je definov´ a na vztahy U U H = U H U = I, coˇz pro re´alnou matici V znaˇc´ı, ˇze V V T = V T V = I. V takov´em pˇr´ı padˇe hovoˇr´ıme o matici ortogon´ aln´ı. Pro involutorn´ı matici B plat´ı B 2 = I.

39

Tedy speci´alnˇe I 2 = I. ˇ x ∈ RN y ∈ RM . Definujeme matici Budte y⊗v kladouce [y ⊗ v] = (x, v)N y.

Jin´ ymi slovy,

y ⊗ v = y vT .

arn´ı doln´ı troj´ uheln´ıkovou indexu k, jestliˇze Natice Ek se naz´ yv´a element´ E k = IN − v ⊗ e k a pˇritom ejT v = 0, j = 1, ..., k. To znaˇc´ı, ˇze 

       Ek =        



1 . . . 1 −vk+1 . . . −vN

1 . . .

Snadno provˇeˇr´ıme, ˇze y ⊗ v splˇ nuje n´asleduj´ıc´ı relaci

1

        .      

2

[y ⊗ v] x = (y, v)N [y ⊗ v] x. takˇze, je - li (y, v)N = 0, a tedy,

2

[y ⊗ v] = 0 σ(y ⊗ v) = {0}

neboli, y ⊗ v je nilpotentn´ı matice. V pˇr´ıpadˇe element´arn´ı matice Ek vid´ıme, ˇze (ek , v)N = 0, takˇze plat´ı Vˇ eta 5.5 Budˇ Ek = IN − vk ⊗ ek , pˇri ˇcemˇz (vk , ej )N = 0, j = k, k + 1, ..., N. Potom existuje inverzn´ı matice Ek−1 a plat´ı Ek−1 = IN + vk ⊗ ek . 40

Vˇ eta 5.6 Nechˇt xT = (ξ1 , ..., ξN ) a nechˇt 0 6= ξk = eTk x. Potom existuje jednoznaˇcnˇe urˇcen´ a element´ arn´ı doln´ı troj´ uheln´ıkov´ a matice indexu k Ek takov´ a, ˇze Ek x = (ξ1 , ..., ξk , 0, ..., 0)T . D˚ ukaz. Hledejme Ek ve tvaru Ek = IN − vk ⊗ ek . Z poˇzadavku elementarity Ek plyne, ˇze (5.15)

ν1 = ... = νk = 0,

kde vkT = (ν1 , ..., νN ). D´ale pak a protoˇze ˇz´ad´ame, aby (5.16)

Ek x = x − (x, ek )N v = x − ξk v ξj − ξk νj = 0, j = k + 1, ..., N,

z pˇredpokladu ξk 6= 0 odvod´ıme vztahy (5.17)

νj =

ξj , j = k + 1, ..., N. ξk

Tedy, Ek existuje a je vztahy (5.15), (5.16), (5.17) urˇcena jednoznaˇcnˇe. ||| Householderovu matici zavedeme pomoc´ı vztah˚ u H = I − 2wwH , wH w = 1, kde w je (obecnˇe komplexn´ı) sloupcov´ y vektor a wH odpov´ıdaj´ıc´ı vektor ˇr´ adkov´ y. Zˇrejmˇe Householderova matice je hermiteovsky zdruˇzen´ a H H = H a t´eˇz unit´arn´ı , neboˇt H H H = H 2 = (I − 2wwH )2 = I. Vˇs´ımnˇeme si, ˇze pro x ∈ C N a y = Hx plat´ı, ˇze y = x − 2wH xw, takˇze y H y = (y, y) = (Hx, Hx) = (H H Hx, x) = (x, x) a (x, y) = (y, x). 41

Lemma 5.2 K dan´emu vektoru v 6= 0 existuje ortogon´ aln´ı matice Q takov´ a, ˇze (5.18)

Qv = −σ||v||e1 ,

kde σ=

š

+1 pro v1 = (v, e1 ) ≥ 0 −1 pro v1 < 0.

D˚ ukaz. Nechˇt u=v+σ a

p (v, v)e1

Q=I −2 Zˇrejmˇe QH = I − 2 a

uuH u . uH u

uuH =Q uH u

QH Q = Q2 = I. D´ale pak Qv = Qu − σ||v||Qe1 .

Avˇsak Qu = −u a Qe1 = e1 − 2σ||v||u(u)1

1 . (u, u)

Odtud plyne, ˇze Qv = −u − σ||v oe1 + 2σ 2 ||v||2

(u)1 u, (u, u)

a protoˇze (u)1 = (v)1 + σ||v no, zjist´ıme, ˇze −u + 2σ||v|| = Tedy

(v)1 + σ||v|| (u, u)

‰ 1 ˆ −2σ 2 (v, v) − 2σ||v||(v)1 + 2σ 2 ||v||2 + 2σ||v||(v)1 = 0. (u, u) Qv = −σ||v||e1 . |||

42

Matice H = (hjk ) takov´a, ˇze hjk = 0 pro j − k > 2 a j ≥ 3, se naz´ yv´a horn´ı Hessenbergova. N´azornˇe je patrno rozm´ıstˇen´ı prvk˚ u na  × × × . .  × × × . .   0 × × . .   0 0 × . .   . . . . . 0 0 0 . .

n´asleduj´ıc´ım sch´ematu  . × × . × ×    . × ×  . . × ×  . . .  . × ×

Lemma 5.3 Ke kaˇzd´e matici A existuje unit´ arn´ı matice T , takov´ a, ˇze A = T HT H , kde matice H je horn´ı Hessenbergova. ˇ yv´ a slabˇe diCtvercov´ a matice A = (ajk ) ajk ∈ R j, k = 1, ..., N , se naz´ agon´ alnˇe dominantn´ı, jestliˇze plat´ı (5.19)

N X

k=1,k6=j

|ajk | ≤ |ajj |, j = 1, ..., N,

pˇri ˇcemˇz existuje alespoˇ n jeden index j0 , 1 ≤ j0 ≤ N , tak, ˇze (5.20)

N X

k=1,k6=j0

|aj0 k | < |aj0 j0 |.

Nastane-li v (5.19) ostr´a nerovnost pro vˇsechny indexy j = 1, ..., N , pak se matice A naz´ yv´a diagon´ alnˇe dominantn´ı. Vˇ eta 5.7 (Gerˇsgorinova) Budˇ A = (ajk ) matice N × N , pˇri ˇcemˇz ajk jsou komlexn´ı ˇc´ısla. D´ ale budˇ λ vlastn´ı hodnota matice A. Potom existuje index j0 , 1 ≤ j0 ≤ N takov´y, ˇze plat´ı vztahy (5.21)

|λ − aj0 j0 | ≤

N X

k=1,k6=j0

|aj0 k |.

ˇ x ∈ RN ⊕ i RN vlastn´ı vektor odpov´ıdaj´ıc´ı vlastn´ı hodnotˇe λ. D˚ ukaz. Bud Tedy, Ax = λx, x 6= 0. 43

Poloˇzme |xj0 | = max{|xj | : j = 1, ..., N }. Jest potom |λ − aj0 | ≤ ≤

PN

k=1,6=j0

PN

k=1,6=j0

|xk | |aj0 k | |x j | 0

|aj0 k |.

Je tak dok´az´ana platnost relace (5.21) a t´ım i vˇeta 5.7. ||| D˚ usledek 5.1 Kaˇzd´ a diagon´ alnˇe dominantn´ı matice A je regul´ arn´ı. D˚ ukaz. Staˇcı si uvˇedomit, ˇze ˇctvercov´ a matice je regul´arn´ı pr´avˇe kdyˇz 0 nen´ı jej´ı vlastn´ı hodnotou. To, ˇze 0 nen´ı vlastn´ı hodnotou dan´e ˇctvercov´e matice vˇsak je bezprostˇredn´ım d˚ usledkem diagon´aln´ı dominance. Kdyby totiˇz 0 byla vlastn´ı hodnotou matice A, platily by vztahy (5.21) pro λ = 0, |aj0 j0 | = |λ − aj0 j0 | ≤

N X

k=1,k6=j0

|aj0 k | < |aj0 j0 |.

Tento spor dokazuje platnost tvrzen´ı. ||| alnˇe dominantn´ı a budˇ E1 = IN − v1 ⊗ e1 Vˇ eta 5.8 Budˇ 0 6= A = (ajk ), diagon´ element´ arn´ı doln´ı troj´ uheln´ıkov´ a matice indexu 1, kde   0  a21 / a11      .   v1 =  . .     . aN 1 / a11 Potom

A(1) = E1 A

je t´eˇz diagon´ alnˇe dominantn´ı. D˚ ukaz.Vyˇsetˇrujme v´ y dle definice, yraz a1jk rovn´ a1jk = ajk −

aj1 a1k , 1, ..., N. a11

44

Jest tedy, 1 |− |ajj

PN

k=2, k6=j

|a1jk |

≥ |ajj | −

1 |a11 | |aj1 a1j |

= |ajj | −

PN

≥ |ajj | −

k=2, k6=j

PN

k=1, k6=j

−

PN

k=2, k6=j

|ajk | −

1 |a11 |

|ajk |

1 |a11 |

|ajk | −

PN

k=1, k6=j

PN

k=2, k6=j

|aj1 a1k |

|aj1 a1k |

≥ 0, j = 2, ..., N. ||| ˇ Ctvercov´ a matice A = (ajk ) nad tˇelesem komplexn´ıch ˇc´ısel se naz´ yv´ a pozitivnˇe definitn´ı, jestliˇze existuje konstanta κ > 0 takov´ a, ˇze (Ax, x)N ≥ κ(x, x)N rm pro vˇsechny vektoryx ∈ C. Jsou - li ajk ∈ R, pak nav´ıc poˇzadujeme, aby A = AT .

5.3

Pˇ r´ım´ e metody ˇ reˇ sen´ı soustav line´ arn´ıch algebraick´ ych rovnic

Pod pˇr´ım´ ymi metodami ˇreˇsen´ı soustav line´arn´ıch algebraick´ ych rovnic rozum´ıme takov´e metody, kter´e jsou zaloˇzeny na algoritmech, jeˇz poskytuj´ı pˇresn´e ˇreˇsen´ı po proveden´ı jist´eho koneˇcn´eho poˇctu aritmetick´ ych operac´ı. Je zˇrejm´e, ˇze pro soustavy, jejichˇz matice jsou troj´ uheln´ıkov´e - horn´ı ˇci doln´ı - se nab´ızej´ı pˇrirozen´e pˇr´ım´e metody. Jsou to t. zv. zpˇetn´ a a pˇr´ım´e dosazen´ı. Vyˇsetˇrujme soustavy (5.22) Ly = c, a (5.23) kde

a

U x = b, 

l11 . . .  l21 . . .  L= . . . . lN 1 . . . . 

u11  0  U = . 0

u12 u22 . 0

. . . . 45

0 0 . lN N



  

 . . u1N . . u2N  . . . .  . . uN N

Nechˇt ljj 6= 0, j = 1, ..., N

a

ujj = 6 0, j = 1, ..., N. Potom (5.24)

c1 l11

y1 =

a

k−1

yk =

(5.25) a podobnˇe

X lkj cj ck − , k = 2, ..., N, lkk j=1 ljj lkk

(5.26)

xN =

a (5.27)

xN −k =

bN −k

uN −k,N −k

−

N X

j=k+1

bN , uN N

uN −k,j bj , k = 1, ..., N − 1. uN −k,N −k ujj

Nejzn´amˇejˇs´ı a pˇritom nejd˚ uleˇzitˇejˇs´ı pˇr´ımou metodou ˇreˇsen´ı soustav lin´arn´ıch algebraick´ ych rovnic je Gaussova eliminaˇcn´ı metoda. Jej´ı z´akladn´ı myˇslenka uheln´ıkovou spoˇc´ıv´a v transformaci p˚ uvodn´ı soustavy s matic´ı A na soustavu s troj´ ˜ Matici A˜ lze z´ıskat postupn´ ym eliminov´ an´ım, jeˇz lze reprezentovat matic´ı A. element´arn´ımi eliminaˇcn´ımi maticemi popsan´ ymi v odstavci 5.2. ˇ tedy A = (ajk ) dan´a regul´arn´ı matice s a11 6= 0 a Bud E1 = I − v1 ⊗ e1 element´arn´ı doln´ı troj´ uheln´ıkov´ a matice indexu 1, kde   0  a12     1   . . v1 =   a11  .   .  aN 1 Jest tedy

a A(1)

A(0) = A 

a11  0 =  . 0

a12 a22 − aa12 a21 11 . aN 2 − aa12 aN 1 11 46

 . . . a1N  . . . a2N − aa1N a21  11 . . . . . a1N . . . aN N − a11 aN 1

Zvl´aˇstn´ı pozornosti si zasluhuje skuteˇcnost, ˇze prvn´ı ˇr´ adek transformovan´e matice A1 je totoˇzn´ y s prvn´ım ˇr´ adkem p˚ uvodn´ı matice A. ˇ 1 < k ≤ N − 1, ak−1 = 0, kde Obecnˇe, bud kk 6   (k−1) (k−1) a11 . . . . . . . a1N   .     .     .     (k−1) A(k−1) =  . . . a(kN )  , 0 . . . akk     . .     . .     . . (k−1) 0 . . . ak−1 . . . aN N Nk pˇri ˇcemˇz ˇr´adky s indexy 1 ≤ j ≤ k, jsou totoˇzn´e pro vˇsechny matice Al , l = j, ..., k − 1. Nechˇt E k = IN − vk ⊗ e k , kde



  vk =   

Potom 

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

(k−1)

(k−1)

ak+1,k /akk . . . (k−1) (k−1) aN k /akk



  .  

A(k) Ek A(k−1) = (k)

a11

0

0

(k)

(k)

. . . a1k A1k+1 . . . (k) akk . . . . . . 0

.

(k)

akk+1

(k)

aN k+1

.

.



a1N ř(k)

(k)

. .

.

akN

.

.

aN N

.

(k)

       ,      

pˇri ˇcemˇz ˇr´adky s indexy 1 ≤ j ≤ k+1 jso totoˇzn´e pro matice A(l) , kde l = j, ..., k. Pr´avˇe proveden´ y indukˇcn´ı krok ukazuje, ˇze ve ”fyziologick´em” pˇr´ıpadˇe, t. j. kdyˇz (k−1) akk 6= 0,

v´ yˇse uveden´ y proces vede po koneˇcn´em poˇctu krok˚ u, k ≤ N − 1, k situaci, kdy A(N −1) je horn´ı troj´ uheln´ıkov´ a matice. 47

Hledan´e ˇreˇsen´ı p˚ uvodn´ı soustavy je evidentnˇe totoˇzn´e s ˇreˇsen´ım soustavy A(N −1) x = b(N −1) , kde

b(N −1) = EN −1 ...E1 b.

Schematicky lze tuto skuteˇcnost postihnout tak, ˇze sestroj´ıme obd´eln´ıkovou matici typu N × N + 1 B = (A, b) a prov´ad´ıme postupnˇe transformace s touto matic´ı, tedy konstruujeme postupnˇe matice B (k) = Ek B (k−1) , k = 1, ..., , pˇri ˇcemˇz

B (0) = B.

Po N − 1 kroc´ıch obdrˇz´ıme  ×  0 B (N −1) =   . 0

matici B N −1 maj´ıc´ı tvar × . . . × . × . . . × . . . . . . . 0 . . . ×. .

. . . . . . . .

 ×  ×  . .  ×

Je vcelku zˇrejm´e, jak uveden´ y postup zobecnit na pˇr´ıpad souˇcasn´eho ˇreˇsen´ı soustav s danou matic´ı soustavy a s v´ıce prav´ ymi stranami b1 , ..., bp , p ≥ 1. V takov´em pˇr´ıpadˇe B(p) = (A, b1 , ..., bp ) s

(k)

(k−1)

B(p) = Ek B(p) kde

,

(0)

B(p) B(p) . Jiˇz v´ıme, ˇze ˇreˇsen´ı soustav Ax = Bj , j = 1, ..., p, se pˇrev´ad´ı na ˇreˇsen´ı spustav (N −1)

A(N −1) x = bj

, j = 1, ...p,

(N −1)

y sloupec matice B(p) . Matice tˇechto soustav jsou vˇsak kde bj je N + j - t´ horn´ı troj´ uheln´ıkov´e a nalezen´ı hledan´ ych ˇreˇsen´ı je z´aleˇzitost´ı rutinn´ı. To vˇse za pˇredpokladu, ˇze (k−1) akk = 6 0. 48

Tento pˇredpoklad je rozhoduj´ıc´ı. Numerick´ a anal´ yza spojen´a s t´ımto pˇredpoladem si vynut´ı posl´eze u ´pravu vlastn´ıho algoritmu ˇreˇsen´ı soustav line´arn´ıch algebraick´ ych rovnic - t. zv. pivotn´ı strategie. Jeˇstˇe neˇz se budeme numerick´ ymi aspekty eliminaˇcn´ıch metod zab´ yvat, ˇ uvedme sloˇzitost v´ ypoˇctu dle v´ yˇse uveden´eho schematu za pˇredpokladu, ˇze proces je algoritmicky provediteln´ y, t.j. (k−1)

akk

6= 0, k = 1, ..., N − 1. (0)

Transformace matice B(p) (= B(p) ) do lichobˇeˇzn´ıkov´eho tvaru stoj´ı asymptoticky N X 1 1 (N − j + p)(N − 1) ∼ N 3 + pN 2 . 3 2 j=1

Je - li na dan´em poˇc´ıtaˇci ˇcasov´ a spotˇreba na jednu aritmetickou operaci ν sekund, pak stanoven´ı ˇreˇsen´ı dan´e soustavy s matic´ı typu N × N pro p prav´ ych stran se realizuje v n´asleduj´ıc´ıch ˇcasov´ ych mez´ıch: N

1 3 3N

(eliminace)

1 2 2N

(zpˇetn´ a substituce)

10

3, 3.10−4 ν sek

5.10−5 ν sek

102

3, 3.10−1 ν sek

5.10−2 ν sek

103

3, 3.105 ν sek

5.10−1 ν sek

106

3, 3.1012 ν sek

5.105 ν sek

Pozn´ amka 5.1 Jistˇe stoj´ı za povˇsimnut´ı, ˇze spotˇreba ˇcasu uveden´ a v posledn´ım ˇr´ adku ˇcin´ı ˇra ´dovˇe ν.109 hodin! To ukazuje, ˇze pro velk´e dimenze jsou metody eliminaˇcn´ı na soudob´ych seri´ aln´ıch poˇc´ıtaˇc´ıch prakticky nepouˇziteln´e. Uveden´a pozn´amka osvˇetluje nejr˚ uznˇejˇs´ı snahy obej´ıt komplikace spojen´e ˇ s ˇcasovou n´aroˇcnost´ı eliminaˇcn´ıch metod. Tendence jsou zˇrejm´e - vyuˇz´ıt bud speci´aln´ı struktury poˇc´ıtaˇce (v´ıceprocesorov´e) a nebo speci´aln´ı struktury dan´e tˇr´ıdy matic ˇreˇsen´ ych soustav (pozitivnˇe definitn´ı a p´asov´e). Pivotn´ı strategie (k−1)

Snadno nahl´edneme, ˇze pˇredpoklad nenulovosti diagon´aln´ıch prvk˚ u akk , y a ˇze algoritmus eliminace se zhrout´ı pˇri jeho nesk = 1, ..., N − 2, je podstatn´ plnˇen´ı. Nepˇr´ıjemn´e je t´eˇz, ˇze tot´eˇz lze ˇr´ıci o numerick´e realizaci, t. j. eliminaˇcn´ı 49

(k−1)

algoritmy se zhrout´ı, pakliˇze uveden´e prvky akk nejsou dostateˇcnˇe velk´e v absolutn´ı hodnotˇe. Demonstrujme takovou situaci na pˇr´ıkladˇe n´asleduj´ıc´ı soustavy. Nechˇt tedy  x3 = 1  x1 + x2 + x1 + x2 + 2x3 = 2 (5.28)  x1 + 2x2 + 2x3 = 1 Matice t´eto soustavy je regul´arn´ı a existuje tud´ıˇz jedin´e ˇreˇsen´ı x1 = −x2 = x3 = 1. Avˇsak po proveden´ı prvn´ıho kroku eliminace obdrˇz´ıme š x3 = 1 (5.29) x2 + x3 = 0, a d´ale eliminovat nelze (5.30)

(1)

a22 = 0.

V obecn´em pˇr´ıpadˇe kdy detA 6= 0 a pˇritom ak−1 kk = 0, mus´ı existovat pro nˇejak´e j > k prvek (k−1)

akj

= 6 0,

jinak by k - t´ y ˇr´adek byl nulov´ y, coˇz vyluˇcuje regularitu matice soustavy. Vhodnou z´amˇenou ˇra´dk˚ u ˇci sloupc˚ u soustavy lze posl´eze doc´ılit toho, aby po takov´e u ´pravˇe diagon´aln ´ı prvek byl roven vhodnˇe vybran´emu prvku k - t´eho ˇr´ adku. Tato u ´vaha m´a opr´avnˇen´ı i v situaci, kdy posuzujeme nenulovost z pohledu dan´eho poˇc´ıtaˇce. Jin´ ymi slovy, permutaci promˇenn´ ych ˇci ˇr´ adk˚ u je nutn´e prov´ adˇet (k−1) nen´ı dostateˇcnˇe velk´ a. kdykoliv absolutn´ı hodnota prvku akk Dokumentujme takov´ y pˇr´ıpad soustavou  x2 + x3 = 1  x1 + x3 = 1 0, 0001x + (5.31) 2  9, 999x3 = 10

Matice t´eto soustavy je regul´arn´ı horn´ı troj´ uheln´ıkov´ a. Provedeme - li zpˇetnou substituci uˇz´ıvaj´ıce tˇr´ı platn´ ych cifer v pohybliv´e desetinn´e ˇc´ arce, zjist´ıme, ˇze x ˜1 = 0, x ˜2 = 0, x ˜3 = 1 splˇ nuj´ı poˇzadovan´e vztahy. ciframi) je d´ano v´ yrazy

Avˇsak pˇresn´e ˇreˇsen´ı (vzat´e se tˇremi polatn´ ymi x1 = 1, −x2 = x3 = 1. 50

ˇ je smutn´ Jak vysvˇetlit tento rozpor tedy numerickou nestabilitu? Odpovˇed a, je to v podstatˇe vˇeci: Gaussova eliminaˇcn´ı metoda je numericky nestabiln´ı. V souvislosti s pˇredchoz´ım pˇr´ıkladem si vˇs´ımnˇeme toho, ˇze souˇcasnou z´amˇenou 2. a 3. ˇr´adku a 2. a 3. sloupce, t. j. y1 = x1 , y2 = x3 y3 = x2 , se obdrˇz´ı soustava y1

+

y2 9, 999y2 y2

+

y3

+ 0, 0001y3

= 1 = 10 = 1.

Zde jiˇz y1 = y2 = −y3 = 1, ych takˇze po zpˇetn´e substituci y → x obdrˇz´ıme pˇresn´e ˇreˇsen´ı v r´amci tˇr´ı platn´ cifer. ˇ asteˇcn´ C´ a pivotace Hledejme pˇrirozen´e ˇc´ıslo t tak, aby t = min{j : k ≤ j ≤ N } a pˇritom

Œ Œ n o Œ (k−1) Œ (k−1) Œatk Œ = max |ajk | : j = 1, ..., N .

´ a pivotace Upln´

Hledejme pˇrirozen´a ˇc´ısla t a s takov´ a, aby t = min {j : k ≤ j ≤ N } s = min {j : k ≤ l ≤ N } a pˇritom

Œ Œ Œ nŒ o Œ (k−1) Œ Œ (k−1) Œ Œats Œ = max Œajl Œ : k ≤ j ≤ N, k ≤ l ≤ N .

D˚ uleˇzit´e doporuˇcen´ı plynouc´ı z praktick´ ych zkuˇsenost´ı prav´ı, ˇze ˇc´ asteˇcn´ a pivot´aˇz je postaˇcuj´ıc´ı; u ´pln´ a pivot´ aˇz jednak v´ ypoˇcet zpomaluje, jednak numerickou stabilitu podstatnˇe nezvyˇsuje. Form´aln´ı v´ yraz pˇredchoz´ıch u ´vah je obsahem n´asleduj´ıc´ıch vˇet. 51

Vˇ eta 5.9 (o u ´pln´e pivotaci) Budˇ A matice typu M × N a nechˇt A = 6 0. Potom existuj´ı cel´ a kladn´ a ˇc´ısla t ≥ 1 a s ≤ min (M − 1, N ), element´ arn´ı arn´ı doln´ı permutaˇcn´ı matice Pj , Ql , j = 1, ..., M − 1, l = 1, ...N a element´ troj´ uheln´ıkov´e matice Ej indexu N − j, j = 1, ..., s, tak, ˇze plat´ı vyj´ adˇren´ı Aˆ = Es Ps ...E1 P1 A Q1 ...Qt “ ’ Aˆ11 Aˆ12 , = 0 0

uheln´ıkov´ a matice typu r × r s pˇri ˇcemˇz Aˆ11 je horn´ı troj´ det Aˆ 6= 0. D´ ale pak r je hodnost matice A. asteˇcn´e pivotaci) Budˇ A matice typu M × N a nechˇt s = min Vˇ eta 5.10 (o ˇc´ (M − 1, N ). Potom existuj´ı element´ arn´ı permutaˇcn´ı matice Pj a element´ arn´ı doln´ı troj´ uheln´ıkov´e adˇren´ı matice Ej indexu M − j, j = 1, ..., s takov´e, ˇze plat´ı vyj´ A˜ = Es Ps ...E1 P1 A, pˇri ˇcemˇz A˜ je doln´ı lichobˇeˇzn´ıkov´ a, t. j. a ˜jk = 0, j > k. Jiˇz v´ıme, ˇze pˇri ˇreˇsen´ı soustav typu Ax = b, kde A je matice typu N × N Gaussovou eliminaˇcn´ı metodou je nutn´e prov´ adˇet pivotaci. ˇ podle sch´ematu u Tedy, bud ´pln´e pivotace A = EN −1 PN −1 ...E1 P1 A Q1 QN −1 nebo uˇzit´ım pivotace ˇc´asteˇcn´e A = EN −1 PN −1 ...E1 P1 A. Obˇe tyto strategie lze bez obt´ıˇz´ı pˇren´est na pˇr´ıpad simultann´ıho ˇreˇsen´ı p soustav typu Ax = bj , j = 1, ..., p. V takov´em pˇr´ıpadˇe se potˇrebn´e transformace prov´ adˇej´ı s matic´ı Bp = (A, b1 , ..., bp ) . 52

Zde je zˇrejmˇe podstatn´ y pˇredpoklad o nez´avislosti vektoru bj na ˇreˇsen´ı soustav s indexy l < j. Nejsou vˇsak ˇr´ıdk´e pˇr´ıpady, kdy bj+1 = bj+1 (x1 , ..., xp ) a v tom pˇr´ıpadˇe naznaˇcen´ y postup selˇze. Vˇ eta 5.11 Budˇ A matice typu N × N a Ak nechˇt oznaˇcuj´ı matici typu k × k tvoˇrenou prvky leˇz´ıc´ımi v pr˚ uniku ˇr´ adk˚ u a sloupc˚ u s indexy 1, ..., k. Je - li (5.32) det Ak 6= 0, k = 1, ..., N − 1, pak existuje pr´ avˇe jedna doln´ı troj´ uheln´ıkov´ a matice L = (ljk ) typu N × N takov´ a, ˇze ljj = 1, j = 1, ..., N, a, ˇze a pr´ avˇe jedna horn´ı troj´ uheln´ıkov´ a matice U = (ujk ) typuN × N takov´ A = L U. D˚ ukaz. Nechˇt N = 1. Potom l11 = 1, u11 = a11 a tvrzen´ı plat´ı. Nechˇt tvrzen´ı plat´ı pro vˇsechny matice ˇr´ adu l ≤ k−1. Vyˇsetˇrujme ’ “ a Ak−1 Ak = , cT akk za pˇredpokladu, ˇze Ak−1 = Lk−1 Uk−1 . Poloˇzme Lk = a Uk =

’

’

Lk−1 lT

Uk−1 0T

0 1 u ukk

“

,

“

,

kde a a c jsou zn´am´e vektory a l je vektor typu (k −1)×(k −1), jenˇz je zapotˇreb´ı urˇcit podobnˇe jako vektor u a ˇc´ıslo ukk . az´ıme z relac´ı Abychom urˇcili souˇcin Lk Uk , vych´ (5.33)

Lk−1 Uk−1 = Ak−1 ,

(5.34)

Lk−1 u = a,

(5.35)

lT Uk−1 = cT 53

a (5.36)

lT u + ukk = akk .

Dle indukˇcn´ıho pˇredpoklagu jsou Lk−1 a Uk−1 urˇceny jednoznaˇcnˇe a nav´ıc det Uk−1

= det Lk−1 det Uk−1 = det Ak−1 6= 0,

takˇze (5.34) je jednoznaˇcnˇe ˇreˇsiteln´ a vzhledem k u. Podobnˇe je tomu s vektorem l v (5.35). D´ıky jednoznaˇcnosti l a u je posl´eze jednoznaˇcnˇe urˇcen i prvek ukk pomoc´ı (5.36). ||| Nen´ı - li pˇredpoklad (5.32) splnˇen, vˇeta 5.11 neplat´ı; na pˇr. ’ “ 0 1 A= . 1 1 Pˇredpokl´adejme, ˇze naopak, A = L U, kde L= takˇze

’

l11 l21

LU = ˇ odkud plyne, ˇze bud

’

0 l22

“

, U=

l11 u11 l21 u11

’

u11 0

u12 u22

l11 u12 l21 u12 + l22 u22

“

“

,

.

l11 = 0 nebo u11 = 0, ˇ ale potom bud L=

’

U=

’

a nebo

0

0

l21

l22

0 0

u12 0

“ “

,

tedy spor. Uk´aˇzeme nyn´ı, ˇze existence LU rozkladu pro A je zaruˇcena vˇzdy, pakliˇze lze Gauss˚ uv algoritmus realizovat bez pivotace. Pˇredpokl´adejme tedy, ˇze A je takov´ a, ˇze Gaussova eliminace projde bez pivotace, to znaˇc´ı, ˇze A(N −1) = EN −1 ...E1 A 54

a tud´ıˇz

−1 (N −1) A = E1−1 ...EN , −1 A

takˇze

A = L U, kde

−1 L = E1−1 ...EN −1

a U =A

(N −1)



  =  

(N −1)

a11

0

(N −1)

. . . a1N . . . −1 . . . aN NN



   . 

Protoˇze zˇrejmˇe ljj = 1,j = 1, ..., N , v´ yˇse uveden´e tvrzen´ı plat´ı. Obecnˇe vˇsak se Gaussova eliminace, byˇt pro symetrick´e matice bez permutace ˇr´adk˚ u ˇci sloupc˚ u prov´adˇet ned´a. To ukazuje pˇr´ıklad ’ “ 0 1 . 1  Pro pozitivnˇe definitn´ı symetrick´e matice vˇsak lze Gaussovu eliminaci prov´ adˇet bez pivotace. Pozitivnˇe definitn´ı matice jsou totiˇz charakterizov´ any skuteˇcnost´ı, ˇze plat´ı Vˇ eta 5.12 (Sylvesterovo kriterium) Matice A = (ajk ) typu N × N je pozitivnˇe definitn´ı pr´ avˇe kdyˇz (5.37) 0 < det Ak , k = 1, ...N, kde Ak jsou tvoˇreny prvky leˇz´ıc´ımi v pr˚ useku prvn´ıch k ˇr´ adk˚ u a sloupc˚ u. D˚ ukaz. Nutnost je trivi´aln´ı. Je-li A pozitivnˇe definitn´ı, potom, podle pˇredpokladu plat´ı nerovnosti κ(x, x)N ≤ (Ax, x)N , κ > 0, takˇze volbou x= obdrˇz´ıme, ˇze plat´ı

’

x ˜ 0

“

, x ˜ ∈ Rk ,

˜)k ≥ κ(x, x)k . (Ak x ˜, x K d˚ ukazu postaˇcitelnosti nechˇt naopak plat´ı (5.37). Pˇredpokl´adejme, ˇze existuje vektor y ∈ RN takov´ y, ˇze (5.38)

(Ay, y)N ≤ 0, 55

Je tedy zaruˇcena existence vlastn´ı hodnoty λ ≤ 0 a tud´ıˇz, 0 = det (λI − A) =

(5.39)

N X

k=0

aN −k λk ,

kde sign (aN −k ) = (−1)k , k = 1, ..., N. Vid´ıme tedy, ˇze rovnost (5.39) je s λ ≤ 0 vylouˇcena. T´ım je d˚ ukaz vˇety 5.12 proveden. ||| Pozitivnˇ e definitn´ı matice ˇ A = (ajk ) symetrick´ Bud a pozitivnˇe definitn´ı matice typu N × N . Uk´aˇzeme si, ˇze proveden´ı Gaussovy eliminaˇcn´ı transformace zachov´ av´ a symetrii tˇech hlavn´ıch submatic, kter´e jeˇstˇe nejsou horn´ı troj´ uheln´ıkov´e. Tedy, je - li  (k)  (k) (k) (k) a11 . . . a1k a1k+1 . . . a1N   .     .     .   (k) A = k( (k) (k) , . . . akk ) . . . akN  akk+1  0   (k) (k)  0 . . . 0 ak+1k+1 . . . ak+1N     .  . . . . . . . . . 0

.

.

(k)

0

.

aN k+1

. .

potom tvrd´ıme, ˇze (5.40)

ajl = alj , j, l = k + 1, ..., N.

Nechˇt tedy (5.41)

ajl

(k)

.

(k)

aN N

(k)

(k−1)

(k−1)

= alj

, j, l = k, ..., N.

Protoˇze (k)

(k−1)

ajl = ajl

(k−1)

−

akl

(k−1) a , (k−1) jk akk

plyne z (5.41), ˇze (k−1)

akjl = alj

(k−1)

−

alk

(k−1) akk

(k−1)

akj

= aklj .

Plat´ı tedy (5.40). Indukce d´av´ a obecn´ y v´ ysledek. Pozn´ amka 5.2 Tato skuteˇcnost zˇrejmˇe sniˇzuje poˇcet operac´ı nutn´ych k proveden´ı eliminaˇcn´ı transformace zhruba na polovinu.

56

Ze Sylvesterova kriteria (vˇeta 5.12) plyne, ˇze akk > 0, k = 1, ..., N a d´ale Vˇ eta 5.13 Je - li A symetrick´ a pozitivnˇe definitn´ı matice, potom |ajk |2 ≤ ajj akk , j, k = 1, ..., N.

(5.42)

Tud´ıˇz maxim´ aln´ı prvek matice A leˇz´ı na jej´ı diagon´ ale. ˇ P permutaˇcn´ı matice. Potom D˚ ukaz. Bud (P AP T )T = P AT P T = P AP T je symetrick´a. Pro pevnou dvojici index˚ u j a k existuje element´ arn´ı permutaˇcn´ı matice a a pozitivnˇe definitn´ı matice. Na z´akladˇe P tak, ˇze P AP T je symetrick´ Sylvesterova kriteria potom plat´ı relace ’ “ ajj ajk 2 0 < det = ajj akk − ajk akj akk a tud´ıˇz (5.42). Kdyby 6 k, ajk = Max {alt : l, t = 1, ..., N } , j =

pak by na jedn´e stranˇe, dle pˇredpokladu,

ajk > akk , ajk > ajj a souˇcasnˇe ajk ≤

spor. |||

√

ajj

√ akk ,

Vid´ıme tedy, ˇze Gaussovu eliminaci lze pro pozitivnˇe definitn´ı matice prov´ adˇet vˇzdy a to bez pivotace. Hroz´ı jen, ˇze se pˇri prov´ adˇen´ı transformac´ı vytvoˇr´ı dˇelen´ım mal´ ymi ˇc´ısly prvky velk´e. Jiˇz v´ıme, ˇze ty vˇsak mus´ı leˇzet na diagon´ale. D´ale pak z vyj´adˇren´ı (k−1)

(k) ajj

=

(k−1) ajj

odvod´ıme, ˇze

−

(k−1)

(k)

(k−1)

ajj = ajj

+

|akj

(k−1) akk

57

akj

(k−1) a (k−1) jk akk

| ŒŒ (k−1) ŒŒ (k−1) Œajk Œ ≤ 2ajj .

Vid´ıme, tedy, ˇze r˚ ust prvk˚ u transformovan´ ych matic je nez´avisl´ y na velikosti prvk˚ u transformovan´ ych matic. To n´as vede k z´avˇeru, ˇze Gaussova eliminace bez pivotace je pro pozitivnˇe definitn´ı matice bezpeˇcn´ a. Analog vˇety o LU rozkladu pro pozitivnˇe definitn´ı matice lze formulovat takto: Vˇ eta 5.14 Budˇ A symetrick´ a pozitivnˇe definitn´ı matice. Potom existuje pr´ avˇe jedna horn´ı troj´ uheln´ıkov´ a matice R takov´ a, ˇze plat´ı A = RT R. D˚ ukaz. Nechˇt A = L U, kde ljj = 1, j = 1, ..., N. Tento rozklad existuje na z´akladˇe vˇety 5.11 a vˇety 5.12. Na z´akladˇe tohoto rozkladu odvod´ıme, ˇze a11 = u11 > 0 a ujj =

det Aj , j = 1, ..., N. det Aj−1

Poloˇzme D = diag{u11 , ..., uN N }. Potom

˜, U ˜ = D−1 U, A = LDD−1 U = LDU

pˇri ˇcemˇz Ze symetrie A plyne, ˇze neboli, (5.43)

˜ = I = diagU ˜. diagL ˜ T DL ˜T , AT = A = U ˜ T = L, DL ˜ T = U. U

Poloˇzme d´ale

1

R = D− 2 U. Potom

1 1 ˜T U RT R = U T D − 2 D 2 U = U

a dle (5.43) tedy RT = RL U = A. ||| 58

Tˇ r´ıdiagon´ aln´ı matice Matice A = (ajk ) se naz´ yv´ a p´ asov´ a, jestliˇze jestliˇze existuj´ı indexy p ≥ 1 a q ≥ 1 takov´e, ˇze ajk = 0 pro j > k + p a j − q < k. Speci´alnˇe pro p = q = 1 se takov´e p´asov´e matice naz´ aln ´ı. yvaj´ı tˇr´ıdiagon´ Je zˇrejm´e, ˇze Gaussova eliminace p´asovost zachov´ av´ a. Obecnˇe plat´ı i zde, ˇze pivotace je nezbytn´ a i v pˇr´ıpadˇe symetrick´ ych matic. Vyj´ımkami jsou matice diagon´alnˇe dominantn´ı a matice pozitivnˇe definitn´ı . Vyˇsetˇrujme pˇr´ıpad tˇr´ıdiagon´ aln´ıch matic. Matici   0 0 0 a1 c1 0 . . .  b2 a2 c2 . . . 0 0 0     . . . . . . . .  A= .   0 0 0 . . . bN −1 aN −1 cN −1  0 0 0 . . . 0 bN aN

se snaˇzme vyj´adˇriti ve tvaru  1 0 0  β2 1 0  A=  . . .  0 0 0 0 0 0  α 1 γ1 0  0 α2 γ2   . . .   0 0 0 0 0 0

. . . . . . . . . .

. 0 0 . . . . βN −1 . 0

. . . . . . . . . .

0 . 0 0 . 0 . . . . 0 αN −1 . 0 0

0 0 . 1 βN

0 0 . 0 1 0 0 .

γN −1 aαN

Odtud vypl´ yvaj´ı vztahy



   ×  

    

α1 = a1 , βk = bk αk−1 , αk = ak − βk ak−1 , k = 1, ..., N.

ˇ sen´ı se pot´e obdrˇz´ı obvykl´ a z nich posl´eze algoritmus ˇreˇsen´ı. Reˇ ym postupem, dopˇrednou eliminac´ı a zpˇetnou substituc´ı.

59

ˇ bj , j = Cviˇ cen´ı 5.1 Budˇ A tˇr´ıdiagon´ aln´ı pozitivnˇe definitn´ı matice. D´ ale budte 1, ...; p, p ≥ 2. Takto formulovanou u ´lohu lze ˇreˇsit na pˇr. tˇemito dvˇema odliˇsn´ymi postupy: ame postupnˇe A−1 bj , j = 1. Napˇred urˇc´ıme inverzn´ı matici A−1 a poˇc´ıt´ 1, ..., p. ˇ s´ıme postupnˇe soustavy 2. Reˇ Ax = bj , j = 1, ..., p. Rozhodnˇete kter´y z tˇechto postup˚ u je v´yhodnˇejˇs´ı z hlediska sloˇzitosti algoritm˚ u.

5.4

Metody ˇ reˇ sen´ı soustav line´ arn´ıch algebraick´ ych rovnic zaloˇ zen´ e na principu optimalizace

Metody uveden´e v z´ahlav´ı tohoto odstavce jsou typick´e pro soustavy s pozitivnˇe definitn´ı matic´ı. ˇ tedy H pozitivnˇe definitn´ı matice typu N × N . Vyˇsetˇrujme soustavu Bud Hx = b, b ∈ RN .

(5.44)

Vyˇsetˇrujme funkcion´al f : RN → R1 (5.45)

f (x) = (Hx, x) − 2(b, x), x ∈ RN

a zkoumejme jeho vlastnosti. Snadno se pˇresvˇedˇc´ıme, ˇze (5.46) kde (5.47)

ˆ ), x − x ˆ ), x ∈ RN , ˆ ) − (b, x f (x) = (H(x − x ˆ = b, Hx

takˇze d´ıky kladnosti v´ yrazu (Hy, y), y ∈ RN , existuje ˆ ‰ ˆ ), min f (x) : x ∈ RN = f (ˆ x) = −(b, x (5.48)

ˆ splˇ kde x nuje (5.47). Snadno nahl´edneme, ˇze u ´loha nal´ezt ˇreˇsen´ı soustavy (5.44) s pozitivnˇe definitn´ı matic´ı je ekvivalentn´ı s u ´lohou nal´ezt (lok´aln´ı) minimum funkcion´ alu (5.45). S podobnou situac´ı jsme se jiˇz setkali pˇri vyˇsetˇrov´ an´ı zobecnˇen´ ych ˇreˇsen´ı v 5.1. Zab´ yvejme se tedy ot´azkou sestrojov´ an´ı extr´em˚ u funkcion´ al˚ u na RN . 60

Definice 5.2 Funkcion´ a kvadratick´ y, jestliˇze je polyal f : RN → R1 se naz´yv´ nomem druh´eho stupnˇe jakoˇzto funkce N promˇenn´ych x1 , ..., xN . y funkcion´ al na RN pr´ Snadno zsjist´ıme, ˇze f je kvadratick´ avˇe kdyˇz existuje symetrick´a matice H typu N × N , vektor x ∈ RN a ˇc´ıslo c ∈ R1 , takov´e, ˇze plat´ı 1 f (x) = (Hx, x) + (b, x) + c, x ∈ RN . (5.49) 2 Pozn´ amka 5.3 Z definice (5.45) je patrn´e, ˇze t´ımto vztahem definovan´y funkcion´ al je kvadratick´y. Symbolem C (m) (S), m ≥ 0, oznaˇcujeme mnoˇzinu funkcion´ al˚ u, kter´e maj´ı spojit´e parci´aln´ı derivace aˇz do ˇr´ adu m na mnoˇzinˇe S), kde S ⊃ {x ∈ RN : kx − x0 k < δ}, pro nˇejak´e δ > 0 a x0 ∈ RN . Definice 5.3 Nechˇt f ∈ C (1) (S). Vektor  ∂f   g(x) =   

∂x1

∂f ∂xN

se naz´yv´ a gradientem funkcion´ alu f . Jestliˇze f ∈ C (2) (S), matice H(x) =

’

(x) . . . (x)



    

“ ∂2f (x) ∂xj ∂xk

se naz´yv´ a Hessi´anem funkcion´ alu f . Je-li f ∈ C (1) (S), pak (5.50)

f (x + h) = f (x) + (h, g(x)) + o(khk),

zat´ımco pro f ∈ C (2) (S) (5.51) kde pˇri ˇcemˇz

f (x + h) = f (x) + (h, g(x)) + (H(x)h, h) + o(khk2 ), o(α) = F (x; α), F : RN × R1+ → R1 , lim

α→0

|F (x; α)| = 0. α 61

Pozn´ amka 5.4 Vˇs´ımnˇeme si toho, ˇze na z´ arn´ı funkcion´ al, akladˇe (5.51) kaˇzd´y neline´ maj´ıc´ı gradient a Hessi´ an, se lok´ alnˇe chov´ a jako funkcion´ al kvadratick´y. ˜ ∈ S se naz´yv´ alu f ∈ Definice 5.4 Prvek x a stacion´arn´ım prvkem funkcion´ C (1) (S), jestliˇze plat´ı g(˜ x) = 0, ˜. kde g(x) znaˇc´ı gradient f v x ˆ prvkem lok´ Vˇ eta 5.15 Nechˇt f ∈ C (1) (S). Je-li x aln´ıho minima funkcion´ alu f , pak f je stacion´ arn´ı v prvku x. ˆ a h = −hg(ˆ D˚ ukaz. V (5.50) poloˇzme x = x a x), kde h ∈ R1 je re´aln´ promˇenn´a veliˇcina z nˇejak´eho intervalu [0, h0 ]. Plat´ı tedy f (ˆ x + h) = f (ˆ x) − hkg(ˆ x)k2 + o(h).

ˆ . Potom g(ˆ x) 6= 0 a pro dostateˇcnˇe mal´a h zjist´ıme, Nechˇt f nen´ı stacion´arn´ı v x ˇze plat´ı −hkg(ˆ x)k2 + o(h) < 0, neboli

x + h) < f (ˆ x) f (ˆ ˆ nem˚ a tud´ıˇz x uˇze b´ yt prvkem lok´aln´ıho minima. ||| ˆ nen´ı vˇsak postaˇcuj´ıc´ı podm´ınkou pro extr´em Stacionarita funkcion´alu v x ˆ. vx ˆ ∈ S. Je-li H(ˆ Vˇ eta 5.16 Nechˇt f ∈ C (2) (S) a nechˇt f je stacion´ x) arn´ı v x ˆ prvkem ostr´eho lok´ aln´ıho minima funkcion´ alu f . pozitivnˇe definitn´ı pak je x D˚ ukaz. Z (5.51) a z podm´ınky g(ˆ x) = 0 odvod´ıme, ˇze f (ˆ x) + x + h) = f (ˆ

1 (H(ˆ x)h, h) + o(khk2 ). 2

x) poskytuje existenci γ > 0 takov´eho, ˇze Pozitivn´ı definitnost H(ˆ (H(ˆ x)h, h) ≥ γ(h, h). Tud´ıˇz (5.52)

1 γ(h, h) + o(khk2 ) 2 a protoˇze v´ yraz v prav´e ˇc´asti vztahu (5.52) je kladn´ y pro dostateˇcnˇe mal´a khk2 , ˆ je prvkem ostr´eho minima funkcion´ x alu f . ||| f (ˆ x) ≥ x + h) − f (ˆ

Situace je podobn´a vyˇsetˇrujeme-li f z hlediska lok´aln´ıch maxim. V tom ˆ pˇr´ıpadˇe totiˇz −f m´a v odpov´ıdaj´ıc´ıch prvc´ıch lok´aln´ı minima, takˇze je-li x 62

ˆ a je-li Hessi´an H(ˆ prvkem lok´aln´ıho maxima, pak nutnˇe f je stacion´arn´ı v x x) ˆ je potom prvkem ostr´eho lok´aln´ıho maxima. negativnˇe definitn´ı , x Pro indefinitn´ı Hessi´an nelze ˇz´ adn´e z´avˇery stran extr´em˚ u obecnˇe zaruˇcit. ˆ prvkem lok´aln´ıho extr´emu, pak Hessi´an nem˚ Nen´ı-li x uˇze b´ yt definitn´ı , tedy H(ˆ x) m´a ve sv´em spektru jak kladn´e tak z´aporn´e vlastn´ı hodnoty. Dokonce, ˆ. x) semidefinitn´ı, f m˚ je-li Hessi´an H(ˆ uˇze ale nemus´ı m´ıt lok´aln´ı extr´em v x D´ale si pˇripomeˇ nme klasickou definici z elemnt´ arn´ı matematick´e anal´ yzy. Definice 5.5 Nechˇt f ∈ C (2) (S), x ∈ S ⊂ RN s y ∈ RN , kyk = 1. V´yraz Œ dm f (x + τ y) ŒŒ (m) (5.53) f (x; y) = Œ dτ m τ =0

alu f ˇr´ adu m v prvku x ve smˇeru y. se naz´yv´ a smˇerovou derivac´ı funkcion´ V´ ypoˇcet smˇerov´ ych derivac´ı se prov´ ad´ı pouˇzit´ım obvykl´ ych pravidel derivov´an´ı. Tak pro f ∈ C (1) (S), x ∈ RN a y ∈ RN s kyk = 1, (5.54) D´ale pak, (5.55)

f (1) (x; y) =

N X ∂f (x1 , ..., xN )yj = (g(x), y) . ∂x j j=1

f (2) (x; y) = (H(x)y, y) .

Ze Schwarzovy nerovnosti odvod´ıme pro f ∈ C (1) (S), ˇze Œ Œ nŒ o Œ Œ Œ Œ Œ ˆ) , ˆ )Œ = (g(x), y max Œf (1) (x; y)Œ : y ∈ RN , kyk = 1 = Œf (1) (x; y

kde

ˆ= y

g(x) . kg(x)k

Z (5.54) je patrn´e, ˇze f (1) (x; y) = 0 pro vˇsechny smˇery y nastane pr´avˇe kdyˇz g(x) = 0, tedy, je-li x stacion´arn´ım prvkem funkcion´ alu f . Z (5.55) odvod´ıme zase, ˇze f (2) (x; y) > 0 pro vˇsechny smˇery y = 6 0 pr´avˇe kdyˇz H(x) je pozitivnˇe definitn´ı . Pˇredchoz´ı vˇety lze posl´eze pˇreformulovat takto. ˆ ∈ S je prvkem lok´ aln´ıho minima pro f . Vˇ eta 5.17 Nechˇt f ∈ C (1) (S), a x Potom f (1) (ˆ x; y) = 0 pro vˇsechny smˇery y ∈ RN . 63

ˆ ∈ S je takov´y, ˇze f (1) (ˆ Vˇ eta 5.18 Nechˇt f ∈ C (2) (S) a x x; y) = 0 pro vˇsechny N smˇery y ∈ R . ˆ je prvkem ostr´eho lok´ Potom x aln´ıho minima pro f v S, jestliˇze f (2) (ˆ x; y) > N 0 pro vˇsechny smˇery y ∈ R s kyk = 1. Pozn´ amka 5.5 Vˇety 5.17 a 5.18 umoˇzn ˇuj´ı snadn´ a zobecnˇen´ı na pˇr´ıpady, kdy vyˇsetˇrovan´e funkcion´ aly operuj´ı na obecnˇejˇs´ıch prostorech neˇz jsou aritmetick´e vektorov´e prostory. Definice 5.6 Pˇredpokl´ an na S ⊂ RN a nechˇt adejme, ˇze funkcion´ al f je definov´ 1 η ∈ range f = {α ∈ R : ∃x ∈ S ⇒ f (x) = α}. Mnoˇzina Lα = {x ∈ S : f (x) = α} se naz´yv´ au ´rovˇ novou plochou funkcion´ alu f pro hodnotu α.

Snadno se verifikuje n´azor, kter´ y napov´ıd´ a , ˇze v okol´ı ostr´eho lok´aln´ıho ˆ s f (ˆ ˆ se u minima x x) = α ´rovˇ nov´e plochy Lα s α bl´ızk´ ymi k k α chovaj´ı jako nadelipsoidy (t. j. obecnˇe jako elipsoidy ve v´ıcerozmˇern´ ych prostorech). Obecnˇe ˇreˇceno, numerick´e metody hled´an´ı ostr´ ych lok´aln´ıch extr´em˚ u se ”chovaj´ı l´epe”, kdyˇz u ´rovˇ nov´e plochy v okol´ı extr´em˚ u jsou bl´ızk´e ploch´ am kulov´ ym a ”h˚ uˇre, jsou-li tyto plochy v´ yraznˇe deformov´ any od tvaru plochy kulov´e. Veliˇcina, kter´a ud´av´ a v´ ystiˇznou m´ıru odchylky u ´rovˇ nov´ ych ploch od tvaru sf´ery, je d´ana v´ yrazem › š sup {kx − yk : x ∈ Lα } (5.56) : y ∈ Sα , δα = inf inf {kx − yk : x ∈ Lα } pˇri ˇcemˇz

ˆ ‰ Sα = IntLα = x ∈ Lα : ∃ρ > 0 such that ∀h ∈ RN khk < ρ ⇒ x + h ∈ Lα . Zˇrejmˇe, δα = 1, α → 0, implikuje, ˇze Lα jsou kulov´e plochy.

Cviˇ cen´ı 5.2 Jednoduˇse urˇc´ıme odchylku u ´rovˇ nov´e plochy Lα od plochy kulov´e pro pˇr´ıpad kvadratick´eho funkcion´ alu charakterizovan´eho pozitivnˇe definitn´ım Hessi´ anem H. Tehdy je totiˇz δα =

λmax , λmin

kde a

λmin = min {λ ∈ σ(H)}

λmax = max {λ ∈ σ(H)} , takˇze δα je rovno spektr´aln´ımu ˇc´ıslu podm´ınˇenosti matice H. V dalˇs´ı ˇc´ asti tohoto ˇcl´ anku jeˇstˇe uvid´ıme, jak´y v´yznam m´ a tato skuteˇcnost v anal´yze konvergence nˇekter´ych numerick´ych metod. 64

Iteraˇ cn´ı metody sestrojov´ an´ı prvk˚ u extr´ em˚ u Velk´e mnoˇzstv´ı numerick´ ych metod hled´an´ı lok´aln´ıho minima funkcion´ al˚ u m´a iteraˇcn´ı charakter a je zpravidla tvaru xk+1 = xk + τk dk ,

(5.57)

v nˇemˇz dk je zkusm´y smˇer zat´ımco τk je parametr zajiˇsˇtuj´ıc´ı lok´aln´ı minimalizaci ˇci alespoˇ n redukci velikosti funkcion´ alu f na dan´e mnoˇzinˇe. Iteraˇcn´ı proces (5.57) vyˇzaduje tedy strategii pro urˇcen´ı dvou typ˚ u objekt˚ u: an´ı f pod´el pˇr´ımky xk + τ dk , −∞ < v´ ybˇer dk a urˇcen´ı τk s ohledem na chov´ τ < +∞. Definice 5.7 Pˇredpokl´ adejme, ˇze pro funkcion´ al f a vektory x a d existuje τ0 > 0 takov´e, ˇze f (x + τ d) < f (x), 0 < τ < τ0 . alu f . Potom naz´yv´ ame smˇer d smˇerem kles´an´ı funkcion´ Vˇ eta 5.19 Nechˇt f ∈ C (1) (RN ) a nechˇt g(x) oznaˇcuje (jako obvykle v tomto nuje relaci odstavci) gradient funkcion´ alu f v x. Jestliˇze vektor d splˇ (d, g(x)) < 0, pak d je smˇerem kles´ an´ı funkcion´ alu f v x. D˚ ukaz. Z (5.50) obdrˇz´ıme, ˇze f (x + τ d) = f (x) + τ (d, g(x)) + o(τ ), a protoˇze podle pˇredpokladu (d, g(x)) < 0, mus´ı platit nerovnost τ (d, g(x)) + o(τ ) < 0 pro dostateˇcnˇe mal´a τ . Pozn´ amka 5.6 Je-li g(x) = 0„ pak nelze stanovit, zda d je ˇci nen´ı smˇerem kles´ an´ı f v x bez dalˇs´ı dodateˇcn´e informace. Vˇ eta 5.20 Nechˇt f ∈ C (1) (RN ), nechˇt (d, g(x)) = 0 a (H(x)d, g(x)) < 0 pro nˇejak´ a x a d (H(x) znaˇc´ı jako obvykle Heessi´ an f ). an´ı f v x. Potom d je smˇerem kles´ D˚ ukaz. Z (5.51) plyne (5.58)

f (x + τ d) = f (x) + (1/2) (H(x)d, g(x)) + o(τ 2 ),

odkud tvrzen´ı plyne. ||| 65

Vˇ eta 5.21 Nechˇt f ∈ C (1) (RN ). Potom mezi vˇsemi zkusm´ymi smˇery d funkcion´ alu a v okol´ı x nejrychleji je d = −g(x). f v x t´ım smˇerem, v nˇemˇz f kles´ D˚ ukaz. Naˇs´ım u ´kolem je minimalizovat smˇerovou derivaci f v x vzhledem ke vˇsem zkusm´ ym smˇer˚ um d. Z (5.54) plyne, ˇze to znamen´a minimalizovat v´ yraz (d, g(x)) pro vˇsechny smˇery d, pro nˇeˇz kdk = 1. Protoˇze |(d, g(x))| ≤ |g((x)|, minimum se nab´ yv´a pro g(x) . d=− kg(x)k |||

D´ale se zab´ yvejme probl´emem urˇcen´ı τk jsou-li zad´any dk a xk . Snaˇzme se minimalizovat f pod´el pˇr´ımky x = xk + τ dk , −∞ < τ < +∞. Pro obecn´e funkcion´aly stanoven´ı parametru τk je pomˇernˇe sloˇzita´ au ´loha a jej´ı ˇreˇsen´ı vyˇzaduje dalˇs´ı iteraˇcn´ı proces. Pro pˇr´ıpad kvadratick´eho funkcion´ alu to vˇsak nen´ı ˇz´adn´ y probl´em. Nechˇt tedy f (x) = (1/2)(Hx, x) + (b, x) + c. Okamˇzitˇe lze ovˇeˇrit, ˇze f (x + τ d) − f (x) = (1/2)τ 2 (Hd, d) + τ (d, g(x)) + c˜, kde konstanta c˜ je nez´avisl´a na τ . Je-li H pozitivnˇe definitn´ı a d = 6 0, potom (Hd, d) > 0 a pro (d, g(x)) τ= (Hd, d) f (x + τ d) nab´ yv´a minima vzhledem k τ ∈ (−∞, +∞). Tedy, nez´avisle na volbˇe dk v (5.57) je pro kvadratick´ y funkcion´ al f € k  d , g(xk ) (5.59) τk = − (Hdk , dk ) vhodn´a hodnota lok´aln´ıho optimalizaˇcn´ıho parametru. Metoda nejvˇ etˇ s´ıho sp´ adu Na z´akladˇe v´ ysledku vˇety 5.21 je pˇrirozen´e kl´ast dk = −g(xk ) v (5.57), takˇze (5.60)

xk+1 = xk − τk g(xk ).

Tento vztah spolu se strategi´ı urˇcov´ an´ı τk definuje metodu nejvˇetˇs´ıho sp´ adu minimalizace obecn´eho funkcion´ alu f .

66

Pro pˇr´ıpad kvadratick´eho funkcion´ alu ve tvaru (5.58) poskytuje (5.59) evidentn´ı zp˚ usob urˇcen´ı τk . Pro (5.58) je metoda nejvˇetˇs´ıho sp´adu urˇcena formulemi (5.61) gk = Hxk − b, (5.62) (5.63)

τk =

(gk , gk ) , (Hgk , gk )

xk+1 = xk − τk gk ,

kde k = 0, 1, ... a gk = g(xk ). Zˇrejmˇe (5.61) lze ps´at ve tvaru (5.64)

gk+1 = gk − τk Hgk ,

a to umoˇzn ˇuje v´ ypoˇcet gradient˚ u v (5.61) - (5.63) rekurentnˇe. Pˇr´ısluˇsn´ y algoritmus se ˇr´ıd´ı formulemi (5.65) (5.66) (5.67)

τk =

(Hgk , gk ) , (gk , gk )

xk+1 = xk − τk gk ,

gk+1 = gk − τk Hgk ,

kde k = 0, 1, ... a g0 = Hx0 − b. Poznamenejme, ˇze algoritmus urˇcen´ y formulemi (5.65) - (5.67) lze implementovat za pouˇzit´ı jedin´eho n´asoben´ı matic´ı H na jednu iteraci na rozd´ıl od standardsn´ıho algoritmu (5.61) - (5.63), kdy n´asoben´ı matic´ı H jsou zapotˇreb´ı dvˇe. Anal´ yza konvergence metody nejvˇetˇs´ıho sp´adu d´av´ a podnˇet k obecnˇejˇs´ı u ´vaze. Vzhledem k existenci obecnˇe nekoneˇcnˇe mnoha norem jest moˇzno kl´ast si ot´azku, kterou normu zvolit k odhadu chyby metody tak, abychom obdrˇzeli v´ ysledek co moˇzno nejpˇrirozenˇejˇs´ı, nejrealistiˇctˇejˇs´ı, nejvhodnˇejˇs´ı, nejjednoduásv ˇs´ı, nejatd. Na prvn´ı pohled by se mohlo zd´at, ˇze v´ ybˇer normy nen´ı nijak d˚ uleˇzit´ y; vˇzdyˇt N vˇsechny normy na R jsou ekvivalentn´ı. Pˇri bliˇzˇs´ım zkoum´ an´ı zjist´ıme, ˇze pro nˇekter´e normy je takov´a anal´ yza dost obt´ıˇzn´ a. Rozhoduj´ıc´ım kriteriem pro volbu normy bude tedy co nejmenˇs´ı analytick´ a obt´ıˇznost vyˇsetˇrov´ an´ı konvergence. Tu se ukazuje, jako velice v´ yhodn´a a vskutku pˇrirozen´ a t.zv. energetick´ a norma. ˇ Odvodme tedy vˇetu o rychlosti konvergence metody nejvˇetˇs´ıho sp´adu v energetick´e normˇe.

67

Definice 5.8 Budˇ H pozitivnˇe definitn´ı matice typu N × N . Norma k.kH definovan´ a pomoc´ı skal´ arn´ıho souˇcinu φH (x, y) = (Hx, y) , kde (x, y) = (Hx, x).

PN

k=1 (x)k (y)k ,

se naz´yv´ a H - energetickou normou; tedy, kxk2 =

a pomoc´ı metody nejvˇetˇs´ıho sp´ adu Vˇ eta 5.22 Nechˇt {xk } je posloupnost z´ıskan´ ˆ znaˇc´ı prvek glob´ aln´ıho minima funkcion´ alu pro kvadratick´y funkcion´ al f . Nechˇt x f. Potom plat´ı odhad k

(5.68)

ˆ kH ≤ kxk − x

(κ(H) − 1)

(κ(H) + 1)

k

kx0 − x ˆkH .

D´ ale nechˇt pro libovoln´e  > 0 p() znaˇc´ı nejmenˇs´ı index takov´y, ˇze kxk − x ˆkH ≤ kx0 − x ˆkH , pro kaˇzd´e x0 ∈ RN , potom plat´ı nerovnost (5.69)

1 1 κ(H)log . 2 

p() ≤

D˚ ukaz. Pˇredevˇs´ım si uvˇedomme, ˇze plat´ı  € ˆ ) = H −1 gk , gk . (5.70) E(x − x Platnost (5.70) plyne odtud, ˇze

1 ˆ ), x − x ˆ) , (H(x − x 2

f (x) − f (ˆ x) =

ˆ ) = 2{f (x) − f (ˆ E(x − x x)} a ˆ ). g(x) = H(x − x Z vyj´adˇren´ı x) = f (x) − f (ˆ totiˇz plyne, ˇze



   g(x) =   

1 ˆ ), x − x ˆ) (H(x − x 2

∂f ∂x1 (x)

. . .

∂f ∂xN

(x) 68



  1  ˆ)  = H(x − x  2 

a tud´ıˇz ˆ = H −1 g(x). x−x

D´ale pak

€  ˆ ) = (H(x − x ˆ ) = H −1 g(x), g(x) . ˆ ), x − x E(x − x

Poloˇz´ıme-li tedy x = xk+1 a pouˇzijeme-li rekurentn´ı formule pro prvky gk , obdrˇz´ıme, ˇze €  ˆ ) = H −1 gk+1 , gk+1 E(xk+1 − x  € = H −1 (gk − τk Hgk ), gk − τk Hgk

 € = H −1 (I − τk H)gk , (I − τk H)gk

€  = H −1 (I − τk H)2 gk , gk

€  = H −1 gk , gk − 2τk (gk , gk ) + τk2 (Hgk , gk ) V´ ysledkem je vztah

€  = H −1 gk , gk −

ˆ ) − E(xk − x ˆ) = E(xk+1 − x

(gk ,gk )2 (Hgk , gk ). (H gk ,gk )2

(gk , gk )2 (Hgk , gk )(H −1 gk , gk )

ˆ ), E(xk − x

odkud pak ˆ) = E(xk+1 − x

š

1−

(gk , gk )2 k (Hg , gk )(H −1 gk , gk )

›

ˆ ). E(xk − x

Aplikujeme-li nyn´ı Bergstromovo lemma na pozitivnˇe definitn´ı matici H, maj´ıc´ı vlastn´ı ˇc´ısla λ1 , ..., λN , podle nˇehoˇz 1≤

(Hx,x) (x,H −1 x) (x,x) (x,x)

≤ obdrˇz´ıme odhad E(x

k+1

(λ1 +λN )2 4λ1 λN ,

x ∈ RN ,

”

• 4λ1 λN ˆ) ≤ 1 − ˆ ), E(xk − x −x (λ1 + λN )2

coˇz vyj´adˇreno pomoc´ı spektr´aln´ıho ˇc´ısla podm´ınˇenosti κ(H) = (λN /λ1 ) poskytuje odhad “2 ’ κ(H) − 1 ˆ ). ˆ) ≤ E(xk − x E(xk+1 − x κ(H) + 1 69

Vid´ıme tedy, ˇze konvergenci metody nejvˇetˇs´ıho sp´adu lze popsat pomoc´ı energetick´e normy vztahem ’ “ κ(H) − 1 k+1 ˆ ˆ kH , − xkH ≤ kxk − x kx κ(H) + 1 kter´ ym je dok´az´ana platnost prvn´ıho tvrzen´ı dokazovan´e vˇety. Abychom odvodili platnost nerovnosti (5.69) uvˇedomme si, ˇze platnost nerovnosti ˆ kH pro k ≥ kˆ je zaruˇcena platnost´ı nerovnosti ˆ kH ≤ kx0 − x kxk − x ’

κ(H) − 1 κ(H) + 1

“k

≤

ˆ pro k ≥ k. Hledejme tedy nejmenˇs´ı k takov´e, aby ’ “k 1 κ(H) + 1 ≤ .  κ(H) − 1 Jest totiˇz log

1 ≤ k log(1 + σ), 

kde σ= Protoˇze

2 . κ−1

1 1 1 1 log ≤ (κ − 1) log , σ  2  nuje nerovnost odvod´ıme snadno, ˇze hledan´ y index p() splˇ p() ≤

1 1 κ(H)log + 1. 2 

T´ım je dok´az´ana platnost vztahu (5.69) a tedy t´eˇz vˇety 5.22. ||| Pˇ redpodm´ınˇ en´ı. Pˇredpokl´adejme, ˇze C je pozitivnˇe definitn´ı N × N matice faktorizovan´ a na tvar C = EE T a ˇze 1 (5.71) f (x) = (Hx, x) + (x, b) + c, 2 ˇ kde H je pozitivnˇe definitn´ı N ×N matice. Zavedme dalˇs´ı funkcion´ al ˜f pˇredpisem (5.72)

˜f (y) = f (E T y) = 1 (Hy, ˜ y) + (y, b) + c˜, 2 70

pˇri ˇcemˇz (5.73)

˜ = E −1 b, c˜ = c. ˜ = E −1 H(E T )−1 , b H

˜ symetrick´a a nav´ıc ze vztahu (Hy, ˜ y) = (Hx, x) > 0 pro vˇsechna Zˇrejmˇe je H ˜ je pozitivnˇe definitn´ı. x = E T y, plyne, ˇze H Podobnostn´ı transformace ˜ T = (E T )−1 E −1 H = C −1 H (E T )−1 HE ˜ a C −1 H maj´ı totoˇzn´ ukazuje, ˇze matice H a spektra. Oznaˇcme vlastn´ı hodnoty ˜ ˇ ıslo podm´ınˇenosti matice H symboly λ˜1 , ..., λ˜N . C´ ˜ = κ(H)

(5.74)

λ˜N λ˜1

˜ z´avis´ı jeˇstˇe na faktorizaci matice C. je urˇceno maticemi C a H pˇresto, ˇze H Aplikujme nyn´ı metodu nejvˇetˇs´ıho sp´adu na funkcion´ al (5.72). Obdrˇz´ıme tak tento postup ˜ ˜ k − b, ˜ k = Hy (5.75) g (5.76) (5.77)

τ˜k =

˜k ) (˜ gk , g , ˜g ˜k ) ˜k , g (H

˜ k , k = 0, 1, ... yk+1 = yk − τ˜k g

kde y0 lze volit libovolnˇe. Jiˇz v´ıme, ˇze plat´ı (5.78)

˜ −1 b, ˆ=H lim yk = y

k→∞

˜ a ˇz rychlost konvergence je urˇcena ˇc´ıslem podm´ınˇenosti κ(H). k T −1 k k k Poloˇzme x = (E ) y a g = Hx −b, k = 0, 1, ... Snadno se pˇresvˇedˇc´ıme, ˇze plat´ı vztahy ˜ k = (E T )−1 gk , τ˜k = g

(gk , hk ) , yk+1 = E T (xk − τ˜k hk ), (Hgk , gk )

kde hk = C −1 gk . Jelikoˇz yk+1 = E T xk+1 , v´ yˇse uveden´e vztahy lze realizovat t´eˇz pomoc´ı sch´ematu (5.79) gk = Hxk − b, (5.80)

(5.81) (5.82)

hk = C −1 gk ,

τk =

(hk , gk ) , (Hhk , gk )

xk+1 = xk − τk hk , k = 0, 1, ... 71

Extr´emn´ı vektory pro f a ˜f jsou d´any v´ yrazy ˜ ˜ −1 b, ˆ = H −1 b, a y ˆ=H x pˇri ˇcemˇz ˆ = ET x ˆ. y D´ale pak ˆ) ˆ = E T (xk − x yk − y a ˆ kH˜ . ˆ kH˜ = kxk − x kyk − y Na z´akladˇe vˇety 5.22 vˇsak plat´ı ˆ kH˜ kyk − y =



˜ κ(H)−1 ˜ κ(H)+1

‘k

≤



˜ κ(H)−1 ˜ κ(H)+1

‘k

ˆ kH˜ kyk − y

ˆ kH . kxk − x

Posloupnost {xk }, kde x0 lze volit libovolnˇe, konstruovan´ a pomoc´ı vztah˚ u ˜ ˆ a jej´ı rychlost konvergence je urˇcena ˇc´ıslem κ(H). (5.79) - (5.82), konverguje k x ˜ < κ(H), pak odhad (5.69) napov´ıd´ Lze-li naj´ıt matici C tak, aby κ(H) a, ˇze rychlost konvergence (5.79) - (5.82) je lepˇs´ı neˇz (5.65) - (5.67). Matice C se naz´ yv´a matice pˇredpodm´ınˇen´ı metody nejvˇetˇs´ıho sp´ adu. Metoda urˇcen´a vzorci (5.75) - (5.78) se naz´ yv´ a transformovanou pˇredpodm´ınˇenou metodou nejvˇetˇs´ıho sp´ adu, zat´ımco metoda urˇcen´ a vzorci (5.79) - (5.82) se naz´ yv´ a netransformovanou pˇredpodm´ınˇenou metodou nejvˇetˇs´ıho sp´ adu. N´azvy jsou odvozeny ˆ = ET x ˆ , pro metodu danou vzorci (5.75) - (5.78) je transodtud, ˇze extr´em y ˆ . Plat´ı t´eˇz, ˇze xk = formovan´ ym extr´emem pro (5.79) - (5.82), tedy pro x T −1 k (E ) y , a tuto transformaci lze ale nen´ı nutn´e prov´ adˇet pro kaˇzd´e k = 0, 1, ... V´ ypoˇctov´ y proces prov´adˇen´ y na z´akladˇe v´ yˇse uveden´ ych algoritm˚ u se budeme snaˇzit ocenit z pozic jejich sloˇzitosti. Potˇrebn´e gradienty lze stanovit rekurzivnˇe. Jmenovitˇe pak (5.75) a (5.79) lze nahradit formulemi ˜g ˜ k−1 − τ˜k−1 H ˜k = g ˜ k−1 g a gk = τk−1 Hhk−1 . ˜ a Hhk−1 jsou k dispozici z pˇredchoz´ıho iteraˇcn´ıho kroku. ˜ gk−1 Vektory H N´as zaj´ımaj´ı pˇr´ıpady, v nichˇz rozmˇer matice H je znaˇcn´ y a matice H sama je ˇr´ıdk´a. Matice E b´ yv´a zpravidla ˇr´ıdk´ a troj´ uheln´ıkov´ a matice. Pˇri realizaci 72

˜ neb´ yv´ a ˇr´ıdk´ a a tud´ıˇz se transformovan´e metody nar´aˇz´ıme na skuteˇcnost, ˇze H explicitn´ımu v´ ypoˇctu snaˇz´ıme vyhnout. Snadno nahl´edneme, ˇze v r´amci jednoho iteraˇcn´ıho kroku je v´ ypoˇcet vektoru ˜g ˜ k v (5.76) nejpracnˇejˇs´ı. Je tud´ıˇz z´ahodno prov´ H adˇet jeho v´ ypoˇcet nepˇr´ımo. ˇ s´ıme proto napˇred soustavu Reˇ E T z = gk vzhledem k z, pot´e poloˇz´ıme z∗ = Hz a ˇreˇs´ıme Ez∗∗ = z∗ ˜ −1 gk . Je patrn´e, ˇze abychom tento vzhledem k z∗∗ = E −1 H(E T )−1 gk = H vektor urˇcili, mus´ıme ˇreˇsit dvˇe soustavy s troj´ uheln´ıkov´ ymi maticemi a prov´ezt jedno n´asoben´ı matic´ı H. D´ale pak se mus´ı prov´ezt transformace xk = (E T )−1 yk a to opˇet vyˇzaduje ˇreˇsen´ı jedn´e soustavy s troj´ uheln´ıkovou matic´ı. Nejpodstatnˇejˇs´ı sloˇzkou v´ ypoˇctu podle netransformovan´e metody (5.79) (5.82) je krok (5.80), kde se poˇc´ıt´ a C −1 gk a krok (5.81) kde se stanovuje Hhk . Vektor C −1 gk se vypoˇc´ıt´av´a tak, ˇze se ˇreˇs´ı Ez = gk vzhledem k z a E T hk = z se ˇreˇs´ıvzhledem k hk . Potom zˇrejmˇe hk = C −1 gk . I v netransformovan´e metodˇe se ˇreˇs´ı dvˇe soustavy s troj´ uheln´ıkov´ ymi maticemi. Mohlo by se zd´at, ˇze netransformovan´ a metoda je v´ yhodnˇejˇs´ı protoˇze ˆ . Existuj´ı vˇsak vhodn´e transformaˇcn´ı poskytuje pˇr´ımo hledan´ y extr´emn´ı vektor x matice, pro nˇeˇz v´ ypoˇctov´a n´aroˇcnost nen´ı vˇetˇs´ı neˇz pro metodu netransformovanou. Toto tvrzen´ı se pokus´ıme ozˇrejmit pro pˇr´ıpad pˇredpodm´ınˇen´e metody sdruˇzen´ ych gradient˚ u, tedy metody znaˇcnˇe efektivnˇejˇs´ı neˇz je metoda nejvˇetˇs´ıho sp´adu. Tvrzen´ı samo je vˇsak pravdiv´e i pro metodu nejvˇetˇs´ıho sp´adu.

5.5

Zobecnˇ en´ e Bergstromovo lemma pro pozitivnˇ e definitn´ı matice

Pˇri vyˇsetˇrov´an´ı konvergence nˇekter´ ych metod ˇreˇsen´ı soustav line´arn´ıch algebraick´ ych rovnic bude uˇziteˇcn´e n´asleduj´ıc´ı lemma. Lemma 5.4 Budˇ H pozitivnˇe definitn´ı matice typu N × N, N ≥ 1. Nechˇt pro jej´ı spektrum σ(H) = {λj : j = 1, ..., s}, s ≤ N plat´ı vztahy λmin = λ1 < λj < λj+1 < λs = λmax , j = 1, ..., s. 73

Potom pro libovoln´y vektor x ∈ RN plat´ı nerovnosti 2

(5.83) (x, x) ≤ (Hx, x) (H

−1

1 x, x) ≤ λmin λmax

’

λmin + λmax 2

“2

(x, x)2 .

Pozn´ amka 5.7 Lemma 5.4 pod´ av´ a elegantn´ı odhad doln´ı hranice ˇc´ısla podm´ınˇenou sti matice H dan´e veliˇcinou V (x) = (Hx, x) (H −1 x, x) pomoc´ı pod´ılu ˇctverc˚ aritmetick´eho a geometrick´eho pr˚ umˇeru jej´ıch extr´emn´ıch vlastn´ıch hodnot.

D˚ ukaz provedeme za pouˇzit´ı n´asleduj´ıc´ıch dvou lemmat. Lemma 5.5 Nechˇt 0 < µ < ν a f (x) = µ

1 µ + νx, a ≤ x ≤ 1, 0 < a, 2 > ν. x a

Potom (5.84)

1 f (x) ≤ µ + νa. a

u D˚ ukaz. Snadno nahl´edneme, ˇze ze vztah˚ f 0 (x) = −

µ 2µ + ν = 0, f 00 (x) = 2 > 0, x2 x

plyne, ˇze max{f (x) : x ∈ [a, 1]} = f (a) = a odtud tvrzen´ı snadno plyne. |||

74

µ + νa a

Definujme dvˇe veliˇciny φ(ζ) =

(5.85)

s X

ζj λj , s = cardσ(T ),

j=1

ψ(ζ) =

(5.86)

s X

ζj λ−1 j .

j=1

asleduj´ıc´ı tvrzen´ı Lemma 5.6 Pro veliˇciny definovan´e v (5.85) a (5.86) plat´ı n´ (a), (b), (c): (a) λmin λmax ψ(ζ) ≤ (λmin + λmax ) − φ(ζ), (b) Existuje λ ∈ [λmin , λmax ] takov´e, ˇze φ(ζ) = λ. (c)

  €  !2 1 λmin + λmax 2€p  . φ(ζ)ψ(ζ) ≤ λmin λmax

Proof. Abychom uk´azali platnost tvrzen´ı (a), vyˇsetˇrujme v´ yrazy φ(ζ) + λmin λmax ψ(ζ) = λmin

s X

ζj

j=1

s X

ζj

j=1

”

λj λmin

+ λmax

s X

λ ζj min = λj j=1

• λmin + λmax λmin . λmin λj λj

Aplikujeme-li Lemma 5.5 s a=

λmin , µ = λmin , ν = λmax , λmax

obdrˇz´ıme, ˇze φ(ζ) + λmin λmax ψ(ζ) ≤

s X j=1

€  ζj λmin + λmax = λmin + λmax .

Tud´ıˇz, (a) plat´ı. Protoˇze plat´ı i (b).

λmin ≤ φ(ζ) ≤ λmax ,

75

Koneˇcnˇe, piˇsme (na z´akladˇe (a) a (b)), €  φ(ζ)ψ(ζ) ≤ λ λmin + λmax − λ ≤ ’ “  2 1€ 1 λmin + λmax λmin λmax 2 T´ım je d˚ ukaz lemmatu 5.5 proveden. ||| D˚ ukaz lemmatu 5.4. Je li s = 1, pak H = λI, a tvrzen´ı plat´ı. Nechˇt tedy s > 1. Vyjdˇeme ze spektr´aln´ıho vyj´adˇren´ı H=

s X

λj Pj ,

j=1

v nˇemˇz

PjT = Pj = Pj2 , s X

Pj Pk = Pk Pj = δjk Pj , j, k = 1, ..., s,

Pj = I.

j=1

Pro 0 6= x ∈ RN jest tedy, V (x) = (Hx, x) (H −1 x, x) =

X

λj (Pj x, x)

j

X

λ−1 k (Pk x, x)

k

a d´ale pak (5.87)

V (x) =

X

2

(Pj x, x) +

j

j

takˇze V (x) ≥

X X ’ λj

λk + λk λj

k>j

“

(Pj x, x) (Pk x, x),

X XX (Pj x, x) (Pk x, x) = (Pj x, x)2 + 2 j

j

 

k>j

2

X (Pj x, x) = (x, x)2 . j

Na druh´e stranˇe, poloˇzme

ζj = (Pj x, x), j = 1, ..., s a aplikujme lemma 5.6. T´ım vˇsak obdrˇz´ıme platnost tvrzen´ı dokazovan´eho v lemmatu 5.4. |||

76

Pozn´ amka 5.8 Je zˇrejm´e, ˇze pro H = λI, se nabude extr´emn´ı doln´ı (i horn´ı) hranice (Hx, x) (H −1 x, x) = (x, x)2 . Podobnˇe, pro H= se nab´yv´ a extr´emn´ı horn´ı hranice

’

λ1 0

0 λ2

“

,

’

“ λ1 λ1 + + + x12 x22 = (Hx, x) (H x, x) = λ2 λ2 “2 ’ 1 1 (x21 + x22 )2 (λ1 + λ2 ) . 2 λ1 λ2 −1

x14

x42

Cviˇ cen´ı 5.3 Stanovte vˇsechny pˇr´ıpady, pozitivnˇe definitn´ıch matic H, pro nˇeˇz nastane nˇekter´ a z extr´emn´ıch situac´ı ve formuli (5.83). Zobecnˇ en´ e Bergstromovo lemma pro oper´ atory na nekoneˇ cnˇ e dimenzion´ aln´ım Hilbertovˇ e prostoru Pˇredpokl´adejm,e, ˇze E je Hilbert˚ uv prostor a ˇze B(E) znaˇc´ı prostor ohraniˇcen´ ych line´arn´ıch oper´ator˚ u zobrazuj´ıc´ıch E do sebe.

ator. Lemma 5.7 Budˇ T ∈ B(E) samoadjungovan´y pozitivnˇe definitn´ı oper´ Potom pro libovoln´y prvek x ∈ E plat´ı vztahy “ ’ λmin + λmax 2 1 (5.88) (x, x)2 ≤ (Hx, x) (H −1 x, x) ≤ (x, x)2 , λmin λmax 2

pˇri ˇcemˇz λmin a λmax znaˇc´ı hranice spektra oper´ atoru T , t. j. σ(T ) ⊂ [λmin , λmax ] a λmin , λmax ∈ σ(T ). D˚ ukaz. Podobnˇe jako formulace tak i d˚ ukaz lemmatu 5.7 je v podstatˇe totoˇzn´ y s d˚ ukazem lemmatu 5.4. ˇ E = E(λ) spektr´aln´ı m´ıra oper´atoru T . Bud ukaz Vzhledem k pˇredpokladu o nekoneˇcnosti dimenze E je nutn´e upravit d˚ ˇc´aasti (a) v lemmatu 5.6, pˇri ˇcemˇz klademe Z λmax Z λmax φ(x) = λ(dE(λ)x, x), ψ(x) = λ−1 (dE(λ)x, x). λ λ min min 77

Plat´ı tedy, φ(x) + λmin λmax ψ(x) =

Z λmax ” λ

min

• λ λ + λmax min (dE(λ)x, x) λmin λmin λ

a tedy t´eˇz  R  φ(x) + λmin λmax ψ(x) ≤ (λmin + λmax ) λλmax (dE(λ)x, x)  min (5.89)   = (λmin + λmax )(x, x).

Nerovnost ve vztahu (5.89) jsme odvodili pomoc´ı lemmatu 5.5, kde jsme poloˇzili λ a = min , µ = λmin , ν = λmax . λmax Dok´azali jsme tak platnost tvrzen´ı analogick´e k (a) v lemmatu 5.6. Ostatn´ı ˇcx´asti d˚ ukaz lemmatu 5.7 jsou shodn´e s odpov´ıdaj´ıc´ımi ˇc´ astmi d˚ ukazu lemmatu 5.6. T´ımto konstatov´an´ım lze ukonˇcit d˚ ukaz lemmatu 5.7. |||

6

Probl´ emy vlastn´ıch hodnot

Jest na m´ıstˇe si uvˇedomit, ˇze u ´loha nal´ezt vlastn´ı ˇc´ıslo a pˇr´ısluˇsn´ y vlastn´ı vektor dan´e ˇctvercov´e matice je u ´loha neline´ arn´ı. Zkuˇsenost n´am napov´ıd´ a, ˇze na rozd´ıl od line´arn´ıch u ´loh, pro kter´e zpravidla jsme schopni nal´ezt finitn´ı metody k nalezen´ı pˇresn´ ych ˇreˇsen´ı, pro neline´arn´ı u ´lohy takov´e metody zn´amy nejsou. Tak je tomu i v pˇr´ıpadˇe hled´an´ı vlastn´ıch ˇc´ısel a vlastn´ıch vektor˚ u. D´a se dokonce ˇr´ıci, ˇze iteraˇcn´ı metody jsou v tomto pˇr´ıpadˇe jedin´ ymi spolehliv´ ymi prostˇredky jak sestrojovat patˇriˇcn´e aproximace hledan´ ych vlastn´ıch prvk˚ u. Z´akladn´ı metodou v problematice vlastn´ıch ˇc´ısel je mocninn´ a metoda.

6.1

Mocninn´ a metoda

Konvergence mocninn´e metody je zaloˇzena na jednoduch´em pozorov´ an´ı, ˇze totiˇz ˇ divergentn´ı a to kdyˇz |λ| ≥ 1, posloupnost {λk }, kde λ je komplexn´ı ˇc´ıslo, je bud λ 6= 1, stacion´ arn´ı, kdyˇz λ = 1 anebo konvergentn´ı k nule, kdyˇz |λ| < 1. Tuto skuteˇcnost reflektuje posloupnost {T k }, kde T je matice typu N × N s komplexn´ımi prvky tjk , j, k = 1, 2..., N . an´ı algoritm˚ u k Protoˇze chov´an´ı posloupnosti mocnin {T k } je pro vyˇstˇrov´ nalezen´ı vlastn´ıch hodnot matice T rozhoduj´ıc´ı, bude uˇzteˇcn´e pˇripomenout si nˇekter´e d˚ uleˇzit´e vlastnosti funkc´ı, jejichˇz argumentem je matice. K tomu n´am poslouˇz´ı zn´am´e tvrzen´ı z line´arn´ı algebry

78

Tvrzen´ı 6.1 Budˇ T ˇc´ısel. Budˇ

= (tjk ) matice typu N × N nad tˇelesem komplexn´ıch σ(T ) = {λ1 , ..., λN }, s ≤ N,

spektrum matice T a nechˇt

r1 , ..., rs ,

s X

rj = N,

j=1

jsou n´ asobnosti vlastn´ıch hodnot λ1 , ..., λs v charakteristick´em mnohoˇcenu χT (λ) = det(λI − T ) =

s Y

j=1

(λ − λj )rj .

D´ ale nechˇt φ oznaˇcuje minim´aln´ı polynom matice T , t.j. φ(λ) =

s Y

(λ − λj )qj ,

j=1

pˇri ˇcemˇz 1 ≤ qj ≤ rj . Potom existuje N × N matice V, detV = a, ˇze 6 0, takov´   J1   .    . V T V −1 =  .    . Js

V tomto vyj´ adˇren´ı je kaˇzd´y z blok˚ u Jjrj s´ am o sobe blokovˇe diagon´ aln´ı mat¯ aln´ı bloky skal´ ice. Pro pˇr´ıpad qj = 1 jsou jej´ı diagon´ arn´ı matice t.j. Jjl = λj I, kde I je jednotkov´ a matice pˇr´ısluˇsn´eho rozmˇeru. Pro qj > 1 mohou nˇekter´e diagon´ aln´ı bloky b´yt skal´ arn´ı matice avˇsak alespoˇ n jeden z blok˚ u m´ a tvar   λj 1   . .     . . Jjl =  ,  . 1  λj dimJjl = tjl

a ˇc´ısla, (tedy nikoliv nula) a plat´ı rovnosti pˇri ˇcemˇz t1 , ..., ts jsou pˇrirozen´ tj X

dimJjl = rj ,

l=1

max{dimJjl : j = 1, ...s; l = 1, .., tj } = qj . 79

Pˇ r´ıklad 6.1 Nechˇt s = N . Potom qj = rj = 1 a vˇsechny bloky Jjl , j = 1, ..., N , jsou nutnˇe diagon´ aln´ı, takˇze V T V−1 = diag{λ1 , ..., λN }.

Pˇ r´ıklad 6.2 Nechˇt T = λI, N > 1. Zˇrejmˇe σ(T ) = {λ}. Tud´ıˇz s = 1 < N a r1 = N , zat´ımco q1 = 1. Pˇ r´ıklad 6.3 Nechˇt N > 1 a



λ 1  .  T =  



  .  1  λ

. . . .

Zˇrejmˇe, σ(T ) = {λ}, s = 1 a r1 = N . Na rozd´ıl od pˇr´ıkladu 6.1, q1 = N . Definice 6.1 Matice T typu N × N m´ a jednoduchou strukturu, jestliˇze qj = 1 pro j = 1, ..., s, kde s znaˇc´ı poˇcet r˚ uzn´ych vlastn´ıch hodnot matice T . Pozn´ amka 6.1 Matice v pˇr´ıkladech 6.1 a 6.2 maj´ı jednoduchou strukturu; matice z pˇr´ıkladu 6.3 jednoduchou strukturu nem´ a. ˇ Zavedme n´asleduj´ıc´ı oznaˇcen´ı.  0     −1  Bj1 = V    



. . Irj . . 0

kde Irj je jednotkov´a matice rozmˇeru rj = s X

Ptj

l=1

     V,    

dimJjl . Vid´ıme, ˇze

Bj1 = I,

j=1

kde I = IN . D´ale nechˇt

Bjk



    −1  =V    



0 . . Jj − λj Irj

. . 0

80

     V, j = 1, ...s, k = 1, 2...    

Snadno se pˇresvˇedˇc´ıme, ˇze plat´ı vztahy Bjk+1 = (T − λj I)Bjk = (T − λj I)k Bj1 , tedy, Bjk = 0 pro k > qj . D´ale pak plat´ı formule Bj1 Bt1 = δjt Bj1 , j, t = 1, ..., s a Bjk Btl = δjt Bjk+l−1 ; speci´alnˇe

2 = Bj1 , j = 1, ..., s, Bj1

Podobnˇe se lze pˇresvˇedˇcit, ˇze T =

s X

(λj Bj1 + Bj2 )

j=1

a tedy, Tk =

s X

k

(λj Bj1 + Bj2 )

j=1

=

qj s X (t−1) X f (λj ) k

j=1 t=1

kde

(t − 1)!

Bjt ,

fk (λ) = λk , takˇze fk0 (λ) = kλk−1 , fk00 (λ) = k(k − 1)λk−2 , ........................................... (t) fk (λ)

= k(k − 1)...(k − t + 1)λk−t .

Typickou pro konvergenci mocninn´e metody je n´asleduj´ıc´ı situace. Nechˇt spektrum σ(T ) = {λ1 , ..., λs } m´ a strukturu danou vztahy (6.1) Potom

|λ1 | > |λj |, j > 1, q1 = 1. ’

1 T λ1

“k

= B11 +

qj s X (t−1) X (λj ) g k

j=2 t=1

81

(qj − 1)!

Bjt ,

kde

’

gk (λ) = Je okamˇzitˇe patrn´e, ˇze lim

k→∞

’

1 T λ1

neboˇt

λ λ1

“k

“k

.

= B11 ,

(t)

lim gk (λj ) = 0

k→∞

pro j = 2, ..., s a t = 0, 1, ..., qj − 1, j > 1. Vˇ eta 6.1 Je - li q1 = 1, pak posloupnost (’ “k ) 1 T λ1 konverguje k B11 s rychlost´ı s jakou skal´ arn´ı posloupnost (Œ Œ ) Œ λj Œk q −1 j Œ Œ max Œ Œ k : j = 2, ..., s λ1

konverguje k nule. Je - li q1 > 1, pak plat´ı (6.2)

1

(q1 − 1)! lim

k→∞ λk−q1 +1 1

k −q1 +1 T k = B1q1 .

Pˇredchoz´ı teorie nab´ız´ı n´asleduj´ıc´ı obecn´ y algoritmus mosninn´e metody. Algoritmus 6.1 Budˇ {x0k } posloupnost line´ arn´ıch funkcion´ al˚ u na C N takov´ych, ˇze x0k (B11 x0 ) 6= 0a lim x0k = x0∞ , (6.3) k→∞

(6.4)

T

(B11 )

Budˇ ε > 0 zadan´ a tolerance. 1. Zvolme x0 . 2. Poloˇzme k = 0. 3. Sestrojme vektor xk kladouce

x0∞

6= 0.

xk = T xk . 4. Urˇceme λ(k) pomoc´ı pˇredpisu λ(k) =

x0k (xk ) x0k (xk−1 ) 82

a xk+1 =

1 k x . λ(k)

5. Provˇeˇrme, zda Œš Œ Œ ≥ε Œ Œλ(k) − λ(k−1) Œ <ε

(6.5)

(NE), (ANO)

5. Nastane - li v (6.4) pˇr´ıpad NE, poloˇzme k + 1 → k a GOTO 3. 6. Nastane - li v (6.4) pˇr´ıpad ANO , GOTO 7. 7. Poloˇzme λ1 = λk a STOP. Potˇrebn´e vlastnosti algoritmu 6.1 jsou obsahem n´asleduj´ıc´ı vˇety, jiˇz formulujeme ve zjednoduˇsen´e verzi Vˇ eta 6.2 Budˇ T matice typu N × N , takov´ a, ˇze plat´ı vztahy (6.1), (6.4) a 0 xk = x0 pro k = 0, 1, ... takˇze (6.4) plat´ı automaticky. Potom plat´ı rovnosti lim xk = x∞ , T x∞ = λ1 x∞ ,

(6.6)

k→∞

a (6.7)

lim λ(k) = λ1 .

k→∞

D˚ ukaz. Snadno ovˇeˇr´ıme, ˇze xk+1 =

k Y 1 k+1 T x0 λ(l) l=0

a

k Y

λ

(l)

=

l=0

k Y x0 (T xl ) l=0

= a tud´ıˇz xk+1 =

x0l (x)

x0 (T k+1 x0 ) x0 (x0 )

x0 (x0 ) T k+1 x0 . x0 (T k+1 x0 )

Na z´akladˇe vˇety 6.1 odvod´ıme, ˇze lim xk = x0 (x0 ) lim

k→∞

k→∞

B1q1 x0 + Zk x0 , x0 (B1q1 x0 ) + x0 (Zk )

kde, vzhledem k (6.2) lim Zk x0 = 0.

k→∞

83

Posl´eze tedy, lim xk =

k→∞

1 x0 (B

1q1 x0 )

B1q1 x0 .

Jako d˚ usledek obdrˇz´ıme vztahy lim λ(k) =

k→∞

x0 (T k x0 ) k→∞ x0 (T k−1 x0 )

= lim = λ1 lim

k→∞

x0 (k −q1 +1 ( λ11 T )k x0 ) x0 (k −q1 +1 ( λ11 T )k−1 x0 )

= λ1 .

T´ım je d˚ ukaz vˇety 6.2 proveden. |||

6.2

Algoritmy LR a QR

Definice 6.2 Budˇ A matice typu N × N nad tˇelesem re´ aln´ych ˇc´ısel. LR rozkladem matice A nazveme dvojici matic typu N × N takovou, ˇze plat´ı: (a) L = (ljk ) je doln´ı troj´ uheln´ıkov´ a s jednotkovou diagon´ alou, t.j. š 1 pro j = k ljk = j, k = 1, ...N, 0 pro j < k (b) R = (rjk ) je horn´ı troj´ uheln´ıkov´ a, t.j. rjk = 0 pro j > k, j, k = 1, ..., N. (c) A = LR. Pozn´ amka 6.2 Postaˇcuj´ıc´ı podm´ınkou pro existenci LR rozkladu N × N matice A je pln´a regularita matice A, t.j. vˇsechny hlavn´ı minory matice A jsou regul´ arn´ı. Algoritmus 6.2 Budˇ A matice typu N × N a ε > 0 dan´ a tolerance. (1) Poloˇzme A = A1 a k = 1. (2) Nechˇt Ak = Lk Rk je LR rozklad matice Ak . (3) D´ ale poloˇzme (k)

Rk Lk = Ak+1 , Rk = (rjl ). ˇ (4) Posudme, zda max

Œ nŒ o š ≥ ε NE Œ (k) Œ Œrjl Œ : j, l = 1, ...N, j > l < ε ANO. 84

(5) Je -li v (4) NE, poloˇzme k + 1 → k a GOTO (2). (6) Je - li v (4) ANO, poloˇzme R∞ = Rk+1 a GOTO (7). (7) STOP. Tvrzen´ı 6.2 Nechˇt LR rozklady Ak = Lk Rk

(6.8)

existuj´ı pro k = 1, 2, ... Potom plat´ı Ak+1 a Ak jsou podobn´e matice; jmenovitˇe (6.9)

Ak+1 = L−1 k Ak Lk ,

(6.10)

Ak+1 = (L1 ...Lk )−1 A1 (L1 ...Lk ).

Matice Tk , kde Tk = L1 ...Lk je doln´ı troj´ uheln´ıkov´ a a matice Uk , Uk = Rk ...R1 je horn´ı troj´ uheln´ıkov´ a, pˇri ˇcemˇz (6.11)

Ak = A1k = Tk Uk .

D˚ ukaz. Vztahy (6.9) a (6.10) jsou zˇrejm´e. Pomoc´ı nich pak odvod´ıme, ˇze L1 ...Lk Ak+1 = A1 L1 ...Lk . D´ale pak, Tk Uk = (L1 ... Lk−1 ) (Lk Rk ) (Rk ...R1 ) = (L1 ...Lk−1 ) Ak (Rk−1 ...R1 ) = A1 (L1 ...Lk−1 ) (Rk−1 ...R1 ) = A1 Tk−1 Uk−1 = ...................... = A1k = Ak . |||

85

Matice Tk a Uk tvoˇr´ı LR rozklady pro Ak , neboˇt zˇrejmˇe diag Tk = I. Naˇs´ım c´ılem je uk´azat, ˇze (6.12)

lim Lk = I

k→∞

a



λ1  0  lim Ak = lim Rk =  . k→∞ k→∞ 0

(6.13)

× λ2 . 0

. . . .

kde λ1 , ..., λN jsou ˇc´ısla.

 . . × . . ×  , . . .  . . λN

nuje Vˇ eta 6.3 Nechˇt A = A1 splˇ 10 . LR rozklady matic Rk Lk existuj´ı pro k = 1, 2, ... 20 . Nechˇt |λ1 | > |λ2 | > ... |λN | > 0. (6.14)

30 . Pˇredpokl´ adejme, ˇze pro matice X a Y , kde A = XDY je Jordanova forma matice A, existuj´ı LR rozklady X = LX RX , Y = LY RY , diagLX = diagLY = I.

(6.15)

Potom plat´ı (6.12) a (6.13). D˚ ukaz. Podle pˇredpokladu existuje D−1 takˇze Ak = XDk Y = LX Dk LY RY = LX RX (Dk LY D−k ) Dk RY . Avˇsak

 ‘ (k) Dk LY D−k = ljt

je doln´ı troj´ uheln´ıkov´a matice, pˇri ˇcemˇz ’

λj λt

“k

ljt =

š

1 pro 0 pro

(k)

ljt = a (k)

Protoˇze odvod´ıme, ˇze

ljt , LY = (ljt ) j=t j
|λt | > |λj | pro j > t, (k)

lim ljt = 0 pro j > t.

k→∞

86

Tud´ıˇz D−k LY Dk = I + Ek , pˇri ˇcemˇz lim Ek = 0.

k→∞

Z pˇredchoz´ıho plyne, ˇze Ak

= LX RX (I + Ek ) Dk RY −1 = LX (I + RX Ek RX ) RX D k RY k = LX (I + Fk ) RX D RY ,

kde

−1 . Fk = RX Rk RX

a k existuj´ı LR rozklady Protoˇze Ek → I, pro vˇsechna dostateˇcnˇe velk´ ˜ ˜ ˜ ˜ I + Fk = Lk Rk , diagLk = I, ljt = 0 pro j < t. D´ıky tomu, ˇze Fk → 0, odvod´ıme snadno, ˇze ˜ k = I = lim R ˜k . lim L k→∞

Posl´eze,

k→∞

˜ k ) (R ˜ k RX Dk RY ). Ak = (LX L

Z jednoznaˇcnosti troj´ uheln´ıkov´ ych rozklad˚ u plyne, ˇze ˜k Tk = L1 ...Lk = LX L a

˜ k RX D k RY . Uk = Rk ...R1 = R

˜ k → I, R ˜ k → I, obdrˇz´ıme, ˇze plat´ı vztahy Protoˇze L lim Tk = LX ,

k→∞

−1 lim Lk = lim Tk−1 Tk = I,

k→∞

limk→∞ Rk

k→∞

−1 = limk→∞ Uk Uk−1

˜ k RX Dk RY R−1 D−(1+k) R−1 R ˜ −1 = limk→∞ R Y X k −1 . = RX DRX

Avˇsak

−1 RX DRX



λ1  0   =    0

× λ2

× ×

. . 0

a t´ım je d˚ ukaz vˇety 6.3 proveden. ||| 87

λN −1 0

 × ×       ×  λN

7

Poˇ c´ ateˇ cn´ı a okrajov´ e probl´ emy pro obyˇ cejn´ e diferenci´ aln´ı rovnice

7.1

Z´ akladn´ı poznatky z teorie obyˇ cejn´ ych diferenci´ aln´ıch rovnic

V t´eto kapitole budeme vyˇsetˇrovati metod pˇribliˇzn´eho ˇreˇsen´ı rovnic typu (7.1)

y 0 = f (x, y), x ∈ [a, b],

(7.2)

y(a) = y0 ,

kde −∞ < a < b < +∞ a soustav takov´ ych rovnic.  0  y1 = f1 (x, y1 , ..., yN ) .................................... , x ∈ [a, b] (7.3)  0 yN = fN (x, y1 , ..., yN ) a

yj = yj0 , j = 1, ..., N.

(7.4)

Poznamenejme, ˇze probl´em (7.3) - (7.4) zahrnuje jako speci´aln´ı u ´lohu z N = f (x, z, z 0 , ..., z N −1 ), x ∈ [a, b],

(7.5) s podm´ınkami (7.6)

z j−1 (a) = zj0 , j = 1, ..., N − 1.

Staˇc´ı poloˇzit z = y1 , z 0 = y2 , ..., z N −1 = yN , takˇze 0 y10 = y2 , y20 = y3 , ..., yN −2 = yN −1

a 0 yN = z (N ) = f (x, z, z 0 , ..., z N −1 ) = f (x, y1 , ..., yN ).

Kromˇe u ´loh typu (7.1) - (7.2) a (7.3) - (7.4) budeme vyˇsetˇrovati t´eˇz u ´lohy typu y 0 = f (x, y), avˇsak m´ısto poˇc´ateˇcn´ı podm´ınky (7.2) budeme poˇzadovati, aby platila rovnost (7.7)

r(y(a), y(b)) = 0,

kde r je vhodn´a funkce dvou promˇenn´ ych.

88

Obecnˇeji budeme vyˇsetˇrovat u ´lohy yj0 = fj (x, y1 , ..., yN ), j = 1, ..., N, s podm´ınkami (7.8) kde (7.9)

0 = r(y(a), y(b)),  r1 (u1 , ..., uN ; v1 , ..., vN ) r = r(u, v) =  ............................................  . rN (u1 , ..., uN ; v1 , ..., vN ) 

´ ´lohami okraUlohy typu (7.1) a (7.7), respektive (7.3) a (7.8) se naz´ yvaj´ı u jov´ymi na rozd´ıl od u ´loh poˇc´ ateˇcn´ıho typu (7.1) a (7.2) respektive (7.3) a (7.4). Numerick´e metody pˇribliˇzn´eho ˇreˇsen´ı v´ yˇse uveden´ ych u ´loh je zaloˇzeno na vlastnostech ˇreˇsen´ı tˇechto u ´loh. Zformulujme nˇekter´ a tvrzen´ı potˇrebn´ a pro dalˇs´ı vyˇsetˇrov´an´ı. Vˇ eta 7.1 Pˇredpokl´ any adejme, ˇze komponenty f1 , ..., fN vektoru f jsou definov´ a jsou spojit´e na oblasti ‰ ˆ S = (x, y) : a ≤ x ≤ b, y ∈ RN . D´ ale nechˇt existuje lipschitzovsk´ a konstanta L tak; ˇze

(7.10)

˜ )|| ≤ L||¯ ˜ ||, ¯ ) − f (x, y y−y ||f (x, y

plat´ı pro vˇsechna x ∈ [a, b] a y ∈ RN (Lipschitzova podm´ınka). avˇe jedna vektor - funkce Potom ke kaˇzd´emu x0 ∈ [a, b] a y0 ∈ RN existuje pr´ y(x) tak, ˇze 1) y je spojit´ a a spojitˇe diferencovateln´ a vektor - funkce v [a, b]; 2) y0 (x) = f (x, y(x)), x ∈ [a, b]; 3) y(x0 ) = y0 . Pozn´ amka. Jiˇz v´ıme, vze Lipschitzova podm´ınka je splnˇena, pakliˇze parci´aln´ı ∂f derivace ∂xj , j = 1, ..., N existuj´ı a jsou spojit´e na S. V praxi se m˚ uˇzeme setkat pˇrev´ aˇznˇe s pˇr´ıpady, kdy fj , j = 1, ..., N , jsou ∂f yt na S neohraniˇcen´e. V takov´em spojit´e na S, avˇsak derivace ∂xj mohou b´ pˇr´ıpadˇe m´a poˇc´ateˇcn´ı u ´loha (7.3) - (7.4) ˇreˇsen´ı, ale to nemus´ı b´ yt definov´ ano na cel´e oblasti. Pˇ r´ıklad 7.1 Nechˇt N = 1 a y 0 = y 2 , y(0) = 1, 0 ≤ x < +∞. 89

Vid´ıme, ˇze y(x) = a to je ohraniˇcen´a funkce jen pro x < 1.

1 1−x

Pro numerick´e poˇc´ıt´an´ı je d˚ uleˇzit´e, ˇze ˇreˇsen´ı z´avisej´ı na poˇc´ ateˇcn´ıch datech spojitˇe. Tvrzen´ı na toto t´ema je obsahem n´asleduj´ıc´ı vˇety. Vˇ eta 7.2 Za pˇredpoklad˚ u vˇety 7.1 plat´ı odhady (7.11)

||y(x, a1 ) − y(x, a2 )|| ≤ eL(x−a) ||a1 − a2 ||,

kde (7.12)

y0 (x, aj ) = f (x, y), y(a) = aj , j = 1, 2.

V´ yznamn´e pro teorii i praxi jsou line´arn´ı u ´lohy typu Y0 = T (x)Y, Y(a) = Y a ,

(7.13)

pˇriˇcemˇz T = T (x) je matice typu N × N . R˚ ust ˇreˇsen´ı je charakterizov´ an v n´asleduj´ıc´ı vˇetˇe. Vˇ eta 7.3 Nechˇt prvky matice T (x) jsou spojit´e funkce na [a, b] a nechˇt κ(x) = ||T (x)||.

(7.14) Potom plat´ı odhad (7.15)

7.2

||Y(x) − I|| ≤ exp

šZ

x a

› κ(t)dt − 1 pro x ≥ a.

Poˇ c´ ateˇ cn´ı u ´ lohy

Jako motivaci vyˇsetˇrujme skal´ arn´ı u ´lohu (7.16)

y 0 = f (x, y), y(x0 ) = y0 , x0 , x ∈ [0, 1], y ∈ R1 .

ˇ h > 0. Nahradme ˇ Bud v (7.16) derivaci diferenc´ı ˇc´ımˇz obdrˇz´ıme y(x + h) − y(x) = f (x, y(x)), h neboli (7.17)

y(x + h) = y(x) + hf (x, y(x)).

Obdrˇzeli jsme tak algoritmus pro ˇreˇsen´ı. 90

Algoritmus 7.1 Budˇ M cel´e kladn´e. Zvolme h=

x − x0 M

a poˇc´ıtejme rekurzivnˇe pomoc´ı formul´ı η0 = y0 , ηk+1 = ηk + hf (xk , ηk ), kde xk+1 = xk + h, k = 0, 1, ..., M − 1. Algoritmus 7.1 popisuje t. zv. Eulerovu polygon´ aln´ı metodu. Vˇs´ımnˇeme si, ˇze v Algoritmu 7.1 plat´ı vztahy η(x0 ; h) = y0 , a η(x0 + h; h) = η(x; h) + hf (x, η(x; h). To nav´ad´ı na myˇslenku definovati algortimus, v nˇemˇz roli f (x; η(x; h)) pˇrevezme vhodn´a funkce Φ = Φ(x, y; h, f ). Algoritmus 7.2 Volme h > 0 a nechˇt š η0 = y0 (7.18) ηk+1 = ηk + hΦ(xk , ηk ; h; f ) emprok = 1, 2, ... kde xk+1 = xk + h. Jiˇz v´ıme, ˇze pro Eulerovu polygon´aln´ı metodu jest Φ(x, y, h, f ) = f (x, y), Φ tedy nez´avis´ı na h. Pro dalˇs´ı u ´ˇcely oznaˇcme pro libovoln´e x ∈ [a, b] a y ∈ RN symbolem z = z(x) pˇresn´e ˇreˇsen´ı u ´lohy (7.19)

z0 (t) = f (x, z(t)), z(x) = y(x).

Definujme funkci ∆ pˇredpisem š z(x+h)−y h (7.20) ∆(x, y; h; f ) = f (x, y) 91

pro h > 0 pro h = 0.

V dalˇs´ım textu nebudeme vyznaˇcovati f v argumentu funkce ∆ a Φ. Pomoc´ı funkce ∆ definujeme t. zv. lok´ aln´ı diskretizaˇcn´ı chybu, tedy veliˇcinu τ (x, y, h) = ∆(x, y; h) − Φ(x, y; h).

(7.21)

V algoritmu 7.2 zaveden´e metody se naz´ yvaj´ı jednokrokov´e. To evidentnˇe proto, ˇze k z´ısk´an´ı hodnoty pˇribl´ıˇzen´ı v bodˇe xk se vyˇzaduje znalost hodnoty v pˇr´ısluˇsn´em bodˇe pˇredchoz´ım xk−1 . Jednokrokov´a metoda se ukazuje jako vhodn´a, jestliˇze lim τ (x, y; h) = 0,

(7.22)

h→0

coˇz, d´ıky tomu, ˇze lim ∆(x, y, h) = f (x, y),

(7.23)

h→0

vede k podm´ınce (7.24)

Φ(x, y, 0) = f (x, y).

ˇ ık´ame, ˇze jednokrokov´a metoda je konzistentn´ı, jestliˇze plat´ı relace (7.24). R´ Pozn´ amka. Eulerova polygon´aln´ı metoda je zˇrejmˇe konzistentn´ı. Pˇ r´ıklad 7.2 Pˇredpokl´ adejme, ˇze N = 1 a f m´ a vyˇsˇs´ı hladkost neˇz je pouh´ a ∂y jsou spojit´e na [a, b]. Nav´ıc nechˇt spojitost, ˇreknˇeme, ˇze ∂x 1 hp (p) z (x + θh), 0 < θ < 1. z(x + h) = z(x) + hz 0 (x) + h2 z 00 (x) + ... + 2 p! Protoˇze

d f (t, z(t)) |t=x dt = fx (t, z(t)) |t=x + fy (t, z(t)) |t=x z 00 (x) =

= fx (x, y) + fy (x, y)f (x, y), a

z 000 (x) = fxx (x, y) + 2fxy (x, y)[f (x, y)]2 + fy (x, y)z 00 (x),

zjist´ıme, ˇze (7.25)

 p  ∆(x, y; h) = z 0 (x) + h2 z 00 (x) + ... + hp! z (p) (x + θh) = f (x, y) + h2 [fx (, y) + fy (x, y)f (x, y)]  +... 92

Na pˇr. pro Eulerovu metodu τ (x, y; h) = ∆(x, y; h) − Φ(x, y; h)

h [fx (x, y) + fy (x, y)f (x, y)] + O(h). 2 a metoda je ˇr´ adu p, jestliˇze plat´ı Obecnˇe, ˇr´ık´ame, ˇze jednokrokov´ =

τ (x, y; h) = O(hp )

(7.26)

pro x ∈ [a, b], y ∈ RN a pro vˇsechny f ∈ [C (p) ([a, b])]N . Z pˇredchoz´ıho pˇr´ıkladu odvod´ıme pomˇernˇe obecn´ y zp˚ usob konstrukce metod vyˇsˇs´ıch ˇr´ad˚ u. Na pˇr. (7.27)

Φ(x, y; h) = f (x, y) +

h [fx (x, y) + fy (x, y)f (x, y)] 2

definuje metodu 2. ˇr´adu. Na tomto m´ıstˇe je nutn´e uv´ezt na pravou m´ıru skuteˇcnost, ˇze metody vysok´ ych ˇr´ad˚ u (> 2) jsou v praxi v´ıcem´enˇe bezcenn´e, ne´ uˇcinn´e, sloˇzit´e. To je zp˚ usobeno t´ım, ˇze v´ ypoˇcet vyˇsˇs´ıch derivac´ı je probl´em sloˇzit´ y v tom smyslu, ˇze je ˇspatnˇe podm´ınˇen´y (non well posed). Pˇripomeˇ nme, ˇze to znaˇc´ı, ˇze hodnota v´ ysledku v´ ypoˇcet nen´ı spojitˇe z´avisl´a na datech. ˇ Uvedme nˇekter´e pˇr´ıklady metod vyˇsˇs´ıch ˇr´ ad˚ u. Poloˇzme (7.28)

Φ(x, y; h) = a1 f (x, y) + a2 f (x + p1 h, y + p2 hy)),

kde a1 , a2 , p1 , p2 jsou vhodn´e konstanty, kter´e jest urˇcit tak, aby platily vztahy urˇcuj´ıc´ı ˇr´ad metody. Pro Φ z (7.28) plat´ı Φ(x, y; h) = (a1 + a2 )f (x, y) + a1 h [p1 fx (x, y) + p2 fy (x, y)f (x, y)] + O(hp ), odkud vypl´ yv´a, splnˇen´ı podm´ınek (7.29)

a1 + a2 = 1, a2 p1 = 1/2, a2 p2 = 1/2,

je postaˇcuj´ıc´ı aby metoda byla 2. ˇr´ adu. Speci´aln´i volbou a1 = a2 =

1 , p1 = p2 = 1 2 93

obdrˇz´ıme metodu Heunovu 1 [f (x, y) + f (x + h, y + hf (x, y))] . 2 Je patrn´e, ˇze tato metoda vyˇzaduje v´ ypoˇcet dvou funkˇcn´ıch hodnot na jeden krok mertody a to je cenou, jiˇz je nutno zaplatit za vyˇsˇs´ı ˇr´ ad. (7.30)

Φ(x, y; h) =

Jinou metodu obdrˇz´ıme volbou a1 = 0, a2 = 1, p1 = p2 = 1/2. Ta je d´ana sch´ematem Collatzov´ym ’ “ h h Φ(x, y; h) = f x + , y + f (x, y) . (7.31) 2 2

Toto sch´ema je opˇet ˇr´adu 2. a vyˇzaduje v´ ypoˇcet dvou funkˇcn´ıch hodnot na jeden krok metody. Dalˇs´ı metody vyˇsˇs´ı ˇr´ad˚ u jsou spojeny se jm´eny Runge a Kutta. Hledejme jenokrokovou metodu pomoc´ı sch´ematu

(7.32)

Φ(x, y; h) =

1 [k1 + 2k2 + 2k3 + k4 ] , 6

kde k1 k2 k3 k4

= f(x,  € y) , = f €x + h2 , y + h2 k1 , = f x + h2 , y + 12 hk3 , = f (x + h, y + hk3 ) . D´a se uk´azati, ˇze Φ(x, y; h) − ∆(x, y; h) = O(h4 ).

Znamen´a to, ˇze tato, d´a se ˇr´ıci, z´akladn´ı metoda Runge - Kuttovy tˇr´ıdy, je ˇr´adu 4. Aplikujm,e nˇekter´e jednokrokov´e metody na pˇr´ıpad u ´lohy y 0 (x) = f (x), y(x0 ). Pˇresn´e ˇreˇsen´ı t´eto u ´lohy je d´ano v´ yrazem Z x y(x) = f (t)dt + y0 . x0

94

Snadno nahl´edneme, ˇze v problematice numerick´e kvadratury, coˇz je vlastnˇe n´ami vyˇsetˇrovan´a u ´loha, Eulerova metoda vede na zn´amou metodu obd´eln´ıkovou zat´ımco Heunova metoda na metodu lichobˇeˇzn´ıkovou. Runge Kuttova metoda poskytuje metodu aproximace integr´ alu, jeˇz je zn´am´ a pod n´azvem metoda Simpsonova, kter´aˇzto metoda integruje polynomy aˇz do ˇr´adu 3. stupnˇe vˇcetnˇe. Ot´azku konvergence jednokrokov´ ych metod zodpov´ıd´ a n´asleduj´ıc´ı vˇeta. Soustˇred#ıme se na skal´arn´ı pˇr´ıpad s t´ım, ˇze odpov´ıdaj´ıc´ı vˇety maj´ı takˇrka shodn´e formulace i d˚ ukazy. Vˇ eta 7.4 Nechˇt f ∈ C([a, b]), x0 ∈ [a, b] a y0 ∈ R1 . Budˇ y = y(x) ˇreˇsen´ı poˇc´ ateˇcn´ı u ´lohy y 0 = f (x, y), y(x0 ) = y0 . Budˇ Φ funkce spojit´ a na G = {(x, y, h) : a ≤ x ≤ b, |y − y(x)| ≤ γ, 0 ≤ |h| ≤ h0 , , h0 > 0, γ > 0} a nechˇt existuj´ı konstanty κ1 , κ2 tak, ˇze |Φ(x, y1 ; h) − Φ(x, y2 ; h)| ≤ κ1 |y1 − y2 | a |τ (x, y(x); h)| = |∆(x, y(x); h) − Φ(x, y(x); h)| ≤ κ2 |h|p , p > 0, plat´ı pro vˇsechna x ∈ [a, b], |h| ≤ h0 . ˆ 0
exp{κ1 |x − x0 |} − 1 κ1

pro vˇsechna x ∈ [a, b] a vˇsechna hk =

x − x0 , k = 1, 2, ..., k

ˆ = h0 . ˆ Pro γ = ∞ je h pro nˇeˇz |hk | ≤ h. Obecnˇejˇs´ı metody pro ˇreˇsen´ı poˇc´ ateˇcn´ıch u ´loh lze shrnout pod obecn´e sch´ema (7.33) ηk+r + ar−1 ηk+r−1 + ... + a0 ηk = hF (xk , ηk+r , ηk+r−1 , ..., ηk ; h; f ). 95

Zde F = F (x, u1 , ...ur+1 ; h; f ) je vhodn´a funkce a a0 , ...ar−1 jsou parametry metody. Je zˇrejm´e, ˇze k nastartov´ an´ı metod typu (7.33) je nutn´e zn´at r poˇc´ateˇcn´ıch hodnot η0 , ..., ηr−1 . Teorie v´ıcekrokov´ych a obecnˇeji, metod typu (7.33), byˇt pˇrirozenˇe sloˇzitˇejˇs´ı neˇz ta, jiˇz jsme naznaˇcili pro metody jednokrokov´e, je do ˇs´ıˇre i do hloubky propracov´ana a dovedena do realizace v podobˇe softwarov´ ych bal´ık˚ u.

7.3

Okrajov´ eu ´ lohy

Elementy matematick´ e anal´ yzy slab´ ych ˇ reˇ sen´ı. V tomto ˇcl´anku se omez´ıme na line´arn´ı okrajov´e u ´lohy s pozn´amkou, ˇze neline´arn´ı u ´lohy lze v principu vyˇsetˇrovat a ˇreˇsit podobn´ ym zp˚ usobem jako v problematice line´arn´ı; rozd´ıl spoˇc´ıv´ a pouze v tom, ˇze v´ ysledn´e diskr´etn´ı soustavy jsou neline´arn´ı. V´ yznamn´ ym n´astrojem vyˇsetˇrov´ an´ı okrajov´ ych u ´loh je pojem slab´e (variaˇcn´ı) formulace okrajov´ ych u ´loh. ´lohy V´ yklad budeme prov´adˇet na pˇr´ıkladˇe typick´e, byˇt akademick´e, u (7.34)

− y 00 = f v Ω = [a, b],

(7.35)

µ1 y 0 (a) + µ0 y(a) = α a ν1 y 0 (b) + ν0 y(b) = β.

O funkci f pˇredpokl´ad´ame, ˇze patˇr´ı do L2 (a, b); prakticky vˇsak vˇzdy je f po ˇc´astech spojit´a, t. j. f m´a v [a, b] koneˇcn´ y poˇcet bod˚ u nespojitiosti. Form´alnˇe vyˇzaduje t. zv. variaˇcn´ı pˇr´ıstup pojem zobecnˇen´e derivace. Funkce g ∈ L2 (a, b), t. j. g, pro n´ıˇz Z b |g(t)|2 dt < +∞, a

a, ˇze m´a v [a, b] zobecnˇenou derivaci, pakliˇze plat´ı, ˇze existuje h ∈ L2 (a, b) takov´ Z b Z b (7.36) g(t)v 0 (t)dt = − h(t)v(t)dt + g(b)v(b) − g(a)v(a). a

a

´du) funkce g. Pak ˇr´ık´ame, ˇze h je zobecnˇenou derivac´ı (1. ˇra D´a se uk´azati, ˇze v naˇsem pˇr´ıpadˇe, kdy vyˇsetˇrujeme oblast Ω ⊂ R1 , funkce g m´a v [a, b] zobecnˇenou derivaci 1. ˇr´ adu pr´avˇe kdyˇz (klasick´ a) derivace g 0 existuje v [a, b] s vyj´ımkou mnoˇziny E ⊂ [a, b] maj´ıc´ı Lebesgueovu m´ıru m(E) = 0. Je zn´amo, ˇze spojit´a funkce g m´ a v [a, b] zobecnˇenou derivaci h pr´ avˇe kdyˇz g je absolutnˇe spojit´ a v [a, b]. Mnoˇzina funkc´ı maj´ıc´ı zobecnˇen´e derivace 1. (1) ˇr´adu integrovateln´e se ˇctvercem na Ω se znaˇc´ı symbolem W2 (Ω) a naz´ yv´ a se arn´ı souˇcin Sobolevov´ym prostorem. V tomto prostoru lze zav´ezti skal´ (7.37)

(u, v)W21 (Ω) = (u0 , v 0 )2 + (u, v)2 , 96

kde pro u, v ∈ L2 (Ω),

(u, v)2 =

Z

u(t)v(t)dt

Ω

a se, ˇze u0 , v 0 ∈ a u0 a v 0 znaˇc´ı zobecnˇen´e derivace funkc´ı u a v. Pˇredpokl´ad´ 2 L (Ω). Pˇredpokl´adejme, ˇze w ∈ L2 (Ω) je funkce splˇ nuj´ıc´ı okrajov´e podm´ınky µ0 w(a) + µ1 w0 (a) = α a ν0 w(b) + ν1 w0 (b) = β.

(7.38)

D´a se uk´azati, ˇze existuje hust´ a ve smyslu normy prostoru L2 (Ω) mnoˇzina tvoˇr´en´a dostateˇcnˇe hladk´ ymi funkcemi splˇ nuj´ıc´ımi (7.38); na pˇr. mnoˇzina vˇsech polynom˚ u splˇ nuj´ıc´ıch (7.38). Poloˇzme v (7.34) y(s) = u(s) − w(s). Jest tedy, −y 00 = w00 + f v Ω a µ1 y 0 (a) + µ0 y(a) = µ0 (u(a) − w(a)) + µ1 (u0 (a) − w0 (a)) = 0, ν1 y 0 (b) + ν0 y(b) = ν1 (u0 (b) − w0 (b)) + ν0 (u(b) − w(b)) = 0. Zn´ame-li y, zn´ame automaticky i hledna´e u. K urˇcen´ı y m´ame tedy vztahy (7.39)

− y 00 = f + w00 v Ω

(7.40)

µ0 y(a) + µ1 y 0 (a) = ν0 y(b) + ν1 y 0 (b) = 0.

Protoˇzee w00 je zn´am´a funkce, lze (7.39) ps´at ve tvaru (7.41)

− y 00 = F v Ω,

pˇri ˇcemˇz F je dan´a funkce. Vˇs´ımnˇeme si, ˇze ve formulaci (7.34) - (7.35), neboli (7.40) - (7.41), vyˇzadujeme splnˇen´ı tˇechto vztah˚ u identicky v [a, b]. To vˇsak vyˇzaduje znalost F ve vˇsech bodech oblasti Ω a obecnˇeji, znalost vˇsech dat v Ω. Toto je jednak ne pˇr´ıliˇs praktick´e a jednak dokonce ne vˇzdy splniteln´e. Ve v´ıcerozmˇern´ ych oblastech zad´an´ı okrajov´ ych podm´ınek na nˇekter´ ych ˇc´ astech hranice m˚ uˇze b´ yt nedefinovateln´e. Na pˇr. nechˇt Ω = [0, 1] × [0, 1] a nechˇt u(s, 0) = 0 pro s ∈ (0, 1), zat´ımco u(0, y) = 1 pro y ∈ (0, 1). Ot´azkou je, jakou hodnotu m˚ uˇze nab´ yvat u v bodˇe (1, 0), t. j. , ˇcemu je rovno u(1, 0). Uvˇedomme si, ˇze u m´a b´ yt spojit´a na Ω.

97

V takov´ ych probl´emech se uk´aˇze jako velice vhodn´a slab´ a (variaˇcn´ı) formulace probl´emu. Oznaˇcme symbolem V podprostor Sobolevova prostoru W21 (Ω) splˇ nuj´ıc´ıch okrajov´e podm´ınky (7.40). Vyn´asobme rovnici (7.41) funkc´ı v ∈ V a integrujme pod´el Ω. Obdrˇz´ıme vztahy Z b Z b 00 − y (t)v(t)dt = F (t)v(t)dt, v ∈ V, a

a

takˇze poˇzadovan´e ˇreˇsen´ı lze hledat pomoc´ı vztah˚ u (7.42)

Z

b 0

0

y (t)v (t)dt =

a

Z

b a

F (t)v(t)dt, v ∈ V,

coˇz plyne z (7.42) protoˇze v splˇ nuje tyt´eˇz okrajov´e podm´ınky (7.40) jako hledan´e ˇreˇsen´ı y (viz Cviˇcen´ı 7.3). ˇ L mnoˇzina vˇsech funkc´ı y maj´ıc´ıch spojit´e derivace 2. ˇr´ adu Pˇ r´ıklad 7.3 Bud nuj´ıc´ıch okrajov´e podm´ınky (7.40). Ukaˇzte, ˇze rovnosti na [a, b] a splˇ Z

b

y 00 (t)v(t)dt =

a

Z

b

u(t)v 00 (t)dt a

plat´ı pro vˇsechny prvky y ∈ L,pˇriˇcemˇz v je funkce maj´ıc´ı na [a, b] spojitou derivaci 2. ˇr´adu, pr´avˇe kdyˇz v ∈ L. Pro jednoduchost poloˇzme a = 0 a ,b = 1 a µ1 = ν1 = 0. V tomto pˇr´ıpadˇe se okrajov´e podm´ınky (7.40) redukuj´ı na (7.43)

y(0) = y(1) = 0

a tedy t´eˇz (7.44)

v(0) = v(1) = 0.

Formulujme nyn´ı naˇsi okrajovou u ´lohu variaˇcnˇe neboli slabˇe, ˇci pˇresnˇeji ˇreˇceno, ve slab´em smyslu. ˜ (1) (0, 1) tak, aby vztahy ´ Uloha (U). Nal´ezt y ∈ W 2 (7.45)

Z

1

y 0 (t)v 0 (t)dt =

Z

0

0

98

1

F (t)v(t)dt

˜ (1) (0, 1), kde W ˜ (1) (0, 1) znaˇc´ı uz´ platily pro vˇsechny prvky v ∈ W aˇer prostoru 2 2 ∞ C0 ([0, 1]) funkc´ı s kompaktn´ım nosiˇcem na [0, 1] maj´ıc´ıch na [0, 1] spojit´e derivace (1) vˇsech ˇr´ ad˚ u. Uz´ avˇer se m´ın´ı ve smyslu metriky prostoru W2 (0, 1). au ´loha D´a se uk´azat, ˇze u ´loha U m´a pr´avˇe jedno ˇreˇsen´ı. M´a-li okrajov´ y 00 = F rm v [0, 1] s podm´ınkami (7.44) klasick´e ˇreˇsen´ı, pak toto klasick´e ˇreˇsen´ı je totoˇzn´e s ˇreˇsen´ım u ´lohy U. Tato tvrzen´ı jsou zaloˇzena na platnosti slavn´eho Lax - Milgramova lemmatu Lemma 7.1 (Lax − Milgram) Budˇ V Hilbert˚ uv prostor. Budˇ B = B(u, v) biline´ arn´ı forma, jeˇz je spojit´ a a V −eliptick´ a, t.j. existuj´ı konstanty c > 0 a c˜ > 0 tak, ˇze |B(u, v)| ≤ c||u|| ||v||, u, v ∈ V (7.46) a (7.47)

|B(u, u)| ≥ c˜||u||2 , u ∈ V.

Potom pro kaˇzd´y prvek f ∈ V existuje pr´ avˇe jeden prvek u∗ ∈ V takov´y, ˇze B(u∗ , v) = (F, v), ∀v ∈ V.

(7.48) Plat´ı pˇritom (7.49)

||u∗ || ≤ κ||F ||,

avis´ı na F . kde κ nez´

˜ (1) (0, 1) a V naˇsem pˇr´ıpadˇe V = W 2 B(u, v) = Protoˇze 2

||u|| =

Z

1 0

Z

1

u0 (t)v 0 (t)dt. 0

‚ ƒ ||u||2 + ||u0 ||2 dt,

platnost poˇzadavku (7.47) je d˚ usledkem Friedrichsovy nerovnosti Z Z ˜ (1) (0, 1), u0 (t)2 dt ≥ τ u2 dt, ∀u ∈ W 2 Ω

Ω

a na u ani u0 . kde τ je kladn´a konstanta nez´avisl´ Poˇzadavek (7.46) lemmatu 7.1 je automaticky splnˇen d´ıky tomu, ˇze skal´ arn´ı ˜ 1 (0, 1) je d´an v´ souˇcin v W y razem (7.37). 2 Na z´akladˇe lemmatu 7.1 tedy plat´ı D˚ usledek 7.1 Existuje pr´ avˇe jedno slab´e ˇrevsen´ı u ´lohy (U). 99

Diskretizace. ˇ Vh ⊂ V = W ˜ 1 (0, 1), kde Bud 2 dimVh = N < +∞. Nechˇt tvoˇr´ı basi prostoru Vh , tedy

{φ1 , ..., φN } Vh = Span{φ1 , ..., φN }.

Tud´ıˇz, uh =

(7.50)

N X

ν k φk ,

k=1

pro kaˇzd´ y prvek uh ∈ Vh .

Diskretizac´ı u ´lohy (U) se rozum´ı u ´loha (Uh ) Nal´ezt u?h ∈ Vh takov´e, ˇze plat´ı (7.51)

B(u?h , vh ) = (f, vh ), ∀vh ∈ Vh .

Vˇ eta 7.5 Splˇ nuje-li forma B = B(u, v) poˇzadavky Laxova - Milgramova lemmatu, pak lze nal´ezt h0 > 0 tak, ˇze pro kaˇzd´e 0 < h ≤ h0 , existuje pr´ avˇe jedno ˇreˇsen´ı u?h u ´lohy (Uh ). Nav´ıc plat´ı odhad ku? − u?h kV ≤ κdh (u?h )kf k, avis´ı ani na h ani na f a kde κ nez´

dh (v) = inf {kv − vh kV : vh ∈ Vh } . Pˇ r´ıklad 7.4 Pro pˇr´ıpad, kdy v probl´emu (Uh ) h=

1 , a xk = kh, k = 0, 1, ..., N, N

a basi aproximaˇcn´ıho podprostoru tvoˇr´ı po ˇc´ astech line´ arn´ı funkce definovan´e formulemi  0 pro 0 ≤ x ≤ xk−1 ,    1 (x − (k − 1)h) pro xk−1 ≤ xk , h φk (x) =  1 − h1 (x − xk ) pro xk ≤ x ≤ xk+1 ,   0 pro xk+1 ≤ x ≤ 1, 100

lze odvodit, ˇze dh (u) = hkukV . ych prvk˚ u pro aproximaci D˚ usledkem je posl´eze odhad chyby v metodˇe koneˇcn´ u ´lohy (U) ve tvaru ku? − u?h |W ˜ 1 (0,1) ≤ κhkf kL2 (0,1) . 2

Vypl´yv´ a odtud speci´ alnˇe, ˇze ku? − u?h kL2 (0,1) ≤ κ ˆ h2 kf kL2 (0,1) . Po dosazen´ı vyj´adˇren´ı (7.50) do (7.51) obdrˇz´ıme vztahy N X

(7.52)

νk B(φk , φj ) = (f, φj ), j = 1, ..., N.

k=1

Vyplnˇen´ı podm´ınek Laxova - Milgramova lemmatu zaruˇcuje, ˇze det (B(φk , φj )) 6= 0 ? T a soustava (7.52) m´a tud´ıˇz pr´avˇe jedno ˇreˇsen´ı ν ? = (ν1? , ...νN ) . Matice (B(φk , φj ))

se d´ıky jej´ı interpretaci v mechanice kontinua naz´ yv´ a matic´ı tuhosti u ´lohy (Uh ). Je-li forma B = B(u, v) symetrick´ a a jsou-li splnˇeny podm´ınky Laxova Milgramova lemmatu, pak matice tuhosti u ´lohy (Uh ) je pozitivnˇe definitn´ı, t.j. (B(uh , uh )) ≥ τ (uh , uh ) ∀uh , uh ∈ Vh , τ > 0. V´ ypoˇcet prvk˚ u matice tuhosti spoˇc´ıv´ a ve stanoven´ı veliˇcin Z

a+ph

a+(p−1)h

φ0k (x)φ0j (x)dx, j, k = 1, ..., N, p = 1, ..., N − 1,

pˇri ˇcemˇz Ω = (a, b), −∞ < a < b < +∞. Pro v´ yˇse uveden´ y pˇr´ıpad jednorozmˇern´e oblasti s Ω a pro ˇc´ astech line´arn´ıch funkc´ı φk lze tyto integr´aly stanovit pˇresnˇe. Obecnˇe je nutn´e prvky matice tuhosti urˇcovat pˇribliˇznˇe pomoc´ı vhodn´ ych kvadraturn´ıch vzorc˚ u.

101