OSTRAVSKÁ UNIVERZITA PŘÍRODOVĚDECKÁ FAKULTA
NUMERICKÁ MATEMATIKA
ZUZANA VÁCLAVÍKOVÁ
OSTRAVA 2004
Na této stránce mohou být základní tirážní údaje o publikaci.
1
OBSAH PŘEDMĚTU
Úvod ........................................................................................................................................... 3 1. Úvodní pojmy....................................................................................................................... 5 1.1. Matematická a numerická úloha ................................................................................... 5 1.2. Zdroje chyb ................................................................................................................... 8 1.3. Podmíněnost úloh........................................................................................................ 11 1.4. Kontrolní otázky ......................................................................................................... 12 2. Řešení nelineárních rovnic ................................................................................................. 15 2.1. Formulace problému ................................................................................................... 15 2.2. Grafická metoda.......................................................................................................... 16 2.3. Iterační metody řešení nelineárních rovnic................................................................. 17 2.3.1. Metody startovací............................................................................................. 18 2.3.2. Metody zpřesňovací ......................................................................................... 25 2.4. Kontrolní otázky, korespondenční úkol ...................................................................... 32 3. Řešení soustav lineárních rovnic........................................................................................ 35 3.1. Formulace problému ................................................................................................... 36 3.2. Metody přímé .............................................................................................................. 37 3.3. Metody iterační ........................................................................................................... 46 3.4. Kontrolní otázky, korespondenční úkol ...................................................................... 49 4. Vlastní čísla a vlastní vektory matic .................................................................................. 53 4.1. Formulace problému .................................................................................................. 53 4.2. Částečný problém vlastních čísel................................................................................ 54 4.3. Úplný problém vlastních čísel..................................................................................... 57 4.4. Kontrolní otázky, korespondenční úkol ...................................................................... 62 Řešení úloh............................................................................................................................... 65
2 Závěr......................................................................................................................................... 71 Literatura .................................................................................................................................. 73
Úvod
3
ÚVOD Milá čtenářko, milý čtenáři, předkládaná výuková opora je oporou ke kurzu Numerická matematika, který je věnován numerickým metodám algebry. Numerická matematika je v současné době rozvoje počítačové techniky nevyhnutnou součástí nejenom samotné matematiky, ale i spousty dalších technických, či ekonomických oborů. Nabízí nám možnost řešit přibližnými metodami řadu problémů, které matematicky nelze řešit vůbec a nebo je to technicky obtížné. Tyto texty Vás seznámí se základními metodami pro řešení nelineárních rovnic, soustav lineárních rovnic a výpočtem vlastních čísel matic. V jednotlivých kapitolách je uváděn přibližný čas potřebný k prostudování učiva, pokuste se proto studovat systematicky, postupně a dodržujte termíny. V závěru každé kapitoly jsou Kontrolní úlohy, které by jste po prostudování měli být schopni vyřešit samostatně. Správnost řešení si můžete ověřit v kapitole Řešené úlohy. Přesto, že dnes již existuje mnoho programů, které zahrnují algoritmy popisovaných metod, doporučuji, aby jste si algoritmy zpracovali samostatně. Z praktického důvodu ponechávám pouze na Vás, jaký software k práci použijete. Kapitoly 2, 3 a 4 obsahují korespondenční úkol, který je nutno zaslat ke kontrole dle harmonogramu studia. Bližší informace o způsobu komunikace Vám sdělí tutor. Z časových důvodů obsahuje studijní opora pouze stručný popis těch nejzákladnějších numerických metod. Podrobnější informace o jednotlivých metodách, jakož i popis dalších metod můžete naleznout v doporučené literatuře, která je uvedena v závěru textů. Přeji Vám mnoho úspěchů ve studiu. Zuzana Václavíková Autorka studijní opory
Čas potřebný k prostudování učiva předmětu: 30 + 30 hodin (teorie + řešení úloh)
4
1. Úvodní pojmy
1. ÚVODNÍ POJMY V této kapitole se dozvíte: •
co to je numerická matematika a co jsou numerické úlohy;
•
co jsou zdroje chyb a jak pracovat s přibližnými čísly;
•
co je podmíněnost úloh a algoritmů.
Budete schopni: •
transformovat reálný problém, případně matematickou úlohu na úlohu numerickou;
•
odhadnou, jaké chyby budou mít vliv na kvalitu numerického řešení a určit chybové zatížení výsledku;
•
určit, zda je takováto úloha dobře/špatně podmíněná, a tedy zda ji numericky lze, nebo nelze úspěšně řešit.
Klíčová slova této kapitoly: numerická úloha, algoritmus, chyba, zdroje chyb, podmíněnost úlohy, stabilita algoritmu . Čas potřebný k prostudování učiva kapitoly: 4 + 4 hodiny (teorie + řešení úloh)
Průvodce studiem. První kapitola je úvodem do numerické matematiky. Ukážeme si jaké úlohy lze numericky řešit, co všechno obnáší výpočet přibližnými metodami a jak pracovat s nepřesnými čísly, abychom dosáhli co nejkvalitnějších výsledků.
1.1 Matematická a numerická úloha Než si vysvětlíme, co vlastně numerická matematika je, zkusme se zamyslet nad tím, jak se dá dospět od skutečného reálného problému až k jeho řešení. V technické praxi tato cesta může být dvojí: experimentální nebo výpočetní. Experimentálním řešením rozumíme proces, kdy necháme studovaný děj proběhnout a na základě výsledku zodpovíme otázky, které nás zajímaly. Řešení výpočetní znamená, že na základě veličin, které mají na daný proces vliv a funkčních vztahů mezi nimi, matematicky vypočteme jak bude děj probíhat. Oba způsoby však mají svoje výhody i nevýhody. Experiment nám může výborně sloužit k ověření našeho „výpočetního řešení“, ale často se provést nedá. Například, pokud nás zajímá, jaká teplota byla v době ledové, nezbývá nám nic jiného, než ji výpočetně odhadnout, protože čas vrátit nemůžeme. Samozřejmě je
5
6 třeba mít neustále na vědomí, že „výpočetní řešení“ v konečném důsledku může, ale nemusí odpovídat skutečnosti, protože během výpočtu se dopouštíme spousty zaokrouhlení a přiblížení. Uvažujme již zmíněný reálný problém: Jaká byla teplota v době ledové? Aby byl problém výpočetně řešitelný, musí být převeden na tzv. matematickou úlohu- což znamená převedení „slovní úlohy“ do matematické podoby. V našem případě vědec (fyzik, geolog, paleontolog,…) jednoznačně určí, co mělo na teplotu vliv a na základě těchto znalostí sestaví rovnice (nebo jiné matematické vztahy), jejichž řešením obdržíme hodnotu teploty. Jsou-li tyto rovnice matematicky řešitelné, můžeme je vyřešit a dostaneme tak teoretické přesné řešení matematické úlohy. Často se však stane, že rovnice řešit buďto nelze vůbec, nebo je to časově zdlouhavé, náročné, atp. Potom nám nezbývá nic jiného, než se uspokojit alespoň s přibližným řešením matematické úlohy. Matematickou úlohu tak převedeme na úlohu numerickou tj. na jednoznačný popis funkčního vztahu mezi konečným počtem vstupních a výstupních dat (resp. dat vyjádřitelných konečným počtem čísel). Numerická úloha je tedy matematický model reálného problému, který lze realizovat na počítači v konečném čase. Jednoznačná specifikace konečné posloupnosti operací, pomocí kterých ze vstupních dat jednoznačně obdržíme výstupní data nazýváme algoritmem. Matematická disciplína, která konstruuje a analyzuje metody, algoritmy a postupy pro realizaci numerických úloh na počítači se nazývá numerická matematika Příklad. Přetransformujte problém na úlohu matematickou a poté na úlohu numerickou: Nalezněte řešení diferenciální rovnice y′′ − 4 y′ + 4 y = 0, y (0) = 1, y′(0 ) = 4 . Řešení: Reálný problém je již ve tvaru matematické úlohy, neboť máme přímo rovnici (lineární diferenciální rovnici 2. řádu s konstantními koeficienty homogenní) jejíž řešením obdržíme odpověď. Úloha však není numerická, protože výstupní data-tj. samotné řešení rovnice, je funkce y ( x ) , která není vyjádřitelná konečným počtem čísel (představuje nespočetně mnoha údajů- funkční hodnoty v každém reálném čísle, a těch je nespočetně.) Abychom úlohu převedli na numerickou, musíme z nespočetně mnoho výstupních dat „vybrat“ konečně mnoho, například tak, že budeme požadovat řešení diferenciální rovnice pouze v bodech x = 0, x = 0.1, x = 0.2, ..., x = 0.9, x = 1 , tedy výstupem bude 11 hodnot, konkrétně hodnoty y (0 ), y (0.1), y (0.2 ),..., y (0.9 ), y (1) . Příklad. Přetransformujte problém na úlohu matematickou a poté na úlohu numerickou: Nepřímým měřením byla zjišťována napínací síla lana. Experiment byl organizován tak, že mezi dvěma body A, B upevněnými ve vzdálenosti c bylo zavěšeno lano délky d . Na toto lano bylo zavěšeno závaží o hmotnosti m tak,
1. Úvodní pojmy aby se nehýbalo. V místě zavěšení závaží se vytvořil pravý úhel. Zjistěte napínací sílu delší části lana. Řešení: Experiment si můžeme načrtnout takto
A
B
c
α x
d-x 90˚
r FB
X α r FA
r r Fg = mg V místě zavěšení se vytvořil pravý úhel, tedy víme, že úhel AXB je 90o , dále označme úhel ABX jako α . Z fyziky víme, že hmotnost závaží vyvolá tíhovou r r r sílu Fg = mg , kde g je tíhové zrychlení, a tato síla se rozkládá na dvě složky r napínacích sil: jedna ve směru od bodu B , označili jsme ji FB a druhá ve směru r od bodu A , označená FA , přičemž se zachová úhel α (viz obrázek). Přeponu pravoúhlého trojúhelníku ABX jsme označili c , což je vzdálenost bodů A, B a odvěsny x, d − x , neboť víme pouze to, že součet těchto stran je d - délka lana. r Zajímá nás velikost síly FB , víme, že r FB = mg sin α a úhel α můžeme spočítat z trojúhelníku ABX jako x sin α = . c Z Pythagorovy věty určíme vzdálenost x , protože x 2 + (d − x )2 = c 2 ,
tedy d ± 2c 2 − d 2 . 2 My jsme si označili x kratší odvěsnu v trojúhelníku (uvědomme si, že značení mohlo být i obráceně), tedy pro nás je vzdálenost x menší ze dvou kořenů kvadratické rovnice, tudíž x=
7
8
d − 2c 2 − d 2 . 2 Nyní již vidíme, že tato úloha vede na matematickou úlohu spočítat r d − 2c 2 − d 2 . FB = mg 2c Matematická úloha je zároveň i numerickou úlohou. Vstupní data jsou čtyři r hodnoty m, g , c, d a výstupní hodnota je jedno číslo FB . x=
1.2 Zdroje chyb.
Při vytváření matematického modelu reálného problému, stejně tak i při samotném řešení numerickou cestou se vždy dopouštíme jisté „idealizace“. Při sestavování matematické úlohy jsme často například nuceni použít fyzikální vztahy, které platí pouze za „ideálních podmínek“- třeba zanedbáváme tření, atp., nebo nejsme schopni zachytit všechny faktory, které mají na daný děj vliv. Rozdíl řešení reálného problému a řešení matematické úlohy nazýváme chybou matematického modelu. Chyba, které se dopouštíme při samotném numerickém řešeni nazýváme chybou metody. Je nutné si uvědomit, že pokud je chyba matematického modelu příliš velká (model špatně popisuje reálnou situaci) nebude výsledek odpovídat skutečnosti, bez ohledu na přesnost samotného řešení. K posouzení kvality výsledku musíme vzít v úvahu i tzv. chyby ve vstupních datech, protože ty jsou často výsledkem měření a každé měření je zatíženo chybou. A nakonec jsou to chyby zaokrouhlovací, způsobené realizací výpočtu na počítači. Označme x skutečnou hodnotu nějakého čísla a ~ x jeho přibližnou hodnotu, tzv. aproximaci. Chybu lze potom vyjádřit dvojím způsobem, jako • absolutní chybu A( x ) = x − ~ x , která vystihuje odchylku, co do velikosti a ~ x−x • relativní chybu R( x ) = ~ , určující odchylku procentuelně. x Příklad. Vezměme číslo x = 2 = 1.41421... a jeho aproximaci ~ x = 1.41 . Určete absolutní a relativní chybu aproximace.
Řešení: Absolutní chyba je A( x ) = x − ~ x = Relativní chyba je
2 − 1.41 = 1.41421... − 1.41 = 0.00421...
x−~ x A( x ) 0.00421... = = = 0.002988... , ~ ~ x x 1.41 tedy odchylka aproximace od skutečné hodnoty činí přibližně 0.2988 % . R( x ) =
1. Úvodní pojmy
Poznámka.
1 ⋅ 10− d , 2 kde d je počet desetinných míst na která zaokrouhlujeme. Tedy například chceme-li zaokrouhlit na dvě desetinná místa, tak číslo 2.984 ≈ 2.98 , ale číslo 2.985 ≈ 2.99 . Zaokrouhlování provádíme tak, aby zaokrouhlovací chyba byla A( x ) ≤
Šíření chyb. Budeme-li pracovat s čísly přibližnými, naskýtá se otázka, jak se chyba aproximace bude během výpočtu kumulovat. Platí následující věta.
Věta. Nechť x1 , x2 jsou přesné hodnoty a ~ x1 , ~ x2 jejich aproximace. Potom platí a) A( x1 ± x2 ) = A( x1 ) + A( x2 ) b) A( x1 ⋅ x2 ) = A( x1 ) ⋅ ~ x2 + A( x2 ) ⋅ ~ x1 x2 + A( x2 ) ⋅ ~ x1 ⎛ x1 ⎞ A( x1 ) ⋅ ~ ⎟ ⎜ . c) A⎜ ⎟ = ~ x22 ⎝ x2 ⎠ Příklad. Obě uvedená čísla jsou zaokrouhlena správně na uvedený počet desetinných míst. Určete nejmenší interval, v němž bude ležet součet, rozdíl, součin a podíl těchto čísel. ~ x = 7.9256, ~ y = 3.29 .
Řešení: Je-li číslo ~ x = 7.9256 správně zaokrouhleno na čtyři desetinná místa, znamená to, že zaokrouhlovací chyba je A( x ) ≤ 0.5 ⋅ 10−4 = 0.00005 a podobně A( y ) = 0.5 ⋅ 10 −2 = 0.005 . Tedy chyba součtu a rozdílu těchto čísel dle předchozí věty bude A( x ± y ) = A( x ) + A( y ) = 0.00505 , tudíž nejmenší interval, ve kterém bude ležet součet těchto čísel je [(x + y ) − A(x + y ); (x + y ) + A(x + y )] tedy [11.2156 − 0.00505; 11.2156 + 0.00505] = [11.21055; 11.22065] . Nejmenší interval, ve kterém bude ležet rozdíl čísel je [(x − y ) − A(x − y ); (x − y ) + A(x − y )] , tedy [4.6356 − 0.00505; 4.6356 + 0.00505] = [4.63055; 4.64065]
9
10
Analogicky postupujeme u součinu a podílu. Chyba součinu je A( x ⋅ y ) = A( x ) ⋅ ~ y + A( y ) ⋅ ~ x = 0.0397925 a interval, v němž
leží součin bude [( x ⋅ y ) − A( x ⋅ y ); ( x ⋅ y ) + A( x ⋅ y )] , tedy [26.0354315; 26.1150165] . y + A( y ) ⋅ ~ x ⎛ x ⎞ A( x ) ⋅ ~ Chyba podílu je A⎜⎜ ⎟⎟ = = 0.00367629 a nejmenší interval, 2 ~ y ⎝ y⎠ ve kterém bude ležet podíl je [2.40532067; 2.41267325] .
Chyba při výpočtu funkčních hodnot. V praxi často potřebujeme určit chybu při výpočtu funkční hodnoty. K odvození použijeme Taylorova rozvoje. r r Nechť f X je funkce n -proměnných X = ( x1, x2 ,..., xn ) a předpokládejme, že r místo přesných hodnot X = ( x1, x2 ,..., xn ) pracujeme s aproximacemi ~r ~r X = (~ x 1, ~ x2 ,..., ~ xn ) . Taylorův rozvoj kolem bodu X je n ~r r 1 n ∂ 2 f ~r ∂f ~r xi ) + ∑ 2 ⎛⎜ X ⎞⎟ ⋅ ( xi − ~ xi )2 + ... f X = f ⎛⎜ X ⎞⎟ + ∑ ⎛⎜ X ⎞⎟ ⋅ ( xi − ~ ⎝ ⎠ i =1 ∂xi ⎝ ⎠ 2 i =1 ∂xi ⎝ ⎠ a zanedbáme-li členy druhého a vyšších řádu dostaneme, že chyba funkční hodnoty je přibližně n ~r r r ∂f ⎛ ⎞ ⋅ A( xi ) . A f X ≈ f X − f⎜X⎟ = ∑ ⎝ ⎠ i =1 ∂xi
( )
( )
( ( ))
( )
Příklad. Určete velikost chyby napínací síly delší části lana (viz příklad z předchozí části), pokud víte, že hmotnost závaží je m = (2.00 ± 0.03)kg , vzdálenost bodů A, B je c = (21.50 ± 0.05)cm , délka lana d = (280.3 ± 0.5)mm a tíhové zrychlení
dosazujeme zaokrouhleno na dvě desetinná místa jako g = 9.81 ms −2 . Řešení: V předchozím příkladě jsme ukázali, že napínací síla delší části lana je r d − 2c 2 − d 2 . FB = mg 2c Hodnoty vstupních veličin jsou zatížené chybami, připomeňme, že zápis m = (2.00 ± 0.03)kg znamená, že hmotnost závaží je m = 2.00 kg a chyba s jakou jsme hmotnost naměřili je A(m ) = 0.03 kg . Hodnoty c, d musíme ještě převést na základní jednotky- tedy na metry, takže c = 0.2150 m , A(c ) = 0.0005 m , d = 0.2803 m a A(d ) = 0.0005 m . Tíhové zrychlení je zaokrouhleno na dvě
desetinná místa g = 9.81 a tak A( g ) = 0.5 ⋅ 10−2 . r r Hodnota FB je funkcí čtyř proměnných m, g , c, d , tedy velikost chyby FB lze spočítat jako
1. Úvodní pojmy
11
r r r r ∂ FB ∂ FB ∂ FB ∂ FB r ⋅ A(d ) . ⋅ A(c ) + ⋅ A( g ) + ⋅ A(m ) + A FB = ∂d ∂c ∂g ∂m
( )
Po spočtení všech příslušných parciálních derivací a dosazení do tohoto vztahu obdržíme
( )
r d − 2c 2 − d 2 d − 2c 2 − d 2 ⋅ A(m ) + m ⋅ A( g ) + A FB = g 2c 2c ⎧ ⎪ 2c 2 2c 2 − d 2 ⎪⎪ + − mg ⎨ ⎪ ⎪ ⎪⎩
(
)
−
1 2
⎡ + ⎢d − 2c 2 − d 2 ⎢⎣ 2c 2
(
mg ⎡ + ⎢1 + d 2c 2 − d 2 2c ⎢ ⎣
(
)
−
1⎤ 2⎥
⎥⎦
1 ⎤⎫ 2⎥⎪
)
⎥⎦ ⎪⎪ ⎬ ⋅ A(c ) + ⎪ ⎪ ⎪⎭
⋅ A(d )
r Dosazením konkrétních čísel a vyčíslením dostaneme FB = 7.414 N a r r A FB = 0.293 N , což znamená, že relativní chyba je R FB = 0.0395 tedy asi
( )
( )
3.95% . V praxi se měření, jehož chyba je do 5% , považuje za dobré.
1.3 Podmíněnost úloh.
Ukázali jsme si, že pokud jsou vstupní data zatížena chybou, potom tato chyba nejenom že „nezmizí“ během výpočtu, ale naopak, může se ještě znásobovat. Je tedy na místě ptát se jak moc se během řešení může chyba kumulovat, neboli jak velkou chybu na výstupu způsobí vstupní chyby. Budeme říkat, že úloha je dobře podmíněná jestliže malá změna ve vstupních datech vyvolá malou změnu na výstupu. Míru této změny vyjadřujeme číslem podmíněnosti c p , které je rovno podílu relativní chyby na výstupu a relativní chyby na vstupu. r r Označíme-li vstupní hodnoty vektorem X arvýstupní Y , potom RY r . cp = RX Čím je číslo podmíněnosti větší, tím je úloha hůře podmíněná ( a tudíž i malá změna na vstupu způsobí velké chybové zatížení výstupu), pro c p ≥ 100 mluvíme o velmi špatně podmíněných úlohách. U takovýchto úloh je vhodné zvážit, zda má vůbec smysl úlohu řešit, nebo zda není lepší pokusit se problém přeformulovat na jinou matematickou či numerickou úlohu. Algoritmus nazveme stabilní je-li dobře podmíněný a numericky stabilní (tedy málo citlivý na zaokrouhlovací chyby během realizace na počítači).
() ( )
12 Příklad. Určete, zda úloha stanovit tan x v okolí bodu 1.57 je dobře podmíněná.
Řešení: Vezměme x = 1.57 a A( x ) = 0.01 . Tato chyba na vstupu vyvolá na výstupu chybu A( y ) = tan x − tan ( x + A( x )) = tan 1.57 − tan 1.58 = 1364.41 . Relativní chyba na
vstupu je R(x ) =
A( x ) = 0.0064 x
a relativní chyba na výstupu je
A( y ) = 1.0865 y a tedy číslo podmíněnosti c p = 170.58 . Úloha stanovit tan x v okolí bodu 1.57 je R( y ) =
velmi špatně podmíněná. Souvisí to s faktem, že v bodě x =
π 2
≈ 1.571 funkce
tan x není definovaná.
1.4 Kontrolní otázky. Kontrolní otázky.
1. Rozhodněte, zda úloha je matematická, resp. numerická. a) Určete všechny kořeny rovnice x5 − 2 x 3 + 4 x 2 − 3x + 2 = 0 . b) Spočítejte neurčitý integrál ∫ cos x dx .
Určete, co jsou vstupní a výstupní hodnoty. 2. Pokuste se popsat nějaký problém ve vaší oblasti studia, převeďte ho na úlohu matematickou a poté na úlohu numerickou. 3. Chceme zjistit, z čeho je vyrobena koule, a to tak, že změříme hustotu materiálu, a poté určíme z tabulek, o jakou látku se může jednat. Hustotu měříme nepřímo pomocí hmotnosti a objemu, protože platí m ρ= , V kde m je hmotnost koule, V je její objem. Hmotnost a průměr koule d můžeme změřit přímo. Převeďte na úlohu na numerickou. Určete, co jsou d vstupní a výstupní hodnoty. (Pomůcka. Objem koule s poloměrem r = je 2 4 V = πr 3 .) 3 4. Předpokládejte, že uvedená čísla jsou správně zaokrouhlena na daný počet desetinných míst. Určete absolutní a relativní chybu a nejmenší interval, ve kterém bude ležet výsledek
1. Úvodní pojmy
13
a) 1.5283+12.7391 b) 2.383-2.384 c) 13.635*2.354 d) 10.5394/0.0025 Povšimněte si nárůst relativní chyby v případě odečítání téměř stejných čísel a nárůst absolutní chyby při dělení malým číslem (číslem s nulami za desetinnou čárkou). Těmto situacím se proto v numerických výpočtech snažíme vyhnout. 5. Spočítejte chybu měření (absolutní i relativní), které jsme popsali v příkladě 3, jestliže průměr koule jsme naměřili d = (5.00 ± 0.05)cm a hmotnost m = (0.5845 ± 0.0005)kg . 6. Rozhodněte, zda úloha určit cos x v okolí bodu
π 2
je dobře podmíněná.
Spočtěte číslo podmíněnosti. 7. Určete, zda numerická úloha určit hustotu koule pomocí hmotnosti a průměru (dle příkladu 3) je dobře podmíněná, vzhledem na změny měření a) hmotnosti b) průměru koule.
Shrnutí kapitoly. Abychom mohli úspěšně používat metody numerické matematiky k výpočtům konkrétních problému, je nutno mít pořád na paměti rizika, která to sebou nese. Je proto nevyhnutné pochopit, co je to přibližná hodnota (aproximace) a její chyba a jak se během výpočtu chovají. Rovněž je potřebné dokázat transformovat úlohu tak, aby byla dobře podmíněná a tudíž numerické řešení vedlo ke „správným“ výsledkům. Pokud jste byli schopni samostatně vyřešit kontrolní otázky, jste na nejlepší cestě k pochopení dalších kapitol.
2.Řešení nelineárních rovnic
2. ŘEŠENÍ NELINEÁRNÍCH ROVNIC V této kapitole se dozvíte:
•
co je iterační proces;
•
jak numericky naleznout kořeny nelineárních rovnic, které analyticky řešit nelze, nebo lze jen velice obtížně;
•
co jsou metody startovací a zpřesňovací a jak je vhodně kombinovat, abychom dosáhli co nejkvalitnějšího numerického řešení.
Budete schopni:
•
úplně numericky vyřešit nelineární rovnici (metodou grafickou, startovací i zpřesňovací);
•
odhadnou chybu řešení;
•
zhodnotit kvalitu numerického řešení a případně navrhnout možnosti zlepšení výpočtu.
Klíčová slova této kapitoly: iterace, iterační proces, konvergence iteračního procesu, rychlost konvergence
Čas potřebný k prostudování učiva kapitoly: 6 + 6 hodin (teorie + řešení úloh)
Průvodce studiem. Většinu nelineárních rovnic nelze řešit jinak než numericky. V této kapitole popíšeme nejzákladnější metody numerického řešení a načrtneme, jak dané metody vhodně kombinovat k dosažení co nejlepšího řešení.
2.1 Formulace problému.
Nechť funkce f : R → R je definována a spojitá na intervalu [a, b] . Hledáme tzv. kořen rovnice f ( x ) = 0 , tedy reálné číslo α ∈ [a, b] , takové, že f (α ) = 0 .
15
16 Poznámka. Pro jednoduchost se budeme zabývat hledáním kořene jedné nelineární rovnice, avšak některé uvedené metody můžeme zobecnit pro soustavu n -obecně r nelineárních rovnic. Uvažujme zobrazení F : R n → R n definováno v oblasti r Ω ⊂ R n . Označíme-li složky funkce F = (F1,..., Fn ) , můžeme soustavu zapsat ve tvaru F1( x1,..., xn ) = 0
F2 ( x1,..., xn ) = 0
... Fn ( x1,..., xn ) = 0 r r r r Hledáme n -tici α = (α1,...,α n ) ∈ Ω , pro kterou F (α ) = 0 , neboli F1 (α1,...,α n ) = 0 F2 (α1,...,α n ) = 0
... Fn (α1,...,α n ) = 0
2.2 Grafická metoda.
Řešíme-li rovnici f ( x ) = 0 je vhodné nejdříve užít tzv. grafické metody řešení, neboli obyčejné načrtnutí grafu, které nám dá jistou představu o počtu kořenů rovnice a jejich lokalizaci.
Příklad. Nalezněte nejmenší kladný kořen rovnice e x − 2 cos x = 0 .
Řešení: Grafické řešení můžeme provést dvojím způsobem. Buďto si načrtneme přímo graf funkce y = e x − 2 cos x , kořeny rovnice budou průsečíky grafu funkce a osy x (graf 1), nebo načrtneme do téhož grafu dvě funkce y1 = e x , y2 = 2 cos x a kořeny budou průsečíky grafů funkcí y1 a y2 (graf 2).
2.Řešení nelineárních rovnic
17
Graf 1. exp(x) - 2*cos(x) y 1.5
1
0.5 0 -1
-0.5
0
0.5
1 x
-0.5
-1
Graf 2. exp(x), 2*cos(x) y 2.5
2
1.5
1
0.5 -1
-0.5
0
0.5
1 x
Červenou barvou je načrtnutý graf funkce y1 = e x , modře graf y2 = 2 cos x . Z obou grafů je vidět, že hledaný kořen rovnice leží v intervalu [0, 1] . 2.3 Iterační metody řešení nelineárních rovnic.
U numerických metod, kterým se budeme v této kapitole věnovat, budeme postupovat tak, že skutečné řešení hledáme jako limitu vhodně zvolené posloupnosti přibližných aproximací řešení. Takovým metodám říkáme iterační. Jednotlivé prvky této posloupnosti nazýváme iterace a nultý prvek iterační posloupnosti nazýváme počátečná aproximace, nebo počátečné přiblížení.
18
U iteračních metod nás budou zajímat především dvě otázky: 1. Zda námi sestavená „vhodná“ posloupnost jednotlivých iterací skutečně konverguje k hledanému řešení 2. a pokud ano, jak rychle konverguje. V praxi se dělíme iterační metody řešení nelineárních rovnic do dvou skupin, na • metody startovací a • metody zpřesňovací. Metody startovací jsou za daných předpokladů vždy konvergentní, nezávisle na „kvalitě“ počátečného přiblížení, ale konvergence je zpravidla velice pomalá. Naopak metody zpřesňovací jsou podstatně rychlejší, ale velmi citlivé na počátečnou aproximaci- ke konvergenci vyžadují, aby počátečná aproximace dostatečně kvalitně „přibližovala“ hledané řešení. Chceme-li využít přednosti obou typů metod, je rozumné je při řešení nelineární rovnice vhodně kombinovat- nejdříve užít metody startovací, pomocí ní určit dostatečně kvalitní přiblížení hledaného kořene a poté ho „zpřesňovat“ metodou zpřesňovací dokud nebude chyba menší než předem zadaná tolerance. Poznámka. V praxi máme vždy představu, jak kvalitně si přejeme výsledek zjistit, neboli jaké maximální možné chyby se můžeme dopustit, aby pro měl nás výsledek ještě smysl. Mluvíme o předem zadané toleranci. U řešení nelineárních rovnic může být náš požadavek dvojího typu: • buď požadujeme, aby α − α~ < ε , tedy aby se aproximace kořene od
•
skutečného kořene lišila maximálně o malou hodnotu ε , nebo požadujeme, aby f (α~ ) < δ , tedy povolíme, že funkční hodnota v aproximaci kořene může být místo skutečné nuly (kořen rovnice) menší než nějaké malé číslo δ .
2.3.1 Metody startovací. Metoda bisekce. Nejjednodušší z metod startovacích je tzv. metoda půlení intervalu, neboli metoda bisekce. Předpokládejme, že hledáme kořen reálné funkce f spojité na intervalu I 0 = [a0 ,b0 ] a navíc předpokládejme, že funkční hodnoty v koncových bodech intervalu mají opačná znaménka, tedy f (a0 ) ⋅ f (b0 ) < 0 . Tyto podmínky nám zaručují, že v intervalu I 0 = [a0 ,b0 ] existuje alespoň jeden kořen, tedy alespoň jedno α ∈ [a0 ,b0 ] pro které f (α ) = 0 . (Spojíme-li bod nad osou x s bodem pod osou x spojitou čarou, musíme nutně projít osou x , tedy bodem α - viz obrázek.)
2.Řešení nelineárních rovnic
19
y = f (x )
a0
α
b0
a0 + b0 , rozpůlí interval na dva další 2 [a0 , s0 ] a [s0 ,b0 ] , přičemž hledaný kořen rovnice bude ležet v tom z nich, pro který mají funkční hodnoty v krajních bodech opět opačná znaménka.
Střed intervalu I 0 = [a0 ,b0 ] , bod s0 =
s0 = b1 a0 = a1
α
Tímto postupem vytvoříme posloupnost intervalů
b0
{I k }∞k =0 ,
kde I k je ten
z intervalů [ak −1, sk −1 ] , [sk −1, bk −1 ] v jehož koncových bodech nabývá funkce f opačná znaménka. Kořen α leží v každém z intervalů I k a vzhledem k faktu, že délka intervalů se zmenšuje vždy o polovinu, budeme se k číslu α konvergentně přibližovat. Délka k -tého intervalu I k je 1 1 1 bk − ak = (bk −1 − ak −1 ) = 2 (bk − 2 − ak − 2 ) = ... = k (b0 − a0 ) 2 2 2
20
Můžeme tedy udělat odhad chyby (neboli určit hodnotu maximální možné chyby, b −a které jsme se mohli dopustit). Pro odhad chyby platí ak − α ≤ 0 k 0 2 b0 − a0 . respektive bk − α ≤ 2k Algoritmus metody bisekce: Nechť ε ≥ 0 je maximální možná chyba, se kterou chceme určit řešení nelineární rovnice f ( x ) = 0 . Nechť interval I 0 = [a0 ,b0 ] je takový, že f (a0 ) ⋅ f (b0 ) < 0 . Vstup: ε , a := a0 , b := b0 Opakuj a+b s := 2 a := s jestliže f (a ) ⋅ f (s ) > 0 b := s jestliže f (a ) ⋅ f (s ) < 0 konec jestliže f (a ) ⋅ f (s ) = 0 dokud b − a ≤ ε .
Výstup: I k = [a, b] .
Poznamenejme ještě, že metoda bisekce je sice pomalá, ale konverguje vždy, nezávisle na vlastnostech funkce f .
Příklad. Určete metodou bisekce nejmenší kladný kořen rovnice e x − 2 cos x = 0 s chybou menší než ε = 0.02 . Jako počáteční aproximaci použijte interval zjištěný grafickou metodou I 0 = [0; 1] , pracujte na 6 desetinných míst.
Řešení: Nejdříve ověříme, zda lze metodu bisekce použít. Funkce f ( x ) = e x − 2 cos x je spojitá funkce, navíc f (0 ) = −1.000000 a f (1) = 1.637677 (funkce nabývá v krajních bodech intervalu hodnoty s opačným znaménkem), tedy metodu bisekce můžeme použít. Výsledky jednotlivých výpočtů podle algoritmu jsou seřazeny do následující tabulky.
2.Řešení nelineárních rovnic krok k 0 1 2 3 4 5 6
21
f ( sk )
ak
bk
sk
0.000000 0.500000 0.500000 0.500000 0.500000 0.531250 0.531250
1.000000 1.000000 0.750000 0.625000 0.562500 0.562500 0.546875
0.500000 0.750000 0.625000 0.562500 0.531250 0.546875 0.539063
<0 >0 >0 >0 <0 >0
bk − ak 1.000000 0.500000 0.250000 0.125000 0.062500 0.031250 0.015625
V šestém kroku algoritmus skončil, protože chyba 6-té iterace je již v toleranci požadované přesnosti, 0.015625 ≤ ε . Za aproximaci kořene zjištěné metodou bisekce s požadovanou přesností ε = 0.02 lze považovat kterékoliv číslo z intervalu I 6 = [0.531250; 0.546875] . Metoda regula-falsi. Velice podobnou metodou jako metoda bisekce je metoda regula-falsi. Mějme tedy reálnou funkci f spojitou na intervalu I 0 = [a0 ,b0 ] a předpokládejme, že funkční hodnoty v koncových bodech intervalu mají opačná znaménka, tedy f (a0 ) ⋅ f (b0 ) < 0 . Přirozená otázka je zda nelze interval místo „půlení“ dělit jiným bodem sk , který by využil vlastnosti samotné funkce f k urychlení konvergence. Takovým bodem může být například průsečík přímky procházející body A = [ a0 , f (a0 ) ] a B = [ b0 , f (b0 ) ] s osou x (viz obrázek).
A
s0 a0
b0
y = f (x )
B
Snadno se dá ukázat, že sk = ak −
f (ak ) (bk − ak ) . f (bk ) − f (ak )
Samotný algoritmus metody regula-falsi je potom analogický algoritmu metody bisekce, jako zastavovací podmínka se zpravidla používá f (sk ) ≤ δ :
22
Vstup: δ , a := a0 , b := b0 Opakuj
a := s
f (a ) (b − a ) f (b ) − f (a ) jestliže f (a ) ⋅ f (s ) > 0
b := s
jestliže
konec
jestliže
s := a −
dokud f (s ) ≤ δ .
f (a ) ⋅ f (s ) < 0
f (a ) ⋅ f (s ) = 0
Výstup: I k = [a, b] .
Metoda regula-falsi je vždy konvergentní metoda, není však vhodná pro užití blízko kořene. Říkáme, že posloupnost {xk } konverguje k číslu α rychlostí řádu r jestliže pro k → ∞ platí r
(
xk +1 − α = C xk − α + o xk − α
(Připomeňme,
že
symbol
g (k ) = o(h(k ))
r
) pro C ≠ 0 .
pro
k →∞
označuje,
že
⎡ g (k ) ⎤ lim ⎢ ⎥ = 0 .) k →∞ ⎣ h(k ) ⎦ Prakticky to znamená: čím větší r , tím rychlejší konvergence. Dá se ukázat, že metoda regula-falsi konverguje rychlostí řádu r = 1 . Mějme však pořád na paměti, že rychlost konvergence vystihuje chování asymptotické pro k → ∞ , což se na několika prvních iteracích se nemusí znatelně projevit. Příklad. Určete metodou regula-falsi nejmenší kladný kořen rovnice e x − 2 cos x = 0 , proces zastavte bude-li f (sk ) ≤ δ = 0.02 . Jako počáteční aproximaci použijte
interval zjištěný grafickou metodou I 0 = [0;1] , pracujte na 6 desetinných míst.
Řešení: Metodu lze použít, neboť podmínky jsou stejné jako u metody bisekce, ty jsme již ověřili v předchozím příkladě. Výsledky výpočtů jsou uvedeny v tabulce. krok k 0 1 2 3
f ( sk )
ak
bk
sk
0.000000 0.379121 0.500260 0.530577
1.000000 1.000000 1.000000 1.000000
0.379121 0.500260 0.530577 0.537668
Hledaný kořen uvedené rovnice s přesností
-0.396981 -0.105766 -0.025118 -0.005801
f (x ) ≤ 0.02 je x = 0.537668 .
2.Řešení nelineárních rovnic Metoda prosté iterace.
Věta. Nechť funkce ϕ je na nějakém intervalu I spojitá a platí 1. ∀x ∈ I : ϕ ( x ) ∈ I (tzv. zobrazení „do sebe“) 2. ∃q ∈ [0;1) : ϕ ( x ) − ϕ ( y ) ≤ q x − y ∀x, y ∈ I (kontraktivní zobrazení). Potom v intervalu I posloupnost
{ }
ke kořenu α .
xk ∞ k =0
existuje právě jeden kořen α rovnice x = ϕ ( x ) a
určena formulí x´k +1 = ϕ ( xk ) konverguje pro každé x0 ∈ I
Metoda prosté iterace vychází z tvrzení této věty. Přepišme rovnici f ( x ) = 0 do tvaru x = ϕ ( x ) tak, aby funkce ϕ splňovala předpoklady věty a sestrojme posloupnost iterací x´k +1 = ϕ ( xk ) pro nějaké x0 ∈ I . Jednotlivé členy posloupnosti pak můžeme považovat za aproximace hledaného řešení α .
Poznámka. Podmínku 2. předchozí věty lze pro diferencovatelnou funkci nahradit podmínkou ∃q ∈ [0;1) : ϕ ′( x ) ≤ q ∀x ∈ I .
Poznamenejme ještě, že podmínky věty jsou ke konvergenci metody prosté iterace pro libovolnou počáteční aproximaci x0 ∈ I postačující, nikoliv však nutné. Můžou tedy nastat případy, že funkce ϕ nesplňuje tyto předpoklady a proces přesto konverguje pro nějakou počáteční aproximaci. Nemáme-li jinou možnost jak rovnici f ( x ) = 0 přepsat, můžeme provést několik iterací a z tabulky ihned vidíme, zda metoda konverguje, nebo nekonverguje. Je-li funkce ϕ diferencovatelná, dá se využitím Taylorova rozvoje ukázat, že rychlost konvergence metody prosté iterace je r = 1 pokud ϕ ′(α ) ≠ 0 , případně r = 2 pokud ϕ ′(α ) = 0, ϕ ′′(α ) ≠ 0 . Algoritmus: Vstup: δ , x0 ∈ I Opakuj xk +1 := ϕ ( xk ) dokud f ( xk +1 ) ≤ δ
Výstup: xk +1 Příklad. Aproximujte metodou prosté iterace e x − 2 cos x = 0 , proces zastavte bude-li
nejmenší kladný kořen rovnice f ( xk +1 ) ≤ δ = 0.005 . Jako počáteční
aproximaci použijte interval zjištěný grafickou metodou I 0 = [0;1] , pracujte na 6 desetinných míst.
23
24
Řešení: V našem případě můžeme použít dva způsoby přepisu rovnice e x − 2 cos x = 0 , a to do tvaru ⎛ ex ⎞ a) x = arccos⎜ ⎟ a ⎜ 2⎟ ⎝ ⎠ b) x = ln (2 ⋅ cos x ) . V případě a) dostaneme iterační předpis ⎛ e xk xk +1 = arccos⎜ ⎜ 2 ⎝
⎞ ⎟, ⎟ ⎠
v případě b) předpis
xk +1 = ln (2 ⋅ cos xk ) . Jednotlivé iterace jsou seřazeny v tabulce. V obou případech jsme vzali za počátečnou aproximaci x0 = 0.5 .
iterace
1 2 3 4 5 6
⎛ e xk xk +1 = arccos⎜ ⎜ 2 ⎝ x0 = 0.5 0.601725 0.421119 0.704658 0.152020
⎞ ⎟ ⎟ ⎠
xk +1 = ln (2 ⋅ cos xk ) x0 = 0.5 0.562563 0.525782 0.548042 0.534792 0.542760 0.537997
Z tabulky vidíme, že v případě iteračního předpisu ⎛ e xk ⎞ ⎟ xk +1 = arccos⎜ ⎜ 2 ⎟ ⎝ ⎠ proces diverguje, zároveň lze snadno ověřit, že funkce ⎛ ex ⎞ ϕ ( x ) = arccos⎜⎜ ⎟⎟ ⎝ 2⎠ nesplňuje na intervalu I = [0,1] předpoklady uvedené věty (nepomůže ani zmenšení intervalu např. na I = [0,0.7 ] ). V případě iteračního předpisu xk +1 = ln (2 ⋅ cos xk ) jsou podmínky věty splněny, metoda konverguje. Hledaný kořen uvedené rovnice s danou přesností je x6 = 0.537997 .
2.Řešení nelineárních rovnic
2.3.2 Metody zpřesňovací. Newtonova metoda (Metoda tečen). Hledáme kořen rovnice f ( x ) = 0 na intervalu I , předpokládejme, že funkce f je na I spojitá a diferencovatelná, tj. pro každé x ∈ I existuje f ′( x ) . Můžeme tedy použít Taylorova rozvoje funkce f v okolí bodu x0 ∈ I : f ( x ) = f ( x0 ) + f ′( x0 )( x − x0 ) + ...
Pro kořen α platí f (α ) = 0 , tedy 0 = f (α ) = f ( x0 ) + f ′( x0 )(α − x0 ) + ...
Vezmeme-li v úvahu pouze první dva členy rozvoje, můžeme za aproximaci kořene považovat x1 f ( x0 ) x1 = x0 − . f ′( x0 ) Nyní použijme Taylorův rozvoj funkce f v okolí bodu x1 , analogicky obdržíme přesnější aproximaci kořene x2 : f ( x1 ) x2 = x1 − . f ′( x1 )
Opakováním tohoto postupu dostaneme iterační posloupnost {xk }∞ k = 0 , kde x0 je zvolená počáteční aproximace a
f ( xk ) . f ′( xk ) Na místě je však otázka, zda bude taková posloupnost vždy konvergovat, a pokud ano tak zda bude konvergovat k hledanému kořenu α . Postačující podmínku konvergence nám dává následující věta. xk +1 = xk −
Věta. Nechť funkce f je na intervalu I = [a, b] spojitá a alespoň dvakrát diferencovatelná. Nechť f ′( x ) ≠ 0 pro každé x ∈ [a, b] a nechť pro každé x ∈ [a, b] f ′′( x ) nemění znaménko. Potom pokud f (a ) ⋅ f (b ) < 0 a zároveň f (a )
f (b ) < b − a , tak Newtonova metoda konverguje pro každé f ′(b )
x0 ∈ I . Samotný algoritmus Newtonovy metody je analogický algoritmu metody prosté iterace: Vstup: δ , x0 ∈ I 0 Opakuj xk +1 = xk − dokud f ( xk +1 ) ≤ δ Výstup: xk +1
f ( xk ) f ′( xk )
25
26
Poznámka. Uvědomíme-li si, že první derivace funkce je směrnicí tečny v daném bodě, potom každá iterace xk +1 je průsečíkem osy x a tečny ke křivce y = f ( x ) , vedené bodem [xk , f ( xk )] (viz graf), proto se metoda často nazývá metodou tečen. Newtonova metoda konverguje rychlostí řádu r = 2 . f ( xk ) Jako zastavovací podmínku lze také použít <ε. f ′( xk )
x0
x1
x2
y = f (x )
Příklad. Aproximujte Newtonovou metodou nejmenší kladný kořen rovnice x e − 2 cos x = 0 , jako počáteční aproximaci použijte hodnotu x0 = 0.5 , pracujte na 8 desetinných míst.
Řešení: Derivace funkce f ( x ) = e x − 2 cos x je funkce jednotlivých iterací jsou seřazeny v tabulce. iterace k
xk
1 2 3 4
0.54082105 0.53978583 0.53978516 0.53978516
f ′( x ) = e x + 2 sin x . Hodnoty
Vidíme, že na uvedený počet desetinných míst se již aproximace nemění od čtvrtého kroku. Přesnost iterace x3 = 0.53978516 je f (0.53978516) = 2.2 ⋅ 10−9 .
2.Řešení nelineárních rovnic
27
Metoda sečen. Na metodu sečen lze pohlížet jako na „zobecnění“ Newtonovy metody pro funkce, které nejsou diferencovatelné. Víme, že derivaci možno aproximovat tzv. diferenčním podílem f ( xk ) − f ( xk −1 ) f ′( xk ) ≈ , xk − xk −1 což tedy znamená, že nevedeme bodem [xk , f ( xk )] tečnu, ale vedeme dvěma body [xk −1, f ( xk −1 )] , [xk , f ( xk )] sečnu- proto název metoda sečen. V praxi to ovšem znamená, že na začátku potřebujeme znát dvě počáteční aproximace x0 , x1 a iterační předpis bude mít tvar xk − xk −1 xk +1 = xk − f ( xk ) . f ( xk ) − f ( xk −1 ) Metoda sečen patří mezi metody interpolační, neboť tento postup můžeme taky chápat jako nahrazení funkce f lineárním interpolačním polynomem určeným hodnotami [xk −1, f ( xk −1 )] , [xk , f ( xk )] a hledáním kořene tohoto polynomu, místo 1 samotné funkce f . Rychlost konvergence je r = 1 + 5 ≈ 1.62 , tedy o něco 2 menší než u metody Newtonovy.
(
)
Algoritmus. Vstup: δ , x0 , x1 ∈ I 0 Opakuj xk +1 = xk − f ( xk ) dokud f ( xk +1 ) ≤ δ
xk − xk −1 f ( xk ) − f ( xk −1 )
Výstup: xk +1
Příklad. Aproximujte metodou sečen nejmenší kladný kořen rovnice e x − 2 cos x = 0 , jako počáteční aproximace použijte hodnoty x0 = 0, x1 = 1 , pracujte na 6 desetinných míst.
Řešení: Hodnoty jednotlivých iterací jsou seřazeny v tabulce. Iterace k
xk
0 1 2 3 4 5
0.000000 0.379121 0.500260 0.544256 0.539672 0.539785
28
6
0.539785
Z tabulky vidíme, že na daný počet desetinných míst se iterace od 5-tého kroku nemění. Přesnost iterace x5 = 0.539785 je f (0.539785) ≈ 0.8 ⋅ 10−6 . Metoda sečen nemusí vždy konvergovat, je nutné, aby počáteční aproximace byly již „dostatečně přesné“. Je proto vhodné kombinovat ji s metodou startovací, a to konkrétně metodou regula-falsi, neboť je principialně stejná. Jak jsme se již zmínili, metody startovací i zpřesňovací mají svoje výhody i nevýhody. Pro dosažení co nejlepšího výsledku při řešení konkrétního příkladu je vhodné oba typy kombinovat. Seznámili jsme se s metodami startovacími i zpřesňovacími, a jsme tedy připraveni komplexně vyřešit úlohu.
Příklad. Hospodář má u domu kruhový trávník o poloměru 1 a jde dát kozu na pastvu. Chce jí připevnit provázkem na obvod trávníku tak, aby mu mohla vyžrat maximálně polovinu trávy. Jak dlouhý musí být provázek?
Řešení: Převeďme nejdříve slovní úlohu na úlohu matematickou, případně numerickou. Celý problém můžeme načrtnout takto:
P k1
k2
ρ =1 S
R x
A
r
Q
Kruh omezený kružnicí k1 , se středem S a poloměrem ρ = 1 , představuje trávník, který chceme rozdělit kružnicí k2 , se středem A, A ∈ k1 a poloměrem r , na dvě obsahově stejné části.
2.Řešení nelineárních rovnic Označili jsme P, Q průsečíky obou kružnic a velikost úhlu ∠PAQ = x . Našim cílem je stanovit velikost poloměru r . Soustřeďme se na dva shodné, rovnoramenné trojúhelníky ΔQSA a ΔPSA . Rovnoramenné proto, že strany QS = AS = r = 1 a podobně PS = AS = r = 1 . Strana SA rozděluje úhel x na polovinu, tedy x ∠PAS = = ∠QAS . 2 Víme, že v rovnoramenném trojúhelníku leží proti stejným stranám stejné úhly, tudíž x x ∠APS = ∠PAS = a ∠AQS = ∠QAS = . 2 2 Součet úhlů v trojúhelníku je 180° neboli π , a tak zbylé úhly v trojúhelnících budou ∠PSA = π − x = ∠QSA . Navíc úsečky PQ a AS jsou na sebe kolmé. Označme dále S k1 obsah kruhu vymezeného kružnicí k1 (obsah trávníku), a S1, S 2 obsahy jednotlivých vrchlíků: S1 vrchlík kružnice k1 vymezený tětivou PQ S 2 vrchlík kružnice k2 vymezený tětivou PQ , 1 přičemž požadujeme, aby S1 + S 2 = S k1 (kružnice k2 dělí kruh k1 na dvě části 2 se stejným obsahem).
Obsah kruhu S k1 = π ⋅ 12 = π .
Obsah vrchlíku S1 (poloměr je ρ = 1 a úhel výseče 2(π − x ) ) bude 1 S1 = ⋅ 12 ⋅ (2(π − x ) − sin (2(π − x ))) 2 a podobně obsah vrchlíku S 2 (poloměr je r a úhel výseče x ) bude 1 S 2 = r 2 ( x − sin x ) . 2 1 S k dostaneme 2 1 1 1 1 ⋅ (2(π − x ) − sin (2(π − x ))) + r 2 ( x − sin x ) = π , 2 2 2 co představuje jednou rovnici o dvou neznámých r a x . Potřebujeme ještě druhou rovnici, vztah mezi r a x . Po dosazení do vztahu S1 + S 2 =
Užijeme-li vlastností pravoúhlých trojúhelníku ΔRSQ a ΔARQ obdržíme RQ sin (π − x ) = 1
29
30
sin tedy r=
x RQ , = 2 r
sin (π − x ) sin x x = = 2 cos . x x 2 sin sin 2 2
Nyní máme dvě rovnice o dvou neznámých 1 1 1 ⋅ (2(π − x ) − sin (2(π − x ))) + r 2 ( x − sin x ) = π 2 2 2 x r = 2 cos 2 které můžeme řešit dvěma způsoby: 1. z druhé rovnice vyjádříme r a dosadíme do rovnice první (tím dostaneme rovnici vzhledem k neznáme x , vyřešíme jí a poté dopočteme neznámý poloměr r 2. nebo z druhé rovnice vyjádříme x , dosadíme do první rovnice a jejím vyřešením dostaneme r . Hodnota x nás nezajímá, tudíž jí zpětně dopočítat nemusíme. Možnost 2. je matematicky i numericky korektnější, neboť počítá hledanou hodnotu přímo a nezatěžuje výpočet chybou průběžného výpočtu „pomocné“ hodnoty x . Tuto cestu ponecháme čtenáři jako domácí cvičení a ukážeme si možnost 1., která je trošku přehlednější. Provedeme-li popsané úpravy a užijeme součtových vzorců, vede rovnice na tvar (π − x ) + (1 + cos x )x − sin x = 1 π , 2 nebo po úpravě
π
+ x cos x − sin x = 0 . 2 Tato rovnice není analyticky řešitelná, jediná možnost je řešení numerické. Řešení budeme hledat ve třech krocích: graficky určíme interval, kde se hledaný kořen rovnice nachází použijeme startovací metodu na získání dostatečně kvalitní aproximace řešení použijeme metodu zpřesňovací. Načrtněme nejdříve graf. Řešíme rovnici
π 2
+ x cos x − sin x = 0 .
2.Řešení nelineárních rovnic
31
1/2*PI - sin(x) + x*cos(x) y
3.75
2.5
1.25
0 0
1.25
2.5
3.75
5 x
-1.25
Záporné hodnoty x nás nezajímají, navíc hodnota x ∈ (0; π ) neboť se jedná o úhel (viz náčrtek úlohy). Jediné řešení rovnice, které má pro naší úlohu smysl leží v intervalu I 0 = [1.25; 2.50] . Použijeme-li jako startovací metodu bisekci pro takto stanovený počáteční interval, obdržíme hodnoty Iterace k
ak
bk
sk
0 1 2 3
1.25 1.875 1.875 1.875
2.50 2.50 2.1875 2.03125
1.875 2.1875 2.03125 1.953125
f ( sk ) >0 <0 <0 <0
Funkce je diferencovatelná, můžeme tedy jako zpřesňovací metodu použít Newtonovou metodu. Derivace funkce f ( x ) =
π
+ x cos x − sin x je f ′( x ) = − x sin x . 2 Hodnoty získané Newtonovou metodou s počátečnou aproximací x0 = 1.953125 jsou uvedeny v tabulce. Iterace k
xk
0 1 2 3 4
1.953125 1.905844186 1.905695731 1.905695729 1.905695729
Přesnost x3 = 1.905695729 je 10−9 .
32
x Nyní dopočteme hledanou hodnotu poloměru r z druhé rovnice r = 2 cos , 2 1 . 905695729 ⎛ ⎞ r = 2 cos⎜ ⎟ ≈ 1.158728473 . 2 ⎝ ⎠ x Označme F ( x ) = 2 cos a A( x ) absolutní chybu, se kterou jsme získali hodnotu 2 x. Potom absolutní chyba poloměru r bude dF A( x ) , A(r ) = dx x tedy A(r ) = − sin A( x ) . Po dosazení vidíme, že řád chyby zůstává stejný, tedy 2 −9 10 . Odpověď. Hospodář musí uvázat kozu na provázek dlouhý 1.158729473 jednotky.
2.4 Kontrolní otázky, korespondenční úkol. Kontrolní otázky.
1. Určete numericky největší záporný kořen rovnice sin 2 x + x = 0 s přesností δ = 10−8 . Použijte kombinaci metody bisekce a Newtonovy metody. Pracujte na 10 desetinných míst. 2. Určete numericky druhý nejmenší kladný kořen rovnice 5ln x sin x = 0 a) metodou bisekce b) metodou regula-falsi. Výsledek zpřesněte na δ = 10−8 c) metodou Newtonovou d) metodou sečen. 3. Nalezněte první kladný kořen rovnice ln x + x = 0 metodou bisekce a prosté iterace. Výsledek zpřesněte Newtonovou metodou.
2.Řešení nelineárních rovnic Část pro zájemce.
1. Pokuste se rozšířit metodu prosté iterace pro řešení soustavy nelineárních rovnic. Řešte metodou prosté iterace soustavu nelineárních rovnic x 2 − x sin y = 3 5 x − y 2 = 2 ln xy (Pomůcka: použijte vyjádření
x = 3 + x sin y y = 5 x − 2 ln xy a jako počátečnou aproximaci vezměte vektor (1,1)T .) 2. Pokuste se rozšířit Newtonovou metodu pro řešení soustavy nelineárních r r r rovnic. (Pomůcka: Zapíšeme-li soustavu ve tvaru F ( x ) = 0 , derivace v
⎛ ∂F ⎞ iteračním předpisu přejde v tzv. Jacobiovu matici ⎜ i ⎟ , kde Fi jsou složky ⎜ ∂x j ⎟ ⎠ ⎝ r r F a 0 je nulový vektor.) Načrtněte řešení soustavy x 2 − x sin y = 3 . 5 x − y 2 = 2 ln xy
Shrnutí kapitoly. Cílem této kapitoly bylo seznámení se základními numerickými metodami pro řešení nelineárních rovnic a jejich vhodným použitím z hlediska výhod a nevýhod, jakož i pochopení principu iteračních metod, které se v numerice velice často objevují. Pokud jste zvládli zpracovat algoritmy všech popsaných metod a vyřešili kontrolní otázky, pokuste se samostatně vypracovat korespondenční úkol.
Korespondenční úkol. Nepřímým měřením se zjišťovala hodnota a pomocí hodnot x, y přičemž funkční
závislost a = a( x, y ) byla a = 2
x2 + e 3 y , kde hodnota x byla získána jako y
kladný kořen rovnice e x − 2( x − 1) = 0. Hodnota y = 0.73721 ± 0.00005 . Diskutujte použití alespoň tří postupů různých metod vzhledem k faktu, že hodnotu a potřebujeme získat s přesností ε . Poté spočtěte hledanou hodnotu a s přesností ε = 10−2 každým z uvedených postupů. Diskutujte možnosti zpřesnění tohoto měření. 2
33
34
2. Řešení soustav lineárních rovnic
3. ŘEŠENÍ SOUSTAV LINEÁRNÍCH ROVNIC V této kapitole se dozvíte:
•
jak numericky řešit soustavy lineárních rovnic;
•
co jsou metody příme, jakou roli hraje výběr hlavního prvku a kdy je nutno ho použít;
•
jak využít tyto metody k numerickému výpočtu inverze matic;
•
jak pracují metody iterační;
•
jak využít speciálního tvaru matice soustavy k numerickému řešení.
Budete schopni:
•
numericky vyřešit soustavu lineárních rovnic metodami přímými (Gaussovou eliminační metodou a metodou LU rozkladu) a rozhodnout, zda je nutno použít metodu s výběrem hlavního prvku, či nikoliv;
•
řešit soustavu lineárních rovnic metodami iteračními;
•
použít tyto metody pro nalezení inverze matice.
Klíčová slova této kapitoly: matice soustavy, rozšířená matice soustavy,hlavní prvek eliminace, hlavní rovnice eliminace, metoda přímá, metoda iterační
Čas potřebný k prostudování učiva kapitoly: 10 + 10 hodin (teorie + řešení úloh)
Průvodce studiem. V této kapitole popíšeme dvě třídy numerických metod řešení soustav lineárních rovnic, a to metody přímé a iterační. Na konkrétních příkladech pak ukážeme možnosti použití jednotlivých metod z hlediska jejich výhod a nevýhod. Tato kapitola předpokládá základní znalosti z lineární algebry, případně teorie matic.
35
36 3.1 Formulace problému.
Nechť
a11x1 + a12 x2 + ... + a1n xn = b1 a21x1 + a22 x2 + ... + a2 n xn = b2 ... am1x1 + am 2 x2 + ... + amn xn = bm
je soustava m lineárních rovnic o n neznámých x1,..., xn . Řešit soustavu lineárních rovnic znamená určit n-tici čísel x1,..., xn vyhovujících všem m rovnicím. r r Označíme-li A = aij , b = (b1, ..., bm )T a x = ( x1, ..., xn )T , můžeme soustavu
( )
zapsat v maticovém tvaru
r r Ax = b . r Matice A se nazývá matice soustavy, sloupcový vektor b nazýváme vektor r pravých stran a matici A b rozšířená matice soustavy.
( )
Připomeňme, že soustava rovnic má řešení právě tehdy, když hodnost r matice soustavy je rovna hodnosti rozšířené matice soustavy h( A) = h A b .
( )
Je-li m = n a matice A regulární (tj. determinant je nenulový), potom pro daný r r r vektor b existuje právě jedno řešení soustavy Ax = b , tzv. teoretické (přesné) řešení r r xt = A−1b . r Je-li h( A) = h A b < n má soustava rovnic nekonečně mnoho řešení, tzv. obecné
( )
řešení, ve tvaru
r r n − h( A) r x = x0 + ∑α i ui , i =1
r r pro libovolné čísla α i , kde x0 je řešení nehomogenní soustavy a ui jsou lineárně nezávislé vektory, tzv. řešení příslušné homogenní soustavy.
Numerické metody řešení soustav lineárních rovnic dělíme do dvou skupin, na metody 1. přímé (finitní) 2. iterační (přibližné). Metody přímé jsou takové, které pokud pomineme chyby zaokrouhlovací, vedou k teoretickému, přesnému řešení. Metody iterační konstruují iterační posloupnost přibližných řešení, po které chceme, aby konvergovala k řešení přesnému. K chybě zaokrouhlovací přistupuje navíc chyba metody. U všech metod nás budou zajímat především tři věci. Rychlost metody- tedy počet potřebných aritmetických operací, požadavky na paměť počítače a přesnost vypočteného řešení.
2. Řešení soustav lineárních rovnic
37
3.2 Metody přímé. Gaussova eliminační metoda. Jedna z nejstarších a nejznámějších metod řešení soustav lineárních rovnic je Gaussova eliminační metoda. Je založená na vlastnosti, že ekvivalentní úpravy matice (tedy násobení řádku nenulovým číslem, přičtení některého řádku k jinému) nemají na řešení vliv- řešení soustavy zůstává stejné.
Metoda postupuje tak, že rozšířenou matici soustavy upravuje ekvivalentními maticovými úpravami až do tvaru trojúhelníkového (nejčastěji horního) a poté substitucí (zpětnou) dopočte řešení. Ukažme si tento postup na konkrétním příkladě. Příklad. Řešte soustavu
x1 − x2 + x3 = 2 2 x1 + x2 − 3 x3 = −5
− x1 + 3x2 − x3 = 2 Řešení: Rozšířená matice soustavy bude ⎛ 1 −1 1 2 ⎞ r ⎜ ⎟ A b = ⎜ 2 1 −3 − 5⎟ , ⎜−1 3 −1 2 ⎟⎠ ⎝ tuto matici budeme upravovat na horní trojúhelníkový tvar tak, že pomocí prvku na pozici hlavní diagonály a11 = 1 vhodným vynásobením prvního řádku a přičtením ke zbylým řádkům vyloučíme prvky na pozicích a21 = 2 a a31 = −1 , čili prvky prvního sloupce ležící pod hlavní diagonálou. První řádek vynásobíme číslem − 2 (tzv. multiplikátorem) a přičteme k řádku druhému
( )
⎛ 1 −1 1 2 ⎞ ⎛ 1 −1 1 ⎜ ⎟ ⎜ − 5⎟ ≈ ⎜ 0 3 −5 ⎜ 2 1 −3 ⎜−1 3 −1 ⎟ ⎜ 2 ⎠ ⎝−1 3 −1 ⎝ a analogicky první řádek vynásobíme multiplikátorem 1 třetímu
2 ⎞ ⎟ − 9⎟ 2 ⎟⎠ a přičteme k řádku
⎛ 1 −1 1 2 ⎞ ⎛1 −1 1 2 ⎞ ⎜ ⎟ ⎜ ⎟ 3 −5 − 9⎟ ≈ ⎜ 0 3 − 5 − 9⎟ . ⎜0 ⎜−1 3 −1 2 ⎟⎠ ⎜⎝ 0 2 0 4 ⎟⎠ ⎝ Nyní se nám povedlo vyloučit- eliminovat- všechny prvky prvního sloupce pod hlavní diagonálou. Provedli jsme tzv. 1 fázi eliminace.
38
Stejným způsobem pokračujeme dále eliminací prvků druhého sloupce ležících pod hlavní diagonálou. Opět pomocí prvku na diagonále a22 = 3 (již naší upravené matice) budeme vylučovat prvky pod diagonálou- v našem případě prvek a32 = 2 . Musíme tedy druhý řádek pronásobit multiplikátorem − 2 a přičíst ke třetímu 3 ⎛1 −1 1 ⎜ ⎜0 3 − 5 ⎜0 2 0 ⎝
2 ⎞ ⎛⎜ 1 − 1 1 ⎟ − 9⎟ ≈ ⎜0 3 − 5 ⎜ 4 ⎟⎠ ⎜ 0 0 10 3 ⎝
2 ⎞⎟ − 9⎟ . ⎟ 10 ⎟ ⎠
Realizací tzv. 2. fáze eliminace jsme matici převedli na horní trojúhelníkový tvar. Samotné řešení obdržíme zpětnou substitucí. Z posledního řádku (představuje třetí rovnici) dostáváme 10 ⋅ x3 = 10 , tedy x3 = 3 , 3 dosadíme do druhého řádku (rovnice) 3 x2 − 5 x3 = −9 , po vyčíslení x2 = 2 , a konečně dosazením do prvního řádku (rovnice) x1 − x2 + x3 = 2 dostaneme x1 = 1 . r Řešení uvedené soustavy tří rovnic o třech neznámých je vektor x = (1, 2, 3)T .
Poznámka. Multiplikátory příslušných fází eliminace jsou v podstatě podíly aij( j −1) mij = − ( j −1) , a jj
kde j značí index fáze eliminace, index i příslušný řádek, ke kterému multiplikátor náleží a horní index ( j − 1) u čísel aij , a jj znamená, že se jedná o čísla v příslušných pozicích matice vzniklé po úpravami se prvky matice mění).
( j − 1)
-fázi eliminace (neboť
V našem případě jsme měli multiplikátory první fáze m21 = −2 , m31 = 1 a 2 multiplikátor druhé fáze m32 = − . 3 Označme L dolní trojúhelníkovou matici s jedničkami na diagonále a s prvky − mij pod diagonálou, v našem případě ⎛1 0 0 ⎞⎟ ⎜ L=⎜ 2 1 0⎟ . ⎟ ⎜ ⎜−1 2 1⎟ 3 ⎠ ⎝
2. Řešení soustav lineárních rovnic
39
Označme dále U horní trojúhelníkovou matici, vzniklou z matice A eliminací, tedy ⎛1 −1 1 ⎞ ⎟ ⎜ U = ⎜0 3 − 5 ⎟ . ⎟ ⎜ ⎜ 0 0 10 ⎟ 3⎠ ⎝ Snadno lze ověřit, že platí LU = A . Takový rozklad (na dolní a horní trojúhelníkovou matici), říkáme mu LU -rozklad matice, lze provést pro každou čtvercovou matici A , pro kterou determinanty det Ak ≠ 0 pro k = 1,..., n − 1 , kde Ak je matice typu (k, k ) sestavená ze společných prvků prvních k-řádků a ksloupců matice A . Budeme-li požadovat, aby matice L měla na diagonále jedničky, je rozklad jednoznačný. Algoritmus Gaussovy eliminační metody můžeme schématicky zapsat takto:
( ) ( ( ) ( )
)
r T Vstup: n , A = aij(0 ) , b = a1(,0n)+1, ..., an(0, n) +1 (tedy rozšířená matice soustavy r A b = aij(0 ) pro i = 1,..., n a j = 1,..., n, n + 1 ).
Opakuj pro
(
k = 1,2,..., n − 1 Opakuj pro i = k + 1, k + 2,..., n a (k −1) mik = − ik(k −1) akk j = k + 1, k + 2,..., n, n + 1 Opakuj pro a (k ) = a (k −1) + m a (k −1) .
)
(
ij
ij
)
Výstup: U = aij(i −1) , y = a1(,0n)+1, a2(1,)n +1, a3(2, n)+1,..., an(n, n−+11) .
ik kj
Gauss-Jordanova metoda. Vraťme se zpátky k samotnému procesu eliminace. U Gaussovy eliminační metody jsme skončili eliminaci ve chvíli, kdy byla matice převedena na horní trojúhelníkový tvar, tj. byly „vynulovány“ prvky pod hlavní diagonálou a samostatně jsme provedli zpětnou substituci. Nic nám však nebrání zahrnout proces substituce, neboli dopočtení samotného řešení do procesu eliminace. Budeme-li pokračovat z horní trojúhelníkové matice eliminací tak, abychom dospěli až k matici jednotkové, potom řešení soustavy budou prvky na pozicích původních pravých stran v rozšířené matici soustavy. Ukažme si to na předchozím příkladě.
Gaussovou eliminační metodou jsme převedli rozšířenou matici soustavy na tvar
40 ⎛1 −1 1 2 ⎞⎟ ⎜ ⎜0 3 − 5 − 9⎟ . ⎜ ⎟ ⎜ 0 0 10 10 ⎟ 3 ⎝ ⎠ Pokračujme nyní v eliminaci dále. Nejdříve „vynulujeme“ prvky druhého sloupce ležící nad hlavní diagonálou. K prvnímu řádku přičteme 1 -násobku řádku 3 druhého ⎛1 −1 1 2 ⎞⎟ ⎛⎜ 1 0 − 2 − 1 ⎞⎟ ⎜ 3 ⎜0 3 − 5 − 9⎟ ≈ ⎜ 0 3 − 5 − 9⎟ . ⎜ ⎟ ⎜ ⎟ ⎜ 0 0 10 10 ⎟ ⎜ 0 0 10 10 ⎟ 3 3 ⎝ ⎠ ⎝ ⎠ Nyní „vynulujeme“ prvky třetího sloupce ležící nad diagonálou, k prvnímu řádku přičteme 1 -násobku třetího řádku a ke druhému přičteme 3 -násobku třetího 5 2 řádku ⎛1 0 − 2 − 1 ⎞⎟ ⎛⎜ 1 0 0 1 ⎞⎟ ⎜ 3 ⎜0 3 − 5 − 9⎟ ≈ ⎜ 0 3 0 6 ⎟. ⎜ ⎟ ⎜ ⎟ ⎜ 0 0 10 10 ⎟ ⎜ 0 0 10 10 ⎟ 3 3 ⎝ ⎠ ⎝ ⎠ Matice je ve tvaru diagonálním, můžeme jí tedy upravit na jednotkovou vhodným vynásobením každého z řádků: první řádek není nutno upravovat, druhý řádek vynásobíme 1 a třetí řádek číslem 3 : 3 10 ⎛1 0 0 1 ⎞⎟ ⎛ 1 0 0 1⎞ ⎜ ⎜ ⎟ ⎜0 3 0 6 ⎟ ≈ ⎜0 1 0 2⎟ . ⎜ ⎟ ⎜ 0 0 10 10 ⎟ ⎜⎝ 0 0 1 3 ⎟⎠ 3 ⎝ ⎠ Matice soustavy je ve tvaru jednotkovém, řešení soustavy se objevilo na pozici vektoru pravých stran, neboť jednotlivé řádky rozšířené matice představují rovnice =1 1 ⋅ x1
1 ⋅ x2
=2 1 ⋅ x3 = 3
Tato modifikace Gaussovy eliminační metody se nazývá Gauss-Jordanova metoda. Přednosti jsou především v tom, že jí lze snadno použít k řešení maticových rovnic AX = B , kde A je matice typu (m, n ) , B je typu (m, r ) a řešení X je matice typu (n, r ) . Taková maticová rovnice představuje řešení r soustav m lineárních rovnic o n neznámých. Tuto metodu lze také použít pro stanovení řešení nedourčené soustavy lineárních rovnic (míň rovnic než neznámých).
2. Řešení soustav lineárních rovnic Příklad. Řešte soustavu rovnic:
x1 + 3 x2 − x3 = 1 x1 − x2 + x3 = 2
41
.
Řešení: r Hodnost matice soustavy a rozšířené matice soustavy jsou h( A) = h A b = 2 , r r r tedy řešení bude ve tvaru x = x0 + α u , pro libovolné α ∈ R . Rozšířenou matici soustavy ⎛1 3 −1 1⎞ ⎜ ⎟ ⎜1 −1 1 2⎟ ⎜0 0 0 0⎟ ⎝ ⎠
( )
musíme upravit, abychom mohli naleznout řešení a to tak, že nulu na diagonále v posledním řádku nahradíme jedničkou a přidáme k vektoru pravých stran další sloupcový vektor (0,0,1)T , tedy ⎛1 3 −1 1 0⎞ ⎜ ⎟ ⎜1 −1 1 2 0⎟ . ⎜0 0 1 0 1⎟ ⎝ ⎠ Tuto soustavu řešíme Gauss-Jordanovou metodou. ⎛ 1 3 − 1 1 0 ⎞ ⎛ 1 3 − 1 1 0 ⎞ ⎛⎜ 1 0 − 1 7 0 ⎞⎟ 4 ⎜ ⎟ ⎜ ⎟ ⎜ 1 − 1 1 2 0⎟ ≈ ⎜ 0 − 4 0 1 0⎟ ≈ ⎜ 0 − 4 0 1 0⎟ ≈ ⎜ ⎟ ⎜ 0 0 1 0 1⎟ ⎜ 0 0 1 0 1⎟ 1 0 1 ⎟⎠ ⎜ 0 0 ⎝ ⎠ ⎝ ⎝ ⎠ ⎛1 0 0 7 0 ⎞⎟ ⎛⎜ 1 0 0 7 4 0 ⎞⎟ ⎜ 4 0⎟ ≈ ⎜ 0 − 4 0 1 0⎟ ≈ ⎜ 0 1 0 − 1 4 ⎟ ⎜ ⎟ ⎜ ⎜0 0 1 0 1⎟ ⎜0 0 1 0 1⎟ ⎝ ⎠ ⎝ ⎠ T
r ⎛7 1 ⎞ Sloupcový vektor ⎜ ,− ,0 ⎟ představuje řešení x0 nehomogenní soustavy a ⎝4 4 ⎠ r vektor (0,0,1)T řešení u příslušné homogenní soustavy. Obecné řešení nedourčené soustavy je ve tvaru T
r ⎛7 1 ⎞ x0 = ⎜ ,− ,0 ⎟ + α ⋅ (0,0,1)T . ⎝4 4 ⎠
Je zřejmé, že obě eliminační metody jsou metody přímé, tedy vedou ke zcela přesnému řešení. V obou uvedených případech jsme pracovali pro jednoduchost s čísly celými a bez zaokrouhlování. V praxi se však setkáváme se soustavami daleko většími a chceme-li pracovat na konečný počet desetinných míst, jsme nuceni zaokrouhlovat. Uvědomíme-li si,
42
kolik aritmetických operací zahrnuje proces eliminace, je na místě otázka, jak se malá chyba ve vstupních datech projeví na výstupu. Předpokládejme, že matice soustavy je A , která má rozklad A = LU . Vlivem ~ ~ zaokrouhlování však eliminací získáme matice L , U odpovídající jiné matici ~ ~~ ~ A = L U . Tedy získané řešení odpovídá soustavě s maticí A , nikoliv soustavě s maticí A . ~ ~ Zajímá nás rozdíl A − A . Předpokládejme, že matice U se liší od matice U o ~ matici chyb E , a podobně matice L od matice L o chybovou matici F , tedy ~ U =U +E ~ L = L+F Potom platí ~ ~~ A − A = LU − L U = LU − (U + E )(L + F ) = UF + LE + EF , což znamená, že pokud budou prvky matice L (multiplikátory) nebo prvky ~ matice U příliš velká čísla (co do absolutní hodnoty), bude rozdíl A − A mnohonásobně větší než vstupní chybové zatížení E nebo F a samotné řešení tudíž špatné. Za extrém lze považovat případ, kdy metodu v tomto tvaru nelze použít vůbec, a (k −1) = 0 , protože to když v některé fázi eliminace je diagonální prvek akk a (k −1) multiplikátor je mik = − ik(k −1) a víme, že nulou dělit nelze (můžeme to chápat i akk tak, že multiplikátor by byl nekonečno). V těchto případech říkáme, že algoritmus je nestabilní a jsme nuceni ho modifikovat. Eliminace s výběrem hlavního prvku. Prvek, jehož prostřednictvím vybíráme multiplikátor k-té fáze (a pomocí něhož nulujeme ostatní prvky v příslušném sloupci) nazýváme hlavní prvek (pivot) k-té fáze. Rovnice (příslušný řádek matice), z níž vybíráme hlavní prvek kté fáze nazýváme hlavní rovnice k-té fáze.
V Gaussově eliminační metodě popsané v předchozím odstavci jsme za hlavní prvek vybrali vždy automaticky prvek diagonální, prováděli jsme tzv. eliminaci bez výběru hlavního prvku. Nyní bude naší snahou vybírat hlavní prvek tak, aby multiplikátory byly co nejmenší čísla, tedy budou to prvky co do absolutní hodnoty největší. Výběr hlavního prvku můžeme provádět 1) v k-té fázi ze všech prvků matice, vzniklé vynecháním řádků a sloupců, které obsahovaly hlavní prvky předchozích fází, tzv. eliminace s úplným výběrem, nebo 2) pouze z některých prvků, tzv. eliminace s částečným výběrem a) řádkovým, pokud v k-té fázi eliminace vybíráme za hlavní prvek ten, který má největší absolutní hodnotu z těch prvků k-tého řádku patřícím
2. Řešení soustav lineárních rovnic zbývajícím n − k + 1 sloupcům, které do k-té fáze neobsahovaly hlavní prvek předchozí fáze b) sloupcovým, pokud v k-té fázi eliminace vybíráme za hlavní prvek ten, který má největší absolutní hodnotu z těch prvků k-tého sloupce patřícím zbývajícím n − k + 1 řádkům, které nebyly do k-té fáze hlavními řádky.
Příklad. Řešte Gaussovou eliminační metodou s výběrem a) řádkovým b) sloupcovým c) úplným soustavu x1 − x2 + x3 = −1
3 x1 + x2 − 3 x3 = 11 − x1 + 3 x2 − 2 x3 = 6 Řešení: a) V 1. fázi vybíráme prvek s největší absolutní hodnotou z prvků 1. řádku, tedy max( 1 , − 1 , 1 ) . Všechny prvky mají stejnou absolutní hodnotu, vybereme tedy ten, který má nejmenší sloupcový index, tj. prvek v pozici a11 = 1 a provedeme eliminaci ⎛ 1 −1 1 − 1⎞ ⎛ 1 − 1 1 − 1⎞ ⎜ ⎟ ⎜ ⎟ 11 ⎟ ≈ ⎜ 0 4 − 6 14 ⎟ ⎜ 3 1 −3 ⎜−1 3 − 2 6 ⎟⎠ ⎜⎝ 0 2 − 1 5 ⎟⎠ ⎝ Ve 2. fázi eliminace vybíráme s největší absolutní hodnotou z prvků 2. řádku, kromě 1. sloupce, protože již obsahoval v předchozí fázi hlavní prvek, tedy max( 4 , − 6 ) . Hlavní prvek bude prvek v pozici a23 = −6 . Provedeme eliminaci 2 ⎞ ⎛⎜ 1 ⎟ 14 ⎟ ≈ ⎜ 0 ⎜ 5 ⎟⎠ ⎜ 0 ⎝ Vhodným přeindexováním sloupců můžeme tvaru. Zpětnou substitucí dostaneme x2 = 2 , x3 = −1 , x1 = 2 . ⎛1 −1 1 ⎜ ⎜0 4 − 6 ⎜0 2 −1 ⎝
2 ⎞⎟ 4 −6 14 ⎟ 8 ⎟⎟ 4 0 3 3⎠ matici přepsat do trojúhelníkového −1
1
b) V 1. fázi vybíráme prvek s největší absolutní hodnotou z prvků 1. sloupce, tedy max ( 1 , 3 , −1 ) , tj. prvek v pozici a21 = 3 a provedeme eliminaci ⎛ 1 −1 1 − 1⎞ ⎛⎜ 0 − 4 2 − 14 ⎞⎟ 3 3 ⎜ ⎟ 1 11 ⎟ 11 ⎟ ≈ ⎜ 3 −3 ⎜ 3 1 −3 ⎜ ⎜−1 3 − 2 29 ⎟⎟ 6 ⎟⎠ ⎜ 0 10 −3 ⎝ 3 3 ⎠ ⎝ Ve 2. fázi eliminace vybíráme s největší absolutní hodnotou z prvků 2. sloupce, kromě 2. řádku, protože již obsahoval v předchozí fázi hlavní prvek, tedy
43
44
(
)
max − 4 , 10 . Hlavní prvek bude prvek v pozici a32 = 10 . Provedeme 3 3 3 eliminaci ⎛0 − 4 4 2 − 14 ⎞⎟ ⎛⎜ 0 0 − 4 ⎞⎟ ⎜ 3 3 5 5 ⎜3 1 11 ⎟ ≈ ⎜ 3 1 − 3 11 ⎟ −3 ⎜ 29 ⎟⎟ ⎜⎜ 0 10 29 ⎟⎟ ⎜ 0 10 −3 −3 3 3 ⎠ ⎝ 3 3⎠ ⎝ Vhodným přeindexováním řádků můžeme matici přepsat do trojúhelníkového tvaru. Zpětnou substitucí dostaneme x3 = −1 , x2 = 2 , x1 = 2 .
c) V 1. fázi vybíráme prvek s největší absolutní hodnotou ze všech prvků max( 1 , − 1 , 1 , 3 , 1 , − 3 , − 1 , 3 , − 2 ) , opět je několik prvků stejných. Vybereme ten s nejmenším indexováním, tj. prvek v pozici a21 = 3 a provedeme eliminaci ⎛ 1 −1 1 − 1⎞ ⎛⎜ 0 − 4 2 − 14 ⎞⎟ 3 3 ⎜ ⎟ ⎜ ⎟ − 3 1 3 11 11 ≈ − 3 1 3 ⎜ ⎟ ⎜ ⎜−1 3 − 2 29 ⎟⎟ 6 ⎟⎠ ⎜ 0 10 −3 ⎝ 3 3 ⎠ ⎝ Ve 2. fázi eliminace vybíráme s největší absolutní hodnotou z prvků matice, kromě 2. řádku a 1. sloupce, protože již obsahovaly v předchozí fázi hlavní prvek, tedy max − 4 , 10 , 2 , − 3 . Hlavní prvek bude prvek v pozici 3 3 a32 = 10 . Provedeme eliminaci 3 ⎛0 − 4 4 2 − 14 ⎞⎟ ⎛⎜ 0 0 − 4 ⎞⎟ ⎜ 3 3 5 5 ⎜3 1 11 ⎟ ≈ ⎜ 3 1 − 3 11 ⎟ −3 ⎜ 29 ⎟⎟ ⎜⎜ 0 10 29 ⎟⎟ ⎜ 0 10 −3 −3 3 3 ⎠ ⎝ 3 3⎠ ⎝ Zpětnou substitucí dostaneme řešení x3 = −1 , x2 = 2 , x1 = 2 .
(
)
Poznámka. Při rozhodování, zda použít Gaussovou eliminační metodu s výběrem hlavního prvku, nebo bez výběru mějme pořád na paměti, že pokud budeme hlavní prvek vybírat, potom kromě aritmetických operací samotné eliminace je počítač nucen provést srovnání prvků a výběr největšího co do absolutní hodnoty, a to v každém kroku eliminace, což je pro větší soustavy časově velice náročné. Je proto vždy potřeba zvážit, zda je nutný výběr prvku, a pokud ano, zda postačí pouze výběr řádkový, resp. sloupcový, nebo je-li nutný výběr úplný.
2. Řešení soustav lineárních rovnic Metoda LU-rozkladu. Již jsme načrtli, že velkou třídu matic lze jednoznačně rozložit na dolní trojúhelníkovou s jedničkami na diagonále a horní trojúhelníkovou matici. Této vlastnosti můžeme využít také k řešení soustav lineárních rovnic a to tak, že místo obecné soustavy r r Ax = b řešíme dvě trojúhelníkové soustavy r r r r Ly = b a Ux = y .
Víme, že A = LU a tak snadno odvodíme algoritmus pro získání prvků matic L,U : Vstup: n,
( )
A = aij
Opakuj pro j = 1, 2, ..., n Opakuj pro i = 1, 2, ..., j i −1
uij = aij − ∑ lir urj r =1
Opakuj pro i = j + 1, j + 2, ..., n
lij =
( )
( )
j −1 ⎞ 1 ⎛⎜ aij − ∑ lisusj ⎟ ⎟ u jj ⎜⎝ s =1 ⎠
Výstup: L = lij , U = uij .
Příklad. Řešte metodou LU-rozkladu soustavu lineárních rovnic x1 − x2 + x3 = −1
3 x1 + x2 − 3x3 = 11 − x1 + 3x2 − 2 x3 = 6
Řešení: Snadno ověříme, že LU rozklad matice A budou matice ⎛1 0 0 ⎞⎟ ⎛1 −1 1 ⎞ ⎜ ⎟ ⎜ L=⎜ 3 1 0⎟ a U = ⎜0 4 − 6⎟ . ⎜⎜ ⎟ ⎜0 0 −1 1 1⎟ 2 ⎟⎠ ⎝ 2 ⎝ ⎠ r r Nyní řešíme dvě trojúhelníkové soustavy, nejdříve Ly = b . Řešení obdržíme r r r přímou substitucí y = (− 1, 14, − 2 )T . Řešíme druhou soustavu Ux = y , zpětnou r substitucí dostaneme řešení původní soustavy x = (2, 2, − 1)T .
45
46 Choleskeho algoritmus. Připomeňme, že matice A se nazývá a) symetrická právě tehdy když AT = A , tedy ∀i, j : aij = a ji a
b) pozitivně definitní právě tehdy když pro r rT r n n T x = ( x1, ..., xn ) platí x Ax = ∑∑ aij xi x j > 0 .
každý
nenulový
vektor
i =1 j =1
Pro symetrickou, pozitivně definitní matici platí, že existuje právě jedna horní trojúhelníková matice U s kladnými diagonálními prvky, že A = U T U . Tedy nehledáme v rozkladu dvě matice L,U , ale pouze jednu matici U . Celý algoritmus metody LU-rozkladu se tím značně zjednoduší, tomuto speciálnímu upravení metody LU-rozkladu pro symetrické, pozitivně definitní matice se říká Choleskeho algoritmus.
3.3 Metody iterační. Jacobiova metoda. Základní myšlenku iteračních metod jsme popsali v kapitole 2. Iterační metody určené pro řešení soustav lineárních rovnic kterým se budeme věnovat r r r r r spočívají v převedení soustavy Ax = b do tvaru x = Hx + q a poté určení iteračního předpisu.
Jacobiova metoda samotný „převod“ realizuje tak, že z i-té rovnice vyjádří i-tou neznámou, tj. ⎞ ⎛ ⎟ n 1 ⎜ ⎜ bi − ∑ aij x j ⎟ pro každé i = 1, 2, ..., n xi = aii ⎜ ⎟ j =1 ⎟ ⎜ j ≠i ⎠ ⎝ r (k +1) r (k ) r = Hx + q , tedy a iterační předpis určí jako x ⎞ ⎛ ⎟ ⎜ n 1 ( ( k +1) k)⎟ ⎜ pro každé i = 1, 2, ..., n . xi bi − ∑ aij x j = aii ⎜ ⎟ 1 j = ⎟ ⎜ j ≠i ⎠ ⎝ T r Je ovšem na místě otázka, jak volit počátečnou aproximaci x (0 ) = x1(0 ),..., xn(0 ) aby metoda konvergovala. Platí následující věta.
(
)
Věta. r r r Iterační posloupnost x (k +1) = Hx (k ) + q konverguje pro libovolnou počátečnou r aproximaci x (0 ) právě tehdy, když spektrální poloměr matice H je ostře menší než 1, tj. ρ (H ) = max λi (H ) < 1 , kde λi (H ) jsou vlastní čísla matice H .
i
2. Řešení soustav lineárních rovnic
Je-li matice A ostře diagonálně dominantní (tj. aii >
47
n
∑ aik
< 1 pro každé
k =1 k ≠i
i = 1,..., n ),
potom Jacobiova metoda konverguje pro libovolnou počáteční r (0 ) aproximaci x .
Příklad. Řešte Jacobiovou metodou soustavu lineárních rovnic 5 x1 + x3 − x4 = 3
5 x2 − x3 + x4 = −3 x1 − x2 + 5 x3 = 2 − x1 + x2 + 5 x4 = 8 Řešení: Matice soustavy je ostře diagonálně dominantní, metoda by tudíž měla r konvergovat pro libovolný počáteční vektor x (0 ) . Pro jednoduchost volme r x (0 ) = (0 ,0 ,0 ,0)T . r Snadno můžeme ověřit, že přesné řešení této soustavy je x = ( 1, − 1, 0, 2 )T . Rekurentní předpis iterace bude 1 x1(k +1) = 3 − x3(k ) + x4(k ) 5 1 x2(k +1) = − 3 + x3(k ) − x4(k ) 5 1 x3(k +1) = 2 − x1(k ) + x2(k ) 5 1 x4(k +1) = 8 + x1(k ) − x2(k ) 5 a jednotlivé hodnoty jsou seřazeny v tabulce.
( ( ( (
Iterace k x1(k ) 0 1 2 3 4 5
0 0.6 0.84 0.936 0.9744 0.98976
)
)
) )
r Složky vektoru x (k ) x2(k ) x3(k ) 0 -0.6 -0.84 -0.936 -0.9744 -0.98976
0 0.4 0.16 0.064 0.0256 0.01024
x4(k ) 0 1.6 1.84 1.936 1.9744 1.98976
Z tabulky vidíme, že jednotlivé iterace složek vektoru řešení skutečně konvergují k přesnému řešení.
48 Gauss-Seidelova metoda. Vylepšením Jacobiovy metody je metoda Gauss-Seidelova, která je postupově stejná, ale využívá nově spočtené iterace ihned k dalšímu výpočtu. Tedy převod soustavy zůstává stejný- z i-té rovnice vyjadřujeme i-tý prvek, ale iterační předpis je pozměněn i −1 n ⎞ 1 ⎛⎜ xi(k +1) = bi − ∑ aij x (jk +1) − ∑ aij x (jk ) ⎟ pro každé i = 1, 2, ..., n . ⎟ aii ⎜ j =1 j = i +1 ⎠ ⎝
Při výpočtu i-té složky (k+1)-iterace xi(k +1) jsou již využity předchozí spočtené složky (až do (i-1)-složky, tedy x (k +1),..., x (k +1) -protože ty jsou již známé), 1
i −1
zatímco ostatní ještě neznáme, tudíž použijeme předchozí k-té iterace, tedy xi(+k1) ,..., xn(k ) . Vzhledem k tomu, že každá nová iterace je přesnější než předchozí (samozřejmě za předpokladu, že proces konverguje), je metoda Gauss-Seidelova o něco rychlejší. Postačující podmínku pro konvergenci Gauss-Seidelovy metody vystihuje následující věta. Věta. Je-li matice A ostře diagonálně dominantní, potom Gauss-Seidelova metoda r konverguje pro libovolnou počátečnou aproximaci x (0 ) . Je-li matice A symetrická a pozitivně definitní, potom Gauss-Seidelova metoda r konverguje pro libovolnou počátečnou aproximaci x (0 ) . Poznamenejme, že pro některé obecné matice A může Jacobiova metoda konvergovat a Gauss-Seidelova metoda konvergovat nemusí, nebo obráceně. Příklad. Řešte Gauss-Seidelovou metodou soustavu lineárních rovnic 5 x1 + x3 − x4 = 3
5 x2 − x3 + x4 = −3 x1 − x2 + 5 x3 = 2 − x1 + x2 + 5 x4 = 8 Řešení: Rekurentní předpis iterace bude 1 x1(k +1) = 3 − x3(k ) + x4(k ) 5 1 x2(k +1) = − 3 + x3(k ) − x4(k ) 5 1 x3(k +1) = 2 − x1(k +1) + x2(k +1) 5 1 x4(k +1) = 8 + x1(k +1) − x2(k +1) 5
( ( ( (
)
)
) )
2. Řešení soustav lineárních rovnic
49
r a jednotlivé hodnoty pro vektor počátečné iterace x (0 ) = (0 ,0 ,0 ,0)T jsou seřazeny v tabulce.
Iterace k x1(k ) 0 1 2 3 4 5
0 0.6 0.936 0.98976 0.9983616 0.999737856
r Složky vektoru x (k ) x2(k ) x3(k ) 0 -0.6 -0.936 -0.98976 -0.9983616 -0.999737856
0 0.16 0.0256 0.004096 0.00065536 0.0001048576
x4(k ) 0 1.84 1.9744 1.995904 1.99934464 1.999895142
Jednotlivé iterace složek vektoru řešení skutečně konvergují k přesnému řešení, matice soustavy byla ostře diagonálně dominantní. I když uvedené věty popisují postačující podmínky konvergence těchto iteračních metod, neříkají nic o rychlosti konvergence. Pro některé matice, přesto že proces konverguje, je velice pomalý, proto je nutné zrychlit ho. Existuje celá řada způsobů jak zrychlení provést. O některých často užívaných, tzv. Aitkenově extrapolační formuli, nebo relaxační metodě, se může čtenář dovědět více např. v RALSTON, A.: Základy numerické matematiky. Praha, Academia, 1978.
3.4 Kontrolní otázky, korespondenční úkol. Kontrolní otázky. 1. Řešte Gaussovou eliminační metodou a) bez výběru b) s výběrem hlavního prvku soustavu 5.2468 x1 − 0.1254 x2 + 2.24296 x3 + 17.3698 x4 = 1.2202
15.4523 x1 + 1.36790 x2 + 7.5241x3 − 2.1511x4 = −3.6541 − 3.2987 x1 + 11.2445 x2 − 0.1167 x3 + 4.5632 x4 = −12.8847 12.2544 x1 − 1.0023 x2 + 4.7758 x3 + 6.2541x4 = 6.1478 Výsledky srovnejte a rozhodněte, zda je eliminace s výběrem prvku nutná, nebo nikoliv. 2. Řešte a) Gaussovou eliminační metodou bez výběru b) s výběrem hlavního prvku (řádkovým, sloupcovým i úplným) c) řešte metodou LU-rozkladu soustavu
50
5.0000 x1 − 0.1200 x2 + 0.0020 x3 = 1.2202 15.000 x1 + 7.0000 x2 − 2.0000 x3 = −3.6541 − 0.0002 x1 + 11.0000 x2 − 0.1167 x3 = −12.8847 Výsledky srovnejte a rozhodněte, zda je eliminace s výběrem prvku nutná, nebo nikoliv. 3. Užijte Gauss-Jordanovy metody k nalezení inverze matice ⎛ 5.0000 − 0.1200 0.0020 ⎞ ⎜ ⎟ ⎜ 15.0000 7.0000 − 2.0000 ⎟ ⎜ − 0.0002 11.0000 − 0.1167 ⎟ ⎝ ⎠ 4. Diskutujte modifikaci metody LU-rozkladu pro třídiagonální matici. 5. Řešte Jacobiovou a Gauss-Seidelovou metodou soustavu lineárních rovnic s rozšířenou maticí soustavy ⎛ 13.50 3.21 ⎞ 1.20 2.78 − 1.36 ⎟ ⎜ − 1.12 ⎟ − 1.12 1.32 5.89 ⎜ 0 ⎜ 0.56 − 11.25 28.15 − 2.74 ⎟ 0 ⎟ ⎜ ⎟ ⎜ − 3.24 1.65 18 . 50 0 . 32 8 . 29 ⎠ ⎝ Výsledky srovnejte. Diskutujte konvergenci metod. 6. Řešte Jacobiovou a Gauss-Seidelovou metodou soustavu lineárních rovnic s rozšířenou maticí soustavy ⎛5 2 0 − 3.5 0 ⎞ ⎟ ⎜ − 3⎟ 1 ⎜0 2 −1 ⎜0 − 6 7 − 7⎟ 0 ⎟ ⎜ ⎟ ⎜1 1 − 3 8 7 ⎠ ⎝ Výsledky srovnejte. Diskutujte konvergenci metod. Část pro zájemce. 1. Upravte algoritmus metody LU-rozkladu pro symetrickou, pozitivně definitní matici. Využijte faktu, že existuje právě jedna horní trojúhelníková matice U
s kladnými diagonálními prvky, že A = U T U . 2. Diskutujte použití Gaussovy eliminační metody k výpočtu determinantu matic. Určete hodnotu determinantu matice ⎛ 7.4 2.2 − 3.1 0.7 ⎞ ⎜ ⎟ ⎜ 1.6 4.8 − 8.5 4.5 ⎟ ⎜ 4.7 7.0 − 6.0 6.6 ⎟ ⎜ ⎟ ⎜ 5.9 2.7 4.9 − 5.3 ⎟ ⎝ ⎠
2. Řešení soustav lineárních rovnic
Shrnutí kapitoly. Jste na konci kapitoly Řešení soustav lineárních rovnic. Po prostudování by jste měli být schopni popsat metody přímé a iterační, případně je modifikovat pro použití ve speciálních případech (speciální typy matic soustavy, řešení maticových rovnic, hledání inverze, řešení nedourčených soustav, atd.). Pokud jste zvládli zpracovat algoritmy všech popsaných metod a vyřešili kontrolní otázky, pokuste se samostatně vypracovat korespondenční úkol.
Korespondenční úkol. Řešte soustavu lineárních rovnic s rozšířenou maticí soustavy ⎛ 15.002 2.133 0.001 − 3.556 0.210 ⎞ ⎟ ⎜ ⎜ 0.115 12.337 − 1.001 1.552 − 3.120 ⎟ ⎜ − 0.025 − 6.550 17.998 0.774 − 7.740 ⎟ ⎟ ⎜ ⎟ ⎜ 1.120 − 18 . 500 1 . 120 3 . 654 17 . 889 ⎠ ⎝ a) Jacobiovou a Gauss-Seidelovou metodou b) Gaussovou eliminační metodou (rozhodněte, zda je nutno použít výběr prvku, či nikoliv) c) Metodou LU-rozkladu. d) Popište použití Gauss-Jordanovy metody k nalezení inverze matice soustavy a určete ji. Opět diskutujte nutnost výběru hlavního prvku.
51
Vlastní čísla a vlastní vektory matic
4. VLASTNÍ ČÍSLA A VLASTNÍ VEKTORY MATIC V této kapitole se dozvíte:
•
co je to vlastní číslo matice a numerické možnosti jeho výpočtu;
•
co je to tzv. částečný problém vlastních čísel a jaké numerické metody lze použít k hledání speciálních vlastních čísel (třeba největšího nebo nejmenšího);
•
co je tzv. úplný problém vlastních čísel, jak využít charakteristického polynomu matice a podobnosti matic k numerickému nalezení všech vlastních čísel matice.
Budete schopni:
•
numericky stanovit dominantní vlastní číslo;
•
určit všechna vlastní čísla užitím LU-rozkladu matice;
•
určit všechna vlastní čísla metodami ortogonálních transformací.
Klíčová slova této kapitoly: Vlastní číslo, vlastní vektor matice, charakteristický polynom matice, podobná matice, ortogonální matice. Čas potřebný k prostudování učiva kapitoly: 10 + 10 hodin (teorie + řešení úloh)
Průvodce studiem. Tato kapitola je věnována numerickému hledání vlastních čísel matic. Je rozdělena na dvě podkapitoly: Částečný problém vlastních čísel a Úplný problém vlastních čísel. Jednotlivé metody jsou demonstrovány na konkrétních příkladech. Tato kapitola rovněž vyžaduje alespoň elementární znalosti z teorie matic.
4.1 Formulace problému.
Nechť A je čtvercová matice řádu n . Číslo λ , pro které má homogenní soustava r r Av = λv nenulové řešení se nazývá vlastní číslo matice A a jemu příslušné nenulové řešení r v vlastní vektor matice A . Soustavu je možné přepsat do tvaru ( A − λI ) vr = or , r kde I je jednotková matice řádu n , a o je nulový vektor. Z algebry víme, že taková soustava má nenulové řešení právě tehdy, když det ( A − λI ) = 0 .
53
54
Tuto rovnici, která představuje polynom n-tého stupně vzhledem k neznámé λ , nazýváme charakteristický polynom matice A . Každé vlastní číslo je tedy kořenem charakteristického polynomu, tudíž čtvercová matice řádu n má právě n vlastních čísel. Navíc jejich polohu můžeme dopředu odhadnout. Věta. Nechť A = aij
( )
je čtvercová matice řádu n , potom její vlastní čísla leží v
komplexní rovině ve sjednocení kruhů z − aii ≤
n
∑
j =1 j ≠i
aij .
Pro nejmenší, co do absolutní hodnoty, navíc platí odhad λmin ≤ n ( r1 ⋅ r2 ⋅ ... ⋅ rn ) , kde 1 2 ⎞2
⎛ ri = ⎜ ∑ aij ⎟ . ⎟ ⎜ j =1 ⎠ ⎝ n
V algebře je někdy postačující určit pouze některé vlastní číslo (nejmenší, největší), tehdy mluvíme o tzv. částečném problému vlastních čísel. Úloha naleznout všechna vlastní čísla dané matice se nazývá úplný problém vlastních čísel.
4.2 Částečný problém vlastních čísel. Mocninná metoda. Mocninná metoda je určená k nalezení tzv. dominantního vlastního čísla, tedy vlastního čísla s největší absolutní hodnotou. Předpokládejme, že matice A má jediné dominantní vlastní číslo a n lineárně nezávislých vlastních vektorů. r Libovolný vektor y (0 ) můžeme zapsat jako lineární kombinaci vlastních vektorů r r r y (0 ) = α1v1 + ... + α nvn r ∞ a sestrojme posloupnost y (k ) k =0 zadanou rekurencí r r y (k ) = A y (k −1) .
{ }
r r r r Protože y (k ) = Ak y (0 ) a Avi = λi vi , bude r r r r r y (k ) = α1 Ak v1 + ... + α n Ak vn = α1λ1k v1 + ... + α n λkn vn .
Označíme-li λ1 dominantní vlastní číslo (které je jediné), bude každé i = 2,..., n a tedy
λi < 1 pro λ1
Vlastní čísla a vlastní vektory matic
55
⎛ n ⎛ λ ⎞k r ⎞ lim ⎜ ∑α i ⎜⎜ i ⎟⎟ vi ⎟ = 0 ⎟ k →∞⎜ i = 2 ⎝ λ1 ⎠ ⎝ ⎠ ( k +1) a y (jk ) pro každé j = 1,..., n budou což znamená, že podíly složek vektorů y j konvergovat k hledanému dominantnímu vlastnímu číslu, tj. y (jk +1) λ1 = lim (k ) . k →∞ y j
V každém kroku tak dostaneme n aproximací téhož vlastního čísla. Z praktického hlediska je proto rozumné za aproximaci považovat jejich aritmetický průměr. Algoritmus mocninné metody: r Vstup: n , A = aij , y (0 )
( )
Opakuj pro k = 1,..., N r r y (k ) = A y (k −1)
y2(k ) y1(k ) yn(k ) ( ) ( ) , λ = , …, λ = 1 n 1 2 y1(k −1) y2(k −1) yn(k −1) (λ ) + (λ1 ) 2 + ... + (λ1 ) n λ1 : = 1 1 n
(λ1 ) 1=
Výstup: λ1
Příklad. Nalezněte mocninnou metodou dominantní vlastní číslo matice 3 ⎞ ⎛ 5.5 − 4 ⎟ ⎜ ⎜ − 4 10.5 − 7 ⎟ ⎜ 3 − 7 10.5 ⎟ ⎠ ⎝
Řešení: r Pro jednoduchost zvolíme počáteční vektor y (0 ) = ( 1,1,1)T . Jednotlivé složky vektorů posloupnosti, aproximace vlastního čísla a jejich aritmetický průměr jsou uvedeny v tabulce. k
y1(k )
y2(k )
y3(k )
(λ1 ) 1
(λ1 ) 2
(λ1 ) 3
λ1
0 1 2 3 4 5
1 4.5 46.25 785.125 14878.0625 286335.7812
1 -0.5 -68.75 -1503.625 -29534.4375 -571150.7812
1 6.5 85.25 1515.125 28789.5625 553665.6562
4.5000 10.2778 16.9757 18.9477 19.2455
-0.5000 137.500 21.8709 19.6422 19.3385
6.5000 13.1154 17.7727 19.0014 19.2315
3.5000 53.6311 18.8731 19.1978 19.2718
Za aproximaci dominantního vlastního čísla můžeme považovat λ1 = 19.2718 .
56 Poznámka. Jsou-li zaručeny předpoklady, bude metoda vždy konvergovat pro libovolně r zvolený počáteční vektor y (0 ) . Nicméně v praxi nemáme předem žádné informace o tom, zda existuje jediné dominantní vlastní číslo nebo zda jsou vlastní vektory lineárně nezávislé, proto se může stát, že metoda konvergovat nebude, nebo bude proces konvergence velice pomalý. Způsoby jak lze proces urychlit můžeme nalézt např. v RALSTON, A.: Základy numerické matematiky. Praha, Academia, 1978.
Metoda Rayleighova podílu. Předpokládejme, že matice A je reálná, symetrická, má jediné dominantní vlastní číslo a n lineárně nezávislých vlastních vektorů. Konstrukce metody Rayleighova podílu je stejná jako metody mocninné, avšak využívá toho, že takováto matice bude mít vlastní vektory ortonormální. Za těchto podmínek lze dominantní vlastní číslo aproximovat jako r T r r Tr y (k ) A y (k ) y (k ) y (k +1) λ1 = lim = lim k → ∞ r (k )T r (k ) k → ∞ r (k )T r (k ) y y y y Metoda Rayleighova podílu zpravidla konverguje rychleji, než metoda mocninná.
Příklad. Nalezněte metodou Rayleighova podílu dominantní vlastní číslo matice 3 ⎞ ⎛ 5.5 − 4 ⎟ ⎜ ⎜ − 4 10.5 − 7 ⎟ ⎜ 3 − 7 10.5 ⎟ ⎠ ⎝
Řešení: r Zvolme počáteční vektor stejně jako v příkladu předchozím y (0 ) = ( 1,1,1)T . r r Výpočet vektorů posloupnosti y (k ) = A y (k −1) zůstává stejný jako u metody mocninné, liší se pouze určení aproximace dominantního vlastního čísla. V r Tr y (k ) y (k +1) . každém kroku aproximujeme λ1 = r (k )T r (k ) y y Hodnoty jsou seřazeny v tabulce. k 0 1 2 3 4 5
r y (k ) (1,1,1) (4.5,-0.5,6.5) (46.25,-68.75,85.25) (785.125,-1503.625,151.125) (14878.0625,-29534.4375,28789.5625) (286335.7812,-571150.7812,553665.6562)
λ1 3.5000 12.6952 19.0226 19.2753 19.2816
Vlastní čísla a vlastní vektory matic
Za aproximaci dominantního vlastního čísla můžeme považovat λ1 = 19.2816 . V porovnání s metodou mocninnou vidíme, že proces konvergence je podstatně rychlejší.
4.3 Úplný problém vlastních čísel.
Metody úplného problému vlastních čísel můžeme rozdělit do dvou kategorií na: 1. metody založené na výpočtu pomocí charakteristického polynomu 2. metody založené na podobnosti matic. První metody plynou z faktu, že vlastní číslo je kořenem charakteristického polynomu det ( A − λI ) = 0 , pro velká n však nastávají komplikace jednak se samotným výpočtem determinantu a poté s výpočtem kořenů. Existují metody, které určují koeficienty charakteristického polynomu iteračně, více o těchto metodách lze nalézt např. v RALSTON, A.: Základy numerické matematiky. Praha, Academia, 1978. Druhá kategorie metod je založená na podobnosti matic. Připomeňme, že dvě matice A, B nazveme podobné existuje-li regulární matice P taková, že P −1 AP = B , respektive A = PBP −1 .
Protože
(
)
det ( A − λI ) = det PBP −1 − λI = det P(B − λI )P −1 = det (B − λI ) tak platí, že podobné matice mají stejná vlastní čísla. (Pozor, z rovnosti vlastních čísel neplyne podobnost matic).
Metoda LU-rozkladu. Hledejme vlastní čísla matice A a předpokládejme, že lze sestrojit její LUrozklad A = LU . Potom matice B = UL bude s maticí A podobná, protože
U = L−1 A a tedy B = L−1 AL . Sestrojme posloupnost matic {Ak } tak, že A0 = A = L0U 0 a Ak = U k −1Lk −1 = LkU k . Dá se ukázat, že pokud součiny L1L2 ⋅... ⋅ Lk pro k → ∞ konvergují k regulární matici, tak matice Ak konvergují k horní trojúhelníkové matici s vlastními čísly na diagonále. Algoritmus metody: Vstup: A A0 : = A Pro k = 0,1, 2,...
57
58
sestroj Ak = LkU k Ak +1 : = U k Lk Výstup: Ak Poznámka. Je-li matice symetrická a pozitivně definitní, použijeme k nalezení LU-rozkladu Choleskeho algoritmus. Příklad. Určete metodou LU-rozkladu všechna vlastní čísla matice 3 ⎞ ⎛ 5.5 − 4 ⎟ ⎜ A = ⎜ − 4 10.5 − 7 ⎟ ⎜ 3 − 7 10.5 ⎟ ⎠ ⎝
Řešení: V tomto případě použijeme k sestrojování LU-rozkladu Choleskeho algoritmus, což je speciálně upravený algoritmus pro symetrické a pozitivně definitní matice. (Podrobně lze najít např. v [6]). Jednotlivé matice posloupnosti {Ak } jsou seřazeny v tabulce. k
Ak
0
3 ⎞ ⎛ 5.5 − 4 ⎟ ⎜ ⎜ − 4 10.5 − 7 ⎟ ⎜ 3 − 7 10.5 ⎟ ⎠ ⎝
1
⎛ 10.045455 − 6.936270 3.082162 ⎞ ⎟ ⎜ ⎜ − 6.936270 10.649156 − 4.213587 ⎟ ⎜ 3.082162 − 4.213587 5.805389 ⎟ ⎠ ⎝
2
⎛ 15.780543 − 6.135374 1.973289 ⎞ ⎟ ⎜ ⎜ − 6.135374 6.601899 − 1.748106 ⎟ ⎜ 1.973289 − 1.748106 4.117558 ⎟ ⎠ ⎝
... 14
... ⎛ 19.281803 − 0.000847 0.000042 ⎞ ⎟ ⎜ ⎜ − 0.000847 4.266010 − 0.154657 ⎟ ⎜ 0.000042 − 0.154657 2.952187 ⎟ ⎠ ⎝
15
⎛ 19.281803 − 0.000399 0.000017 ⎞ ⎟ ⎜ ⎜ − 0.000399 4.271617 − 0.128534 ⎟ ⎜ 0.000017 − 0.128534 2.946580 ⎟ ⎠ ⎝
Za aproximace vlastních čísel můžeme považovat diagonální prvky, tedy λ1 = 19.281803 , λ2 = 4.271617 , λ3 = 2.946580 .
Vlastní čísla a vlastní vektory matic Protože matice mají konvergovat k horní trojúhelníkové matici (v našem případě z důvodu symetrie přímo k diagonální matici), v praxi zpravidla zastavujeme iterační proces ve chvíli, kdy největší prvek pod diagonálou je menší než předem zadaná tolerance (tudíž ho se zadanou tolerancí můžeme považovat za nulový). Z tabulky vidíme, že proces konverguje velice pomalu a uvědomíme-li si, kolik operací je nutno provést v každém kroku, tak pro velké matice (velká n ) je použití metody dost nevýhodné. Pro obecný typ matice se navíc proces může chovat nestabilně.
Metody ortogonálních transformací. Jacobiova diagonalizace. U metod ortogonálních transformací budeme opět sestrojovat posloupnost navzájem podobných matic {Ak } a to tak, že
A0 = A a Ak +1 = QkT Ak Qk , kde Qk jsou ortogonální matice speciálně vybírané tak, aby tato posloupnost konvergovala k matici, pro kterou je snadné určit vlastní čísla (např. trojúhelníkové, nebo diagonální). Připomeňme, že reálná ortogonální matice je matice, jejíž sloupce jsou ortonormální vektory, resp. pro kterou QT = Q −1 . Popíšeme metodu určenou pro symetrické matice- Jacobiovu diagonalizaci, která vychází z následující věty. Věta. Nechť A je reálná symetrická matice. Potom existuje ortogonální matice Q taková, že QT AQ = Λ je diagonální matice s vlastními čísly matice A na diagonále. Matici Q sestrojujeme pomocí posloupnosti matic Q pq (α ) , kde ⎛1 ⎞ ⎜ ⎟ 0 0 ⎟ ⎜ 1 ⎜ ⎟ 1 ⎜ ⎟ c ... − s ⎜ ⎟ Q pq (α ) = ⎜ 0 ... 0 ⎟ ⎜ ⎟ s ... c ⎜ ⎟ 1 ⎜ ⎟ ⎜ 0 0 1 ⎟⎟ ⎜ ⎜ 1 ⎟⎠ ⎝ je matice lišící se od jednotkové pouze v pozicích ( p, p ) , ( p, q ) , (q, p ) , (q, q ) přičemž c = cosα a s = sin α . Úhel α vybíráme tak, abychom nulovali prvek b pq = bqp matice B = QTpq (α ) AQ pq (α ) .
59
60
Požadavek b pq = bqp = 0 vede k podmínce tg 2α =
2a pq a pp − aqq
. Odtud pak
dostaneme
1 a pp − aqq + 2 2r a 1 pp − aqq s 2 = sin 2 α = − 2 2r a pq sc = r 2 r = a pp − aqq 2 + 4a 2pq .
c 2 = cos 2 α =
(
)
Při konkrétním výpočtu není nutno počítat úhel α , ale lze odvodit přímý vztah pro výpočet prvků matice B (volíme α ≤
bpp = bqq =
π
4
):
a pp + aqq + r
(a
2 a − a 2pq
pp qq
bpp
)=b
pp
−r
bip = aip c + aiq s = b pi biq = − aip s + aiq c = bqi bij = aij v ostatních pozicích. Volbu znaménka u sin α a cos α provádíme tak, aby sign sin α = sign a pq .
Celý algoritmus Jacobiovy diagonalizace probíhá následovně: Vstup: první člen posloupnosti A0 = A Pro každé k = 0,1, 2,... • určíme nediagonální pozici ( p, q ) kterou chceme nulovat. Můžeme nulovat vždy největší prvek, co do absolutní hodnoty (algoritmicky náročnější, neboť počítač musí největší prvek naleznout), nebo nulujeme postupně. • stanovíme prvky další matice posloupnosti, matice B dle uvedených vztahů • přeznačíme Ak +1 = B a pokračujeme pro další k Výstup: diagonální matice Ak s vlastními čísly matice A na diagonále. Prvky, které jsme nulovali v předchozích krocích obecně nezůstanou nulové, ale dá se ukázat, že se zmenšují. Proces zastavujeme ve chvíli kdy největší nediagonální prvek je menší než předem daná tolerance. Metody ortogonálních transformací se chovají na rozdíl od metody LU-rozkladu numericky stabilně.
Vlastní čísla a vlastní vektory matic
Příklad. Určete Jacobiovou metodou všechna vlastní čísla matice 3 ⎞ ⎛ 5.5 − 4 ⎟ ⎜ A = ⎜ − 4 10.5 − 7 ⎟ ⎜ 3 − 7 10.5 ⎟ ⎠ ⎝
Řešení: Budeme nulovat vždy největší (co do absolutní hodnoty) nediagonální prvek. V prvním kroku bude tedy pozice ( p, q ) = (2,3) , neboť prvek a23 = −7 je v absolutní hodnotě největší nediagonální prvek. Pozice, které nulujeme, stejně jak hodnoty r , c, s a tvary matic Ak jsou uvedeny v tabulce. matice Ak
3 ⎞ ⎛ 5.5 − 4 ⎟ ⎜ A0 = ⎜ − 4 10.5 − 7 ⎟ ⎜ 3 − 7 10.5 ⎟ ⎠ ⎝
krok k 0 matice Ak
c s pozice ( p, q ) r (2,3) 14 0.707107 -0.707107 5.5 − 4.949747 − 0.707107 ⎞ ⎛ ⎟ ⎜ A1 = ⎜ − 4.949747 17.5 0 ⎟ ⎟ ⎜ − 0.707107 0 3.5 ⎠ ⎝
krok k 1 matice Ak krok k 2 matice Ak ... matice Ak
c s pozice ( p, q ) r (1,2) 15.556349 0.338091 -0.941114 0 − 0.239066 ⎞ ⎛ 19.278174 ⎟ ⎜ A2 = ⎜ 0 3.721826 − 0.665467 ⎟ ⎟ ⎜ − 0.239066 − 0.665467 3.5 ⎠ ⎝ c s pozice ( p, q ) r (2,3) 1.349293 0.763021 -0.646374 ⎛ 19.278174 0.154526 − 0.182412 ⎞ ⎟ ⎜ A3 = ⎜ 0.154526 4.285560 0 ⎟ ⎜ − 0.182412 0 2.936266 ⎟⎠ ⎝ 0 − 0.000018 ⎞ ⎛ 19.281802 ⎟ ⎜ A5 = ⎜ 0 4.283968 0.001725 ⎟ ⎜ − 0.000018 0.001725 2.934230 ⎟ ⎠ ⎝
61
62
Za aproximaci vlastních čísel považujeme diagonální prvky λ1 = 19.281802 , λ1 = 4.283968 , λ1 = 2.934230 . Vidíme, že Jacobiova diagonalizace je podstatně rychlejší než metoda LU-rozkladu.
4.4 Kontrolní otázky, korespondenční úkol. Kontrolní otázky.
1. Odhadněte polohu vlastních čísel a nalezněte mocninnou metodou dominantní vlastní číslo matice ⎛ 0 5 0⎞ ⎟ ⎜ ⎜ − 2 7 2⎟ . ⎜ −1 5 3⎟ ⎠ ⎝ Dále určete všechna vlastní čísla matice jako kořeny charakteristického polynomu (pokud neumíte řešit rovnici 3. stupně, použijte k řešení některé z metod pro řešení nelineárních rovnic). Určete příslušné vlastní vektory jako r r řešení homogenních soustav ( A − λi I ) v = o , srovnejte vlastní vektor příslušný r dominantnímu vlastnímu číslu s vektory y (k ) iterační posloupnosti mocninné metody. 2. Stanovte dominantní vlastní číslo matice ⎛10.29 2.36 0.02 ⎞ ⎟ ⎜ ⎜ 1.56 2.13 0.21 ⎟ . ⎜ 7.15 2.54 10.11⎟ ⎠ ⎝ Určete vlastní čísla metodou LU-rozkladu. Dále určete všechna vlastní čísla matice jako kořeny charakteristického polynomu a určete příslušné vlastní vektory. Výsledek srovnejte. 3. Odhadněte polohu vlastních čísel a stanovte mocninnou metodou a metodou Rayleighova podílu dominantní vlastní číslo matice ⎛ 6 2 2⎞ ⎟ ⎜ ⎜ 2 6 2⎟ . ⎜ 2 2 6⎟ ⎠ ⎝ Určete vlastní čísla metodou LU-rozkladu a určete všechna vlastní čísla matice jako kořeny charakteristického polynomu. Výsledek srovnejte. 4. Určete všechna vlastní čísla matice ⎛6 4 1⎞ ⎟ ⎜ ⎜ 4 6 4⎟ ⎜1 4 6⎟ ⎠ ⎝ jako kořeny charakteristického polynomu. Stanovte vlastní čísla Jacobiovou diagonalizací a metodu LU-rozkladu. Výsledky srovnejte.
Vlastní čísla a vlastní vektory matic
Část pro zájemce.
1. Pokuste se poupravit mocninnou metodu tak, aby hledala nejmenší vlastní číslo (co do absolutní hodnoty). 2. Diskutujte Jacobiovu diagonalizaci pro matice typu (2,2 ) . Jedná se o iterační proces? 3. Jak lze u mocninné metody a metody Rayleighova podílu aproximovat vlastní r vektory užitím vektorů y (k ) konstruované iterační posloupnosti?
Shrnutí kapitoly. Dostáváte se na konec poslední kapitoly studijní opory. Ukázali jsme si, jak je možné vypočítat vlastní čísla matic numerickými metodami. Pokud jste zvládli zpracovat algoritmy všech popsaných metod a vyřešili kontrolní otázky, pokuste se samostatně vypracovat korespondenční úkol.
Korespondenční úkol. Odhadněte polohu vlastních čísel matice ⎛ 7.25 − 2.06 5.12 ⎞ ⎟ ⎜ A = ⎜ − 2.06 7.58 2.13 ⎟ ⎜ 5.12 2.13 12.50 ⎟⎠ ⎝
a určete vlastní čísla této matice jako kořeny charakteristického polynomu. (Pokud neumíte řešit rovnici 3. stupně, použijte k řešení některé z metod pro řešení nelineárních rovnic.) Nalezněte dominantní vlastní číslo této matice metodou mocninnou a metodou Rayleighova podílu a obě metody srovnejte. Nalezněte všechna vlastní čísla metodou LU-rozkladu a metodou Jacobiovy diagonalizace. Srovnejte rychlost konvergence obou metod.
63
Řešení úloh
ŘEŠENÍ ÚLOH Řešení kontrolních otázek. Část 1.4
Příklad 1: a) Určit všechny kořeny algebraické rovnice 5. stupně a5 x 5 + a4 x 4 + a3 x3 + a2 x 2 + a1x + a0 = 0 je úloha matematická i numerická, neboť je dána konečným počtem vstupů a výstupů. Vstupní hodnoty jsou koeficienty rovnice (5 hodnot a5 = 1, a4 = 0, a3 = −2, a2 = 4, a1 = −3, a0 = 2 ), výstupní hodnoty jsou kořeny (v komplexních číslech 5 kořenů). b) Spočítat neurčitý integrál je úloha matematická, ale ne numerická, neboť výstupem je nekonečně mnoho funkcí (lišících se o konstantu). Příklad 3: Změříme-li hmotnost a průměr koule, bude hustota m 3m 6m ρ= = = 3. 3 V 4πr πd Numerická úloha je vyčíslit hodnotu ρ (výstup) ze vztahu 6m ρ= 3, πd pokud známe vstupní hodnoty m, d a konstantu π . Příklad 4: a) Označme
x = 1.5283
a
y = 12.7391 . Absolutní chyby tedy budou
A( x ) = 0.5 ⋅ 10−4 a A( y ) = 0.5 ⋅ 10−4 .
Výsledek je x + y = 14.2674 , absolutní chyba A( x + y ) = 10−4 a relativní R( x + y ) ≈ 0.7 ⋅ 10−5 . Výsledek bude ležet v intervalu chyba [14.2673, 14.2675]. Analogicky obdržíme výsledky ostatních operací: b) x − y = −0.001 , absolutní chyba A( x − y ) = 10−3 . Při odečítání téměř stejných čísel dochází k tzv. ztrátě platných cifer a velkému nárůstu relativní chyby R( x − y) = 1 . c) x * y = 32.09679 , absolutní chyba A(x − y ) ≈ 0.00799 , relativní chyba R( x * y ) ≈ 0.00025 . d) x / y = 4215.76 , absolutní chyba A( x / y ) = 84.2952 , relativní chyba R( x * y ) ≈ 0.02 . Dochází k velkému nárůstu absolutní chyby, nevhodné v numerických výpočtech.
65
66
Příklad 5: Absolutní chybu spočteme jako chybu funkční hodnoty. Víme, že 6m ρ= 3, πd tedy ⎛ 6m ⎞ ⎛ 6m ⎞ ⎛ 6m ⎞ ∂⎜ 3 ⎟ ∂⎜ 3 ⎟ ∂⎜ 3 ⎟ ´ A( ρ ) = ⎝ πd ⎠ A(m ) + ⎝ πd ⎠ A(d ) + ⎝ πd ⎠ A(π ) . ∂π ∂d ∂m
Příslušné parciální derivace budou ⎛ 6m ⎞ ∂⎜ 3 ⎟ ⎝ πd ⎠ = 6 ∂m πd 3 ⎛ 6m ⎞ ∂⎜ 3 ⎟ ⎝ πd ⎠ = 6m (− 3)d − 4 π ∂d m 6 ⎛ ⎞ ∂⎜ 3 ⎟ ⎝ πd ⎠ = 6m (− 1)π − 2 . ∂π d3 Předpokládejme, že za π dosazujeme číslo π = 3.14 , tedy absolutní chyba bude A(π ) = 0.5 ⋅ 10−2 . Po dosazení hodnot (nezapomeňme převést průměr z centimetrů na metry) dostaneme, že hustota materiálu je ρ = 8935.03 kg ⋅ m −3
A( ρ ) = 286.75 kg ⋅ m −3 R( ρ ) = 0.03 = 3% . Z tabulek zjistíme, že hustota odpovídá hustotě mědi, tedy koule byla patrně vyrobena z mědi. Relativní chyba nepřímého měření je 3% , což se v praxi považuje za přijatelný výsledek.
Příklad 6: Zvolme hodnotu x = 1.57 a A( x ) = 0.01 . Pro takto zvolené hodnoty blízko čísla
π 2
dostaneme, že
R( y ) ≈ 1971.53 , R(x ) čili úloha je velmi špatně podmíněná. Příklad 7: Po vyčíslení relativní chyby na vstupu a na výstupu obdržíme, že číslo podmíněnosti vzhledem na změnu měření hmotnosti je c p ≈ 1 , tedy úloha je dobře podmíněná. Vzhledem na měření průměru koule dostaneme c p ≈ 2.94 , úloha je také dobře podmíněná. Pokud bychom chtěli zlepšit přesnost je výhodnější zpřesňovat měření průměru, neboť má větší vliv na výsledek.
Řešení úloh
67
Řešení kontrolních otázek. Část 2.4
Příklad 1: Grafickou metodou můžeme stanovit jako počátečný interval pro metodu bisekce interval [0, − 2.5] .
x + sin(2^x) y 5
2.5
0 -5
-2.5
0
2.5
5 x
-2.5
Použijeme-li metodu bisekce, po 4-tém kroku dostaneme interval I 4 = [− 0.625, − 0.46875] . Jako startovací hodnotu pro Newtonovou metodu možno použít jakoukoliv hodnotu z intervalu I 4 , použijeme-li hodnotu x0 = −0.625 , stačí k dosažení požadované přesnosti dva kroky x2 = −0.6094965404 , f ( x2 ) = 7.63 ⋅ 10−12 .
Příklad 2: Z grafu odhadneme polohu kořenů: jeden kořen je číslo 0, první kladný kořen leží v intervalu [2.5, 5] a náš hledaný druhý kladný kořen leží v intervalu I 0 = [5, 7] . Metodou bisekce obdržíme I 4 = [6.25, 6.375] , metodou regula-falsi při stejném počtu kroků I 4 = [6.2828, 7] . Na požadovanou přesnost postačí 3 kroky Newtonovy metody, vezmeme-li za x0 = 6.25 , tedy x3 = 6.283185307 , chyba je f ( x3 ) = 2.57 ⋅ 10−15 . Metodou sečen, pro startovací hodnoty získané metodou regula-falsi obdržíme po 3 krocích x3 = 6.283185307 , s přesností 0.42 ⋅ 10−9 . Příklad 3: Grafickou metodou získáme startovací interval pro metodu bisekce I 0 = [0.1, 1] . Po čtyřech krocích bisekce dostaneme I 4 = [0.55, 0.60625] , metodou x0 = 0.55 obdržíme x4 = 0.5671432904 s chybou Newtonovou s f ( x4 ) = 5.42 ⋅ 10−20 . Metodou prosté iterace pro vyjádření x = e − x a x0 = 0.55 je
x5 = 0.5681550202 .
68 Řešení kontrolních otázek. Část 3.4
Příklad 1: Řešení soustavy je (2.917888465, − 0.3882588605, − 6.384972676, 0.07915259566)T . Příklad 2: Řešení soustavy je (0.2160, − 1.1785, − 0.6776 )T . LU-rozklad matice je
1 0 0⎞ 0.002 ⎞ ⎛ ⎛ 5 − 0.12 ⎟ ⎟ ⎜ ⎜ L=⎜ 3 1 0 ⎟ , U = ⎜ 0 7.36 − 2.006 ⎟ . ⎜ − 0.00004 1.4945645 1 ⎟ ⎜0 0 2.881397 ⎟⎠ ⎠ ⎝ ⎝
Příklad 3: ⎛ 0.19977 0.00008 0.00213 ⎞ ⎟ ⎜ Inverze matice A = ⎜ 0.01651 − 0.00550 0.09459 ⎟ . ⎜ 1.55610 − 0.51869 0.34705 ⎟ ⎠ ⎝ −1
Příklad 5: r Jacobiovou metodou pro počátečnou aproximaci x0 = (0,0,0,0)T je řešení po 5ti krocích r x5 =(0.661179, − 0.855741, − 0.446439, 2.666860)T . r Gauss-Seidelovou metodou pro x0 = (0,0,0,0)T je r x5 =(0.681117, − 0.880017, − 0.462579, 2.690817 )T . r Skutečné řešení soustavy je x =(0.682682, − 0.881473, − 0.463193, 2.691742)T . Gauss-Seidelova metoda konverguje rychleji. Příklad 6: r Jacobiovou metodou pro počátečnou aproximaci x0 = (0,0,0,0)T je řešení po 5ti krocích r x5 =(1.288571, − 3.318513, − 3.800583, − 0.170962)T . r Gauss-Seidelovou metodou pro x0 = (0,0,0,0)T je r x5 =(1.182009, − 3.313995, − 3.840567, − 0.198531)T . Skutečné řešení soustavy r je x =(1.188329, − 3.323607, − 3.848806, − 0.201591)T . Gauss-Seidelova metoda konverguje rychleji.
Řešení úloh
69
Řešení kontrolních otázek. Část 4.4
Příklad 1: r Pro jednoduchost zvolíme počáteční vektor y (0 ) = ( 1,1,1)T . Jednotlivé složky vektorů posloupnosti, aproximace vlastního čísla a jejich aritmetický průměr jsou uvedeny v tabulce. k
y1(k )
y2(k )
y3(k )
(λ1 ) 1
(λ1 ) 2
(λ1 ) 3
λ1
0 1 2 3 4 5
1 5 35 265 2015 15285
1 7 53 403 3057 23167
1 7 51 383 2899 21967
5.000000 7.000000 7.571429 7.603774 7.585608
7.000000 7.571429 7.603774 7.585608 7.578345
7.000000 7.285714 7.509804 7.569191 7.577440
6.333333 7.285714 7.561669 7.586191 7.580464
Za aproximaci dominantního vlastního čísla můžeme považovat λ1 = 7.580464 . Další dvě vlastní čísla jsou komplexní 1.2116 − 1.0824i, 1.2116 + 1.0824i . Příklad 2: Vlastní čísla jsou {1.6878, 9.6142, 11.2280} . Mocninnou metodou po 5ti krocích r s počátečnou aproximací y0 = (1,1,1)T je λ1 ≈ 11.6234 . Metodou LU-rozkladu po 10ti krocích dostaneme aproximace vlastních čísel {1.6879, 9.7771, 11.0650}, největší prvek pod diagonálou je a21 = 0.07446 . Příklad 3: Vlastní čísla jsou {4,4,10} , mocninnou metodou i metodou Rayleighova podílu vychází aproximace dominantního vlastního čísla pro počátečnou aproximací r y0 = (1,1,1)T zcela přesně. Metodou LU rozkladu s použitím Choleskeho algoritmu dostaneme po 10ti krocích aproximace vlastních čísel {4.0009, 4.0003, 9.9987} . Největší nediagonální prvek je a21 = 0.07523 . Příklad 4: ⎧ 13 129 13 129 ⎫ , + Vlastní čísla jsou ⎨5, − ⎬ ≈ {5.0, 0.8211, 12.1789}. 2 2 2 ⎭ ⎩ 2
Závěr
ZÁVĚR Distanční opora, kterou jste právě dočetli, obsahuje nezbytný základ pro splnění požadavků k úspěšnému absolvování kurzu Numerická matematika. Pokud jste pracovali samostatně, zpracovali algoritmy všech popisovaných metod ve Vámi zvoleném softwaru a zvládli zodpovědět kontrolní otázky a korespondenční úlohy, dá se předpokládat, že jste učivo zvládli a můžete se přihlásit ke kontrolnímu testu. Pokud Vám některé pasáže dělaly problémy, pokuste se je prostudovat znovu a znovu samostatně vypracujte kontrolní otázky. Při jakýchkoliv nejasnostech kontaktujte tutora, který Vám zajisté poradí, jak při opakování postupovat, případně diskutuje s ostatními studenty o dané problematice. Naleznete-li v textu chyby či nejasnosti, autorka Vám bude vděčná, budete-li ji kontaktovat.
71
72
Další zdroje
LITERATURA 1. BURIAN, K.: Algebra 3: Lineární algebra, Ostrava 1973 2. DĚMIDOVIČ, B. P., MARON, J. A.: Základy numerické matematiky. Praha, SNTL, 1966 3. FIEDLER, M.: Speciální matice a jejich použití v numerické matematice, SNTL Praha, 1981 4. HAVRDA, J., MÍKA, S., PŘIKRYL, P.: Numerické metody a matematická statistika, Praha, ČVUT, 1980 5. KOSTRA J., POMP M.: Teorie matic, Ostravská univerzita, PřF, Ostrava, 1997 6. MÍKA S.: Numerické metody algebry, SNTL Praha, 1982 7. POLÁK, J. : Přehled středoškolské matematiky, 8. vyd., Prometheus, Praha, 2003 8. RALSTON, A.: Základy numerické matematiky. Praha, Academia, 1978. 9. RIEČANOVÁ Z. a kolektív: Numerické metódy a matematická štatistika, ALFA, Bratislava, 1987
73
75
POZNÁMKY
