VYSOKÁ ŠKOLA POLYTECHNICKÁ JIHLAVA Katedra elektrotechniky a informatiky Obor Aplikovaná informatika

VYSOKÁ ŠKOLA POLYTECHNICKÁ JIHLAVA Katedra elektrotechniky a informatiky Obor Aplikovaná informatika

Numerické řešení vlastních funkcí d i f e re n c i á l n í h o o p e r á t o r u bakalářská práce

Autor: Luboš Klouda Vedoucí práce: Doc. RNDr. Ivana Pultarová PhD. Jihlava 2014

Abstrakt Práce se věnuje problému nalezení aproximace nejmenšího vlastního čísla a příslušné vlastní funkce eliptického diferenciálního operátoru dvou proměnných na konvexní polygonální oblasti s homogenními Dirichletovými okrajovými podmínkami s důrazem na řešení náročných praktických případů. Použitou numerickou metodou je metoda konečných prvků pro aproximaci úlohy a metoda inverzní iterace pro nalezení nejmenšího vlastního čísla a vlastního vektoru matice. Součástí práce je porozumění teorii, sestavení programu v prostředí Matlab, Scilab nebo v jazyce C s využitím nástrojů pro úsporné uložení dat a snaha o co nejefektivnější algoritmus. Některé získané numerické výsledky budou porovnány se známým přesným řešením. Autor práce se pokusí o vylepšení tohoto algoritmu např. metodami typu multigrid, kde lze očekávat úsporu iterací vzhledem k tomu, že hledaný vlastní vektor odpovídá nejmenšímu vlastnímu číslu, tedy nízké frekvenci.

Klíčová slova Pozitivně-definitní symetrický diferenciální operátor, variační metoda, problém na nejmenší vlastní číslo, nastavení Visual Studia, základní programovací metody pro lineární numeriku.

Abstract The thesis is dealing with the solution of the least eigenvalue problém including proper eigenfunction of an elliptic differential operator in two variables on a convex polygonal domain with homogenous Dirichlet boundary conditions with emphasis on several practical examples. Numerical method used are: the finite element method for approximation of the problem and inverse iteration method for finding the smallest eigenvalue

and

its

corresponding

eigenvector.

The

thesis

contains

also

the understanding of the theory, writing of program in Matlab, Scilab or in C language with use of tools for optimal data storage and an attempt for as much effective algorithm as possible. Several numerical results will be compared with known exact solutions. The author will try to improve the algorithm, e.g. with methods of multigrid type, where iteration can be saved due to the fact that the desired eigenvector corresponds to the smallest eigenvalue – low frequencies.

Key words Positive definite symmetric differential operator, variatonal method, smallest eigenvalue problem, Visual Studio settings and basic programming methods for linear numerics.

Poděkování Na tomto místě bych rád poděkoval svému vedoucímu práce doc RNDr. Ivaně Pultarové, Ph.D. za poskytnutí tématu a možnost vytvářet ho pod jejím vedením.

Obsah 1

Úvod.......................................................................................................................... 9

2

Předmět práce ......................................................................................................... 11 2.1

3

Formulace úlohy............................................................................................... 11

Funkcionální přístup ............................................................................................... 12 3.1

Pomocná tvrzení ............................................................................................... 13

3.2

Variační metody ............................................................................................... 13

3.2.1

Definice (zobecněného) slabého řešení lineárních diferenciálních rovnic 13

3.2.2

Existence řešení lineární úlohy ................................................................. 14

3.2.3 Aplikace variačních principů na existenci slabého řešení rovnice a problému vlastních čísel, Galerkinova aproximace ............................................. 15 3.2.4 4

Galerkinova aproximace ........................................................................... 16

Metoda konečných prvků........................................................................................ 17 4.1.1

Konečné prvky .......................................................................................... 17

4.1.2

Konstrukce prostoru konečných prvků ..................................................... 17

4.1.3 Odhady, odhady chyb a konvergence lineární úlohy a úlohy o vlastních číslech .................................................................................................................. 18 5

Lineární algebra ...................................................................................................... 19 5.1

Soustava lineárních rovnic ............................................................................... 20

5.1.1

Gaussova eliminace .................................................................................. 20

5.1.2

Choleskiho metoda ................................................................................... 20

5.1.3

Iterační metody ......................................................................................... 21

5.2

Problém vlastních čísel v lineární algebře ....................................................... 21

5.2.1

6

Definice a výpočet vlastních čísel ............................................................ 22

5.3

Podmíněnost matice ......................................................................................... 23

5.4

Multigridní metoda........................................................................................... 23

Zvolené postupy ...................................................................................................... 24 6.1

Řešení soustavy rovnic v inverzní iteraci pro problém vlastních čísel ............ 24

6.1.1

Gaussova eleminace .................................................................................. 24

6.1.2

Choleskiho metoda ................................................................................... 25

6.1.3

Metoda sdružených gradientů ................................................................... 25

6.1.4

Použití multigridu ..................................................................................... 25

6.2

Postup řešení úlohy .......................................................................................... 29

6.2.1

Příprava a ověření předpokladů ................................................................ 29

6.2.2

Vlastní provedení metod ........................................................................... 32

6.2.3

Vlastní aplikace metod na některé konkrétní případy .............................. 32

7

6.2.4

Numerické výsledky experimentů a postupy při programování a ladění . 32

6.2.5

Doby běhů pro pevně stanovený počet iterací .......................................... 39

6.2.6

Parametry a výsledky při iterování do dosažení přesnosti 0.00001 .......... 41

Závěr ....................................................................................................................... 46

Seznam použité literatury ............................................................................................... 49 Seznam obrázků .............................................................................................................. 52 Seznam použitých zkratek .............................................................................................. 53 Přílohy............................................................................................................................. 55 1

Obsah přiloženého DVD ......................................................................................... 55

1 Úvod Předkládaná práce se bude zabývat aplikací metody konečných prvků a pokusí se přiblížit variační metody a jejich aplikace na řešení rovnic a hledání vlastních vektorů a vlastních čísel

a to i těch nejmenších, přiblížit numerické experimenty

v programátorském prostředí Visual Studio 2010 a 2013, včetně využití jejich nastavení a voleb. K problematice jsem se dostal při samostudiu a vyhledávání na internetu, kdy různé problémy z techniky a inženýrské praxe v atomech apod. byly různými metodě konečných prvků víceméně podobnými metodami (meshfree, kolokace, variační) řešeny nebo činěny pomocí výsledků těchto výpočtů závěry. Sama tato metoda tvoří rozsáhlou partii matematiky a jenom pro tuto metodu existuje několikero obměn a zobecnění pro různá specifika úloh. Užití variačních metod a slabých a obecnějších formulací při řešení úloh souvisejících s diferenciálními rovnicemi je rozsáhlá oblast matematiky zahrnující lineární algebru, numerickou matematiku, aplikace funkcionální analýzy, teorie Sobolevových prostorů a různých variant a obměn podobných metodě konečných prvků. Existence řešení konkrétní rovnice nebo úlohy na vlastní čísla není samozřejmá, daří se obecně dokázat jen pro určitou třídu úloh, tedy jen pro určité formulace úlohy s určitými diferenciálními operátory důležitými pro vlastnosti objektů figurujících v úloze a definování nových objektů se zajímavými vlastnostmi. Pro naprogramování úlohy bude používáno výhradně Visual Studio 2010 kombinováno s Visual Studio 2013, dále některé knihovny pro rozšíření jeho funkcí. Pro účely práce jsem předstudoval, ale nevyužil možnosti dalších programovacích nástrojů. Pro aplikace metody jsou třeba pole o velikosti rozměru sítě přibližně čtvrté mocniny rozměru čtvercové oblasti, na níž je úloha definována. Takový požadavek na paměť je pro řešení náročných praktických úloh nepoužitelný. Pro vlastní řešení bude použita Galerkinova metoda, vlastnosti báze v obecném prostoru (že z ní lze vygenerovat libovolný vektor prostoru) a iterační metody pro hledání nejmenšího vlastního čísla z lineární algebry. Pro zobrazování výsledků budu používat výstup textového proudu do souboru a tento načítat a zobrazovat programem Scilab jako matice. Budu zde používat poznatky o operačních systémech a programovacích prostředích z přednášek bakalářského studia. 9

Pokusím se porovnat jednotlivé programátorsko–překladově–optimalizační a numerické a zefektivňující přístupy pomocí měření doby běhu programu časem z procesoru (dostupným v include souboru) a rychlosti konvergence pomocí téhož a počtu iterací pro dosažení přesnosti. Multigridní metody a jejich principy budou také součástí této práce, což jsou dnes současné jedny z aktuálních problémů v tomto odvětví, a na toto téma se stále objevuje řada nových článků. Cílem této práce je seznámit se s metodami a obecnými postupy a aplikovat je na základní (zápis zjednodušující podmínka homogenity v okrajových podmínkách) konkrétní úlohu. Zjistit rychlost konvergence a časovou a operační náročnost jednotlivých přístupů pro konkrétní zvolené metody a pokusit se na různých úrovních o efektivní algoritmus (konání pokusů s různými podmínkami/options). Další oblasti, které tato práce pokrývá, jsou ve stručnosti teoretická hlediska a některé vývody z teorie těchto metod (na existenci řešení a způsoby jeho aproximace) důležité pro jejich používání.

10

2 Předmět práce Problém najít vlastní funkci operátoru se objevuje v mnoha oblastech výpočetní matematiky. Tato úloha se objevuje např. v zpracování signálu, úlohách v teorii obsluhy, při řešen stacionárních stavů různých procesů popsaných diferenciálními operátory, jako např. stacionární Schrödingerova rovnice a z ní odvozené úlohy, stacionární stavy chemických reakcí, radiační transportní úloha v jaderném reaktoru. Tyto úlohy mohou být v diskretizované podobě reprezentovány velkými soustavami lineárních rovnic. Pro diskretizovanou úlohu nalezení nejmenšího vlastního čísla a vlastního vektrou diferenciálního operátoru použijeme metodu inverzní iterace a některé její další úpravy. V praktických výpočtech budeme uvažovat pozitivně definitní symetrický diferenciální operátor druhého řádu na dvourozměrné oblasti. Nejzajímavější bude část o aplikaci některých metod, např. multigridu.

2.1 Formulace úlohy Naším úkolem je nalézt nejmenší reálné číslo λ takové, že existuje spojitá funkce 𝑣 tak, že na oblasti Ω, která je částí 𝑅 2 , platí 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 (𝑎1 (𝑥, 𝑦) ) + (𝑎2 (𝑥, 𝑦) ) + (𝑎2 (𝑥, 𝑦) ) + (𝑎3 (𝑥, 𝑦) ) 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦 = 𝜆𝑣 (1a). a 𝑣(𝑥, 𝑦) = 0 na hranici oblasti Ω, kde funkce a1, a2 a a3 jsou spojité na Ω a matice 𝑎̃(𝑥, 𝑦) = (

𝑎1 (𝑥, 𝑦) 𝑎2 (𝑥, 𝑦)

𝑎2 (𝑥, 𝑦) ) 𝑎3 (𝑥, 𝑦)

je pozitivně definitní pro všechna (𝑥, 𝑦) ∈ Ω a existuje kladná konstanta 𝑐0 tak, že vlastní čísla matice 𝑎̃(𝑥, 𝑦) jsou větší než 𝑐0 pro všechna (𝑥, 𝑦) ∈ Ω. Budeme uvažovat pouze tzv. variační formulaci této úlohy, tedy nalézt funkci v ve vhodném vektorovém prostoru V tak, že pro všechny funkce 𝑢 ∈ 𝑉 platí (1b)

𝑎(𝑢, 𝑣) = 𝜆(𝑢, 𝑣), kde: 𝜕𝑣 𝜕𝑢

𝜕𝑣 𝜕𝑢

𝜕𝑣 𝜕𝑢

𝜕𝑣 𝜕𝑢

𝑎(𝑢, 𝑣) ≔ ∫Ω 𝑎1 (𝑥, 𝑦) 𝜕𝑥 𝜕𝑥 + 𝑎2 (𝑥, 𝑦) 𝜕𝑦 𝜕𝑥 + 𝑎2 (𝑥, 𝑦) 𝜕𝑥 𝜕𝑦 + 𝑎3 (𝑥, 𝑦) 𝜕𝑦 𝜕𝑦 𝑑𝑥𝑑𝑦, (𝑢, 𝑣) = 𝜆 ∫Ω 𝑣𝑢 𝑑𝑥𝑑𝑦 (1c) Hledané číslo λ se nazývá vlastní číslo úlohy a příslušná funkce v se nazývá vlastní funkce této úlohy (funkce 𝑎.. lze indexovat i indexy parciálních derivací). 11

Budeme hledat přibližné řešení, tedy přibližnou hodnotu nejmenšího λ a přibližnou funkci v splňující uvažovanou rovnici. Zvolíme vhodný konečně dimenzionální podprostor 𝑉𝑛 prostoru V a budeme hledat 𝜆𝑛 a 𝑣𝑛 ∈ 𝑉𝑛 splňující uvedenou rovnici pro všechna 𝑢𝑛 ∈ 𝑉𝑛 . Aplikací tohoto postupu dostaneme úlohu najít nejmenší vlastní číslo a příslušnou vlastní funkci symetrické a pozitivně definitní matice. Tuto úlohu lineární algebry ve tvaru Ax = λQx budeme řešit opět přibližně numerickými iteračními metodami. Zvolíme metodu inverzní iterace a použijeme několik jejích variant. Nejzajímavější variantou bude multigridní přístup, tedy úprava poskytující urychlení konvergence pomocí využití vlastností púvodní úlohy. Speciálně bude využita skutečnost, že hledaná vlastní funkce odpovídá nízké „frekvenci“, a tedy lze očekávat její dobrou aproximaci již na velmi hrubé úrovni. Poznamenejme, že použitá multigridní metoda je oproti tradičním multigridním metodám pro úlohu o vlastních číslech zcela novou, v literatuře dosud nepublikovanou metodou.

3 Funkcionální přístup Definice (vlastní číslo): Řekneme, že 𝜆 ∈ 𝐶 je vlastní číslo zobrazení 𝐴: 𝑉 ⟶ 𝑉, kde 𝑉 je lineární vektorový prostor, právě když existuje 𝑣 ∈ 𝑉\{0} a platí 𝐴𝑣 = 𝜆𝑣, vektor 𝑣 nazveme potom vlastním vektorem. Skalární součin a normu lze definovat několika způsoby. Pro 𝑘 ≥ 0, 1 ≤ 𝑝 < ∞ 1/𝑝

𝛼 definujeme normu ‖𝑓‖𝑊𝑝𝑘(Ω) : = (∑|𝛼|≤𝑘‖𝐷weak 𝑓‖𝑝𝐿𝑝 (Ω) )

, prostor 𝑊𝑝𝑘 (Ω) jako prostor

funkcí, které mají příslušnou normu konečnou. Pro 𝑝 = ∞ potom ‖𝑓‖𝑊∞𝑘(Ω) ≔ 𝛼 max|𝛼|≤𝑘 ‖𝐷weak 𝑓‖𝑝𝐿∞ (Ω) a znovu pro funkce mající tuto hodnotu konečnou,

pro něž dává smysl ji definovat, prostor 𝑊∞ 𝑘 (Ω). Zde 𝛼 jsou multiindexy a to včetně opakujících se variant a mají symbolizovat indexy parciálních derivací (viz [1]). Tyto prostory splňují definici vektorového normovaného prostoru. Pro 𝑝 = 2 lze analogicky 𝛼 𝛼 definovat také skalární součin (𝑢, 𝑣)𝑊2𝑘 (Ω) = ∑|𝛼|≤𝑘 ∫ 𝐷weak 𝑢. 𝐷weak 𝑣 𝑑Ω.

Vektorový prostor se skalárním součinem, který je úplný (každá cauchyovská posloupnost je v něm konvergentní) se nazývá Hilbertův prostor [2]. Banachovým prostorem nazveme úplný normovaný prostor.

12

3.1 Pomocná tvrzení Věta: Nechť 𝑇 je totálně spojitý (ohraničenou množinu zobrazuje na množinu, kde lze z každé posloupnosti vybrat podposloupnost konvergentní v nadprostoru) operátor v separabilním ℋ [2], 𝜀 libovolné kladné číslo. Pak lze najít takové operátory 𝑇1 a 𝑇2 , že platí 𝑇 = 𝑇1 + 𝑇2 , přičemž operátor 𝑇1 má konečnou dimenzi a pro normu operátoru platí ‖𝑇2 ‖ < 𝜀. Více např. v [2]. Pro některé závěry je i důležitá věta o pevném bodu (apod. viz např.[3]). (1)

Věta: Identita z 𝑊2 (Ω) do 𝐿2 (Ω) je totálně spojitá. Důkaz lze najít např. v [4].

3.2 Variační metody Vlastnosti úlohy na vlastní čísla souvisejí s vlastnostmi řešení variačních formulací lineárních diferenciálních rovnic.

3.2.1 Definice (zobecněného) slabého řešení lineárních diferenciálních rovnic Definice: Mějme oblast 𝐺

s Lispchitzovskou [2] hranicí Γ, operátor

𝐴=

∑𝑖,𝑗≤𝑘(−1)|𝑖| 𝐷𝑖 (𝑎𝑖𝑗 𝐷 𝑗 ) s omezenými měřitelnými koeficienty a bilineární formu, (𝑘) pro niž platí |(𝑣, 𝑢)| ≤ 𝐾‖𝑣‖𝑊 (𝑘) (𝐺) ‖𝑢‖𝑊 (𝑘) (𝐺) , 𝑢, 𝑣∈𝑊2 (𝐺), 𝐾 > 0, funkce 𝑓 ∈ 2

𝐵11 𝐿2 (𝐺), operátory … 𝐵𝑟1

2

… 𝐵1𝜇 … … na jednotlivých částech Γ1 , … , Γ𝑟 hranice Γ; prostor … 𝐵𝑟𝜇𝑟

(𝑘) 𝑉 = {𝑣; 𝑣 ∈ 𝑊2 (𝐺), 𝐵11 (𝑣) = 0, … , 𝐵𝑟𝜇𝑟 (𝑣) = 0 na Γ };

funkce

𝑔𝑝1 ∈

(𝑘)

𝐿2 (Γ𝑝 ), … , 𝑔𝑝𝜇𝑝 ∈ 𝐿2 (Γ𝑝 ), 𝑝 ∈ 1, … , 𝑟, a funkce 𝑤∈𝑊2 (𝐺), pro niž platí na Γ𝑝 (𝑝 = 1, … , 𝑟) 𝐵𝑝1 𝑤 = 𝑔𝑝1 (𝑆), … , 𝐵𝑝𝜇𝑝 𝑤 = 𝑔𝑝𝜇𝑝 (𝑆), funkce ℎ𝑝𝑙 (𝑆) ∈ 𝐿2 (Γ𝑝 ), 𝑝 = 1, … , 𝑟, (𝑘)

𝑙 = 1, … , 𝑘 − 𝜇𝑝 . Funkci 𝑢 ∈ 𝑊2

nazveme slabým řešením problému s okrajovými

podmínkami daného v této větě zavedenými symboly, jestliže 𝑢 − 𝑤 ∈ 𝑉, 𝑡 𝜕 𝑝𝑙 𝑣

𝑘−𝜇

((𝑢, 𝑣)) = (𝑣, 𝑓) + ∑𝑟𝑝=1 ∑𝑙=1 𝑝 ∫Γ

𝑝

𝑡 𝜕𝑣 𝑝𝑙

ℎ𝑝𝑙 𝑑𝑆

platí

pro

každé

𝑣 ∈ 𝑉.

13

Poznámka: pro specifické třídy operátorů lze dokázat, že zobecněné slabé řešení je i tzv. řešením regulárním, tj. klasickým. Definice je převzata z [2], více obecněji v [3].

3.2.2 Existence řešení lineární úlohy Věta (Rieszova): Každý lineární omezený [2] funkcionál v Hilbertově prostoru lze vyjádřit ve tvaru 𝜑(𝑢) = (𝑢, 𝑣), kde 𝑣 je prvek tohoto prostoru, který je tímto funkcionálem jednoznačně určený. Nechť je lineární operátor 𝐴 pozitivně definitní na husté [2] podmnožině v Hilbertově prostoru. Zkonstruujeme-li nový prostor se skalárním součinem ((𝑢, 𝑣)) = (𝑢, 𝐴𝑣), pak funkcionál 𝐹(𝑢) = ((𝑢, 𝑢)) − 2(𝑓, 𝑢) nabývá v tomto novém prostoru svého jednoznačného minima (viz [2]). Věta (Laxova-Milgramova, varianta 1): Nechť (𝑋, ‖. ‖𝑋 ) je Banachův prostor, 𝑙 ∈ 𝑋 ′ je spojitý lineární funkcionál a 𝑎: 𝑋 × 𝑋 ⟶ 𝑅 je spojitá bilineární symetrická (tj. 𝑎(𝑣, 𝑢) = 𝑎(𝑢, 𝑣) pro libovolné 𝑢, 𝑣) forma, která je koercivní (V-eliptická), tj. existuje 𝛼 > 0 takové, že 𝑎(𝑢, 𝑢) ≥ 𝛼‖𝑢‖𝑋 1 2

2

∀𝑢 ∈ 𝑋, pak funkcionál 𝐽(𝑢) ≔

𝑎(𝑢, 𝑢) − 𝑙(𝑢) má jedinou minimalizující hodnotu pro 𝑢∗ ∈ 𝑋. Tato minimalizující

hodnota je jediné řešení

𝑎(𝑢∗ , 𝑣) = 𝑙(𝑣) ∀𝑣𝜖𝑋, přitom je ‖𝑢∗ ‖ ≤ ‖𝑙‖/𝛼.

Tvrzení této a následujících tří vět lze nalézt v [5], důkaz obdobného tvrzení můžeme najít např. v [2]. Věta (Laxova-Milgramova, varianta 2): Předpoklady a značení jsou jako v Laxově |𝑎(𝑢,𝑣)| . 𝑋 ‖𝑣‖𝑋

Milgramově větě varianatě 1 výše. Označme 𝐴 ≔ ‖𝑎‖𝐶←𝑋×𝑋 ≔ 𝑠𝑢𝑝𝑢,𝑣∈𝑋\{0} ‖𝑢‖

Nechť 𝑆 ⊆ 𝑋 je konečně dimenzionální podprostor 𝑋. Jednoznačná minimalizující funkce označená 𝑢𝜖𝑋, resp. 𝑢𝑆 𝜖𝑆 funkcionálu 𝐽 v 𝑋, resp. v 𝑆, splňuje ‖𝑢 − 𝑢𝑆 ‖𝑋 ≤ 𝐴 𝛼

𝑖𝑛𝑓𝑣∈𝑆 ‖𝑢 − 𝑣‖𝑋 . Nechť dodatečně je 𝐻 Hilbertův prostor se skalárním součinem

(. , . )𝐻 a normou ‖. ‖𝑋 tak, že 𝑋 je spojitě a hustě vnořen do 𝐻. Pro 𝜑𝜖𝐻 označme 𝑢𝜑 ∈ 𝑋

jediné

řešení

𝑎(𝑣, 𝑢𝜑 ) = (𝜑, 𝑣)𝐻

pro

∀𝑣 ∈ 𝑋.

Pak

‖𝑢 − 𝑢𝑆 ‖𝐻 ≤ A‖𝑢 −

𝑢𝑆 ‖𝑋 𝑠𝑢𝑝𝜑∈𝐻\{0} 𝑖𝑛𝑓𝑣∈𝑆 ‖𝑢𝜑 − 𝑣‖ / ‖𝜑‖𝐻 . 𝑋

Věta (Laxova-Milgramova, varianta 3): Nechť je (𝑋, ‖. ‖𝑋 ) Banachův prostor a 𝑎0 , 𝑎1 : 𝑋 × 𝑋 ⟶ 𝐶 dvě spojité (viz 2.5.2 v [5]) bilineární formy. O 𝑎0 předpokládáme, že je hermitovská a koercivní, dále navíc u 𝑎1 požadujeme (i to co pro 𝑎0 ) koercivitu 14

a existenci konstanty 𝐴 ∈ 𝑅 + \{0}, že |𝑎1 (𝑢, 𝑣)| ≤ 𝐴‖𝑢‖𝑋 ‖𝑣‖𝑋 ∀𝑢, 𝑣 ∈ 𝑋. Dále položme 𝑎 ≔ 𝑎0+𝑎1 jako rovnost pro každé dva body. Pro všechny 𝑢 ∈ 𝑋\{0} předpokládáme 𝑎(𝑢, 𝑢) ≠ 0. Pak mají 𝑎(𝑢, 𝑣) = 𝑙(𝑣) pro ∀𝑥 ∈ 𝑋 a 𝑎(𝑣, 𝑢) = 𝑙(𝑣) pro ∀𝑥 ∈ 𝑋 jednoznačné řešení pro libovolný spojitý lineární funkcionál 𝑙 ∈ 𝑋 ′ . Taktéž viz [3]. Věta (Laxova-Milgramova, varianta 4): Pro Hilbertův prostor (ℋ, (. , . )) a spojitou, bilineární formu 𝑎(. , . ), pro kterou existuje konstanta 𝜉 ∈ 𝑅 + \{0}, že |𝑎1 (𝑢, 𝑣)| ≥ 𝜉‖𝑢‖𝑋 ‖𝑣‖𝑋 ∀𝑢, 𝑣 ∈ 𝑋, a spojitý lineární funkcionál 𝐹 ∈ ℋ ′ , existuje jediné 𝑢 ∈ ℋ takové, že 𝑎(𝑢, 𝑣) = 𝐹(𝑣) ∀𝑣 ∈ ℋ. Tuto větu jsem převzal z (viz 2.5.2) [5], [4]. Věta: Převezměme značení z varianty 3 Laxovy-Milgramovy věty. Dále nechť 𝑆 ⊆ 𝑋 je konečně dimenzionální podprostor. Pak 𝑎(𝑢𝑆 , 𝑣) = 𝑙(𝑣) pro ∀𝑣 ∈ 𝑆 má jednoznačné řešení 𝑢𝑆 ∈ 𝑆 pro libovolné 𝑙 ∈ 𝑋 ′ . Nechť je 𝑎 navíc koercivní, tj. existuje 𝛽 > 0 2

s 𝑎(𝑢, 𝑢) ≥ 𝛽‖𝑢‖𝑋 pro ∀ 𝑢 ∈ 𝑋. Pak jednoznačná řešení 𝑢 a 𝑢𝑆 splňují odhad chyby 𝐴

‖𝑢 − 𝑢𝑆 ‖𝑋 ≤ 𝑖𝑛𝑓𝑣∈𝑆 ‖𝑢 − 𝑣‖𝑋 , kde 𝐴 ≔ ‖𝑎‖𝐶←𝑋×𝑋 . Navíc nechť jsou 𝐻, 𝜑, 𝑢𝜑 jsou 𝛼 jako v předešlém. Nechť dodatečně je 𝐻 Hilbertův prostor se skalárním součinem (. , . )𝐻 a normou ‖. ‖𝑋 , tak jako ve variantě 3 Laxovy-Milgramovy věty, v její části o hustém vnoření 𝑋 v 𝐻. Pro 𝜑𝜖𝐻 označme 𝑢𝜑 ∈ 𝑋 jediné řešení 𝑎(𝑣, 𝑢𝜑 ) = (𝜑, 𝑣)𝐻 pro ∀𝑣 ∈ 𝑋. Pak ‖𝑢 − 𝑢𝑆 ‖𝐻 ≤ A‖𝑢 − 𝑢𝑆 ‖𝑋 𝑠𝑢𝑝𝜑∈𝐻\{0} 𝑖𝑛𝑓𝑣∈𝑆 ‖𝑢𝜑 − 𝑣‖ / ‖𝜑‖𝐻 . 𝑋

Toto tvrzení lze nalézt v [3].

3.2.3 Aplikace variačních principů na existenci slabého řešení rovnice a problému vlastních čísel, Galerkinova aproximace Uvažujeme úlohu na vlastní čísla (1b). Funkcionálem z předchozí kapitoly Laxovy– Milgramovy věty varianty 4 bude pro následující odstavec (𝑓, . ). Řekneme, že operátor 𝐾 z normovaného prostoru do normovaného prostoru je kompaktní, jestliže v obrazu omezené množiny lze z každé poslouponosti vybrat konvergentní posloupnost. To je ekvivalentní tvrzení pro

Hilbertovy separabilní prostory, že obraz každé slabě

konvergentní posloupnosti konverguje v normě- tj. silně. Z kapitoly o variačním počtu a z [2], [3] víme, že operátor 𝐾: 𝑓 ⟶ 𝑢 je kompaktní a spojitý. Existence řešení úlohy na vlastní čísla pro symetrickou a V - eliptickou bilineární formu vyplývá 15

z Riesz-Schauderovy věty a Fredholmovy alternativy (viz [2] a [5]). Pro naši úlohu z pozitivní definitnosti 𝑎 výplývá, že vlastní čísla mohou být pouze kladná a jediným hromadným bode je nekonečno. V literatuře se dovozuje, že když je forma 𝑎 hermitovská je 𝐾 samoadjungovaný (diskuze v [5]). Dále viz [6], podobně existence a nabývání extrému v [7]. Dále diskuze v [8], [5]. (1)

Při ověřování vlastnosti kompaktnosti 𝐾 je důležité tvrzení, že identita z 𝑊2 (Ω) do 𝐿2 (Ω) je totálně spojitá (viz [4]) (dále viz např.[2], [3], s návodem k dalšímu [7]).

3.2.4 Galerkinova aproximace Z tohoto variačního přístupu nebo pomocí tzv. Galerkinovy metody lze pro konkrétní bázi prostoru funkcí odvodit soustavu rovnic Ax = λQx, která poskytuje aproximaci přesného

řešení.

Vhodnou

volbou

posloupností

těchto

prostorů

(už i konečnědimenzionálních) dostaneme posloupnost hodnot a funkcí, které tedy konvergují k přesnému řešení (viz odkazy z předchozí kapitoly). Pro účel této kapitoly uvažujme dále jen bilineární formy. (𝑛)

Věta: Nechť je forma 𝑎 spojitá symetrická koercitivní (V-eliptická) a nechť 𝜆1 ≤ (𝑛)

(𝑛)

𝜆2 ≤ … ≤ 𝜆𝑛

jsou

aproximace

vlastních

čísel

získané

řešením

konečně

dimenzionální variační úlohy na úplný systém vlastních čísel. Pak pro každé vlastní (𝑛)

číslo 𝑎 platí 𝜆𝑖 = lim 𝜆𝑖 𝑛→∞

≥ 0 (viz varianta pro jednonásobná vl.č. v [2], nebo [3]).

Dále teorie viz [9]. Co se týče konvergence vlastních hodnot, tu lze vyvodit z toho, že když rezolventa přesného problému obepíná 𝑚 hodnot, bude je v limitě a tedy i od nějakého indexu obepínat i rezolventa aproximace vlastních čísel na konečně dimenzionálním prostoru Galerkinovy aproximace (v této teorii je v jednom úseku potřeba stejnoměrná omezenost

rezolvent, kterou lze zajistit,

či

izolovanost

vlastního

čísla).

Pro jednonásobná vlastní čísla nebo vlastní čísla přesného operátoru, která mají nenulovou vzdálenost od zbytku spektra aproximace, lze najít odhad aproximace vlastního čísla a tedy i konvergenci metody (viz [5]).

16

4 Metoda konečných prvků 4.1.1 Konečné prvky Řešení variační úlohy není prakticky dosažitelné. Jako jedna z variant prostoru pro Galerkinovu konečně-dimenzionální aproximaci řešení variační úlohy lze volit funkce definované pomocí konečných prvků.

4.1.2 Konstrukce prostoru konečných prvků Definice: Nechť 𝐾 ⊆ 𝑅 𝑛 je oblast s po částech spojitou hranicí (oblast prvku), 𝑃 je konečně dimenzionální prostor funkcí v 𝐾 (tvarové funkce), 𝒩 = {𝑁1 , 𝑁2 , … , 𝑁𝑘 } je báze pro duální prostor 𝑃′ (uzlové funkce nebo proměnné). Pak trojici (𝐾, 𝑃, 𝑁) nazveme konečný prvek. V doporučené literatuře je řada příkladů takových prostorů (viz [4], [5]). Definice (lokální interpolant): Mějme konečný prvek (𝐾, 𝑃, 𝒩), nechť dále množina {𝜙𝑖 : 1 ≤ 𝑖 ≤ 𝑑} je báze duální k 𝒩. Je-li pro 𝑣 funkce, pro kterou jsou všechny 𝑁𝑖 ∈ 𝒩, 𝑖 = 1, … , 𝑑 definovány, pak definujeme lokální interpolant pomocí: ℐ𝐾 𝑣: = ∑𝑑𝑖=1 𝑁𝑖 (𝑣)𝜙𝑖 . Tento lok. interpolant má vlastnosti: je lineární na 𝑃, na 𝑃 se chová jako identita. Definice (zjemnění neboli podrozdělení): Podrozdělením oblasti Ω je konečný počet ̅𝑖 = Ω ̅. otevřených množin {𝐾𝑖 }, takové, že 𝐾𝑖 ∩ 𝐾𝑗 = ∅ pro 𝑖 ≠ 𝑗 a ∪ 𝐾 Poznámka -definice (interpolant): obsahují- li uzlové proměnné derivace do řádu 𝑚, pak pro funkce diferencovatelné do řádu 𝑚 definujeme globální interpolant po jednotlivých

oblastech

jako

lokální

interpolanty.

Poznámka: Existuje teoretický koncept s využitím obecného integrálu, že lze interpolant definovat i pro nespojité funkce. Konečné prvky jsou několikero typů trojúhelníkových, které se liší 𝒩. Tvrzení: Vezměme vhodný prostor konečných prvků tak, že je možné vzít vrcholy tak, že obraz funkce globálním interpolantem je vždy podle typu prvku několikrát (𝑟) diferencovatelný, navíc je dokonce z prostoru 𝑊∞𝑟+1 . Náznak důkazu: dokážeme spojitost interpolantu pro každou hranu a zbytek se dále 17

ověří už podle definice. Učiňme pro tuto chvíli předpoklad, že body prvku (mj. tam kde beru funkcionály 𝑁) jsou pro oba hraničící lokální interpolanty stejné (mimo svoji oblast je lze dodefinovat jako polynomy a to proto, že je to konečná suma), a definujme rozdíl obou interpolantů na hranici jako funkci 𝑤. A protože 𝑤 je polynom konečného stupně a jeho zúžení na hranu má uzlové proměnné nulové, tudíž musí být sám nulový. Tedy je interpolant spojitý podél hrany. Ověření spojitosti derivace se učiní využitím per partes metody. Lze definovat několik relací ekvivalence mezi prvky: –afinní a –interpolační. Nebo kombinace obou. Definice: Řekneme, že omezená oblast Ω je hvězdicová vzhledem ke kouli 𝐵 ⊆ Ω, právě když (∀𝑥 ∈ Ω)(konvexní obal(𝑥, 𝐵) ⊆ Ω). Tato vlastnost je důležitá při konstrukci podoby zbytku zobecněného Taylorova polynomu a odhadů (závisí na robustnosti a ta na poloměru maximální vepsané kružnice danému trojúhelníku příslušné triangulace) jeho zbytku (Bramble- Hilbertova lemmatu), které závisí na diam(Ω). Tvrzení (Sobolevova nerovnost): Nechť oblast Ω ⊆ 𝑅 𝑛 , diam(Ω) = 𝑑, hvězdicová vůči nějaké kouli 𝐵, 𝑢 ∈ 𝑊𝑝𝑚 (Ω), 𝐶0∞ (Ω), pak 𝑢 ∈ 𝐶(Ω) a platí ‖𝑢‖𝐿∞(Ω) ≤ 𝐶𝑚,𝑛,𝛾,𝑑 ‖𝑢‖𝑊𝑝𝑚 (viz [4], a v [5] kde se definuje v 4.2.16 i parametr 𝛾). Odhady odchylky interpolantu jsou funkcí diam(𝐾), na kterém jsou činěny, tzv. lokální chyba interpolace.

4.1.3 Odhady, odhady chyb a konvergence lineární úlohy a úlohy o vlastních číslech Definice (lokální indikátor chyby): Pro Galerkinovské 𝑢ℎ definujme ℰ𝑒 (𝑢ℎ ) ≔ ∑𝑇⊆ 𝑇𝑒 ℎ2𝑇 ‖𝑓 + ∇. (𝛼∇𝑢ℎ )‖2 𝐿2 (𝑇) + ℎ𝑒 ‖[𝛼𝑛. ∇𝑢ℎ ]𝑛 ‖2𝐿2 (𝑒) , kde oboje ℎ značí míru 𝐾 umocněnou na 1/dim(𝐾), 𝑒 je hrana a 𝑛 je k ní normálový vektor a 𝛼 je definováno vztahem 9.0.2 v [5] - tj. udává tvar operátoru. Tento vztah je podstatný pro odhad chyby interpolace Galerkinovského řešení od přesného shora pomocí této veličiny a minima kvadratické formy.

18

Podobně můžeme obdržet odhad chyby pro řešení zobecněného řešení zdola pro dvě oblasti sdílející hranu pomocí maxima kvadratické formy. Dále lze odvodit ještě odhad pro vnitřek funkce v lokálním indikátoru chyby, ale pod maximem a ne pod integrálem. V literatuře je pro adaptivně zjemňovanou mřížku odvozen odhad pro lokální indikátor 2

chyby pro zjemňovanou mřížku: ℰ𝑒 (𝑢𝑗 ) ≤ 𝐶 ∑𝑇∈𝑇𝑒 (‖|𝑢𝑗+1 − 𝑢𝑗 |‖ 𝑇 + 𝑜𝑠𝑐𝑗 (𝑇)2 ) , kde osc je definován jako úměrný normě rozdílu integrovaného vnitřku indikátoru chyby a její projekce na prostor počástech-polynomů vzhledem k triangulaci, kde se tvrzení dokazuje na trojúhelníku, kterému přibude vnitřní bod. Podobně můžeme najít odhad pro nový osc pomocí starého a normy rozdílu dvou výsledků pro dvě různá zjemnění. Navíc máme odhad i pro rozdíl od přesného řešení v každé iteraci. V konečněprvkové bázi lze pro bilineární formu najít odhad pomocí složek bazických vektorů a z nich znovu pomoci integrálních odhadů a Hölderovy nerovnosti zeshora i největší vlastní číslo normou 𝑎 a nejmenšího vlastního čísla jako převrácenou hodnotu normy𝑎 a součinu s 2/𝑛 mocninou počtu trojúhelníku v triangulaci. Podmíněnost matice A lze odhadnout při splnění 9.6.2 z [5] i jako 𝜘(𝑎) ≤ 𝐶. 𝑁(1 + 2

|𝑙𝑜𝑔(𝑁 ℎ𝑚𝑖𝑛 (𝑁)2 )|), tj. pro nejmenší vlastní číslo 𝜆𝑚𝑖𝑛 ≥ (𝑁ℎ𝑚𝑖𝑛 (𝑁)2 )𝑝 /(𝐶𝑝𝑁) pro 𝑝 > 2, pro podmíněnost pak už jenom stačí zvolit vhodné 𝑝 a získáme odhad 1

(kde ‖𝑎‖ určuje podobným způsobem 𝐶 a pří odhadu se používá koercivita 𝑎). Konvergence v aproximačním prostoru lze zvolit tak, že |𝜆 − 𝜆(𝑗) | < 𝐷. ℎ𝑟 , kde 𝑟 je stupeň polynomu v konečných prvcích a 𝐷 je vhodná konstanta. Tato tvrzení jsou podložena důkazy v [3], [6], [2], dále též [9], více teoreticky v [10].

5 Lineární algebra V následujících kapitolách popíšeme řešení dosud získané konečněrozměrné úlohy o vlastních číslech Ax = λQx, kde 𝐴 a 𝑄 jsou reálné symetrické pozitivně definitní matice.

Nejdříve

uvedeme

metody

pro

řešení

soustav

lineárních

rovnic

a potom ukážeme jejich využití při řešení úlohy na vlastní čísla.

19

5.1 Soustava lineárních rovnic V praxi se používají přímé metody (Gaussova eliminace, Choleskiho rozklad) a iterační metody (např. metoda sdružených gradientů).

5.1.1 Gaussova eliminace Gaussova eleminace je univerzální metoda řešení soustavy lineárních rovnic. Z hlediska implementace na počítači ji můžeme doplnit částečnou nebo úplnou pivotaci. Tato známá metoda převede řešení soustavy na soustavu s horní trojúhelníkovou maticí se stejnou množinou řešení. Pro každý sloupec provedeme následující úpravu: první řádek zbytku matice pro další výpočty podělíme prvkem jeho příslušného sloupce. Pak aktuální řádek násobíme pro každý jednotlivý řádek mínus hodnotou jednotlivého řádku v příslušném sloupci a jako vektor ho k tomuto řádku přičteme. V dalším cyklu bereme další sloupec. Je patrné, že metoda kromě jediného prvku v každém sloupci vyrábí nuly. Tyto využijeme při zpětném chodu hledání řešení rovnice. Existuje varianta Gaussovy eliminace se zpětným chodem, který se realizuje obdobně, jen se postupuje odspoda nahoru, tedy nenulové zůstanou jenom diagonální prvky. Tato metoda je používána pro určení inverzní matice pro menší případy.

5.1.2 Choleskiho metoda Pro pozitivně definitní symetrické matice existuje metoda, která umožňuje matici rozložit na součin dvou trojúhelníkových matic dolní krát horní, tatáž matice. Výhodou této metody je, že na rozdíl od Gaussovy eliminace stačí pro její použití pro řešení

soustavy

s trojúhelníkovou

lineárních

maticí,

ale

rovnic dvakrát.

pouze Tato

algoritmus

řešení

soustavy

metoda

bude

popsána

z numericko-programátorského hlediska popsána v části zabývající se praktickými aplikacemi. Pro zmiňovaný typ matic (pozitivně definitní symetrické, mohla by fungovat i pro symetrické, pokud budou mít podmínky v tělese hodnoty) existuje tento rozklad

20

2

2 jako 𝐴 = 𝐿. 𝐿𝑇 . Matici sestrojíme podle vztahu 𝐿𝑟𝑟 ≔ √𝐴𝑟𝑟 − ∑𝑟−1 𝑠=1 (𝐿𝑟𝑠 ) a pro 𝑘 = 𝑟 + 1

1, … , 𝑛 platí pro další členy vztah 𝐿𝑘𝑟 ≔ 𝐿 (𝐴𝑘𝑟 − ∑𝑟−1 𝑠=1 𝐿𝑟𝑠 𝐿𝑘𝑠 ). 𝑟𝑟

5.1.3 Iterační metody Metoda sružených gradientů Pro pozitivně definitní symetrickou matici lze použít následně popsanou

iterační

metodu sdružených gradientů [11], [12]. Mějme maticovou rovnici 𝐴𝑥 = 𝑏. V každém kroku (v prvním 𝑟0 = 𝐴𝑥 − 𝑏, 𝑝0 = 𝑟0 ). Pro 𝑘 − 1 ⟶ 𝑘 tedy na každém dalším kroku: 𝑇 α𝑘 = (𝑟𝑘−1 𝑝𝑘−1 )/(𝑝𝑘−1 𝑇 𝐴𝑝𝑘−1 ), 𝑥𝑘 = 𝑥𝑘−1 + α𝑘 𝑝𝑘−1 , 𝑟𝑘 = 𝑟𝑘−1 − α𝑘 𝐴𝑝𝑘−1 , 𝛽𝑘 = 𝑇 (𝑟𝑘𝑇 𝑟𝑘 )/(𝑟𝑘−1 𝑟𝑘−1 ) = −(𝑟𝑘𝑇 𝐴𝑝𝑘−1 )/(𝑝𝑘−1 𝑇 𝐴𝑝𝑘−1 ), 𝑝𝑘 = −𝑟𝑘 + 𝛽𝑘 𝑝𝑘−1 .

Vztah pro každé 𝛼𝑘 lze dovodit z při známém směru 𝑝𝑘 z podmínky minimalizovat chybu 𝑥𝑘 od přesného řešení - tedy vzdálenost měřenou v normě - (toto lze učinit i pro obecnou matici 𝐴 ostatní operace už nikoliv). Ve směru daném směrovým vektorem minimalizuje kvadratický funkcionál. Za předpokladu, že jsou 𝑝𝑘 mezi sebou dolů vzájemně 𝐴 ortogonální, se podaří odvodit vztah pro 0 = (𝑝𝑘−1 , 𝐴𝑝𝑘 ) a 𝑝𝑘 = −𝑟𝑘 + 𝛽𝑘 𝑝𝑘−1 ⇒ 𝛽𝑘 . 𝑝𝑘+1 𝑇 𝐴𝑝𝑖 = −𝑟𝑘+1 𝑇 𝐴𝑝𝑖 + 𝛽𝑘 𝑝𝑘 𝑇 𝐴𝑝𝑖 obecně první člen díky postupu a vlastnostem nově vytvářeného prostoru – nově vytvářená 𝑝𝑘 i 𝑑𝑘 nepatří do obalu předchozíchdruhý díky platnosti indukčního předpokladu. Je možno odvodit druhý vztah pro 𝛽𝑘 a to také z platnosti vztahů mezi vytvářenými prostory (je to ten druhý v definici než používáme pro výpočet). Poznámka: existuje i metoda bikonjugovaných gradientů, kde se pracuje s maticí a s ní adjungovanou, které jsou zdrojem vedle sebe jdoucích iterací dvojice vektorů, tato metoda [11] však v pestré paletě článků má řadu úprav [12], [13] a vylepšení, řadu komentářů. Další metody viz [3].

5.2 Problém vlastních čísel v lineární algebře V konečně dimenzionální lineární algebře platí tvrzení, že matice je diagonalizovatelná v bázi vlastních vektorů, právě když se všechny geometrické a algebraické násobnosti pro každé z vlastních čísel shodují. Pro kvadratické formy platí

tvrzení, že počet kladných a záporných a nulových 21

koeficientů v diagonální bázi je stejný a neměnný s bazí, tedy je pro každou tuto bázi stejný.

5.2.1 Definice a výpočet vlastních čísel Tvrzení (pro konečně dimenzionální prostory ekvivalentní definice): Mějme konečně dimenzionální prostor, vlastním číslem nazveme každé řešení rovnice det(𝐴 − 𝜆𝐸) = 0, kde 𝐸 je jednotková matice. Alternativně viz [13]. Řešení problému vlastních čísel

Úplná úloha na vlastní čísla Řešit úlohu na nalezení všech vlastních čísel je pro velké matice náročná úloha, většina metod přitom vychází právě z výše uvedeného tvrzení o nalezení kořenů polynomu obecně takového řádu, jako je dimenze prostoru, na kterém úlohu řešíme. Pro opravdu velké

matice

se

na

toto

používají

iterační

metody.

Další obměnou téhož je úloha najít dokonce aproximace příslušných vlastních vektorů. Kde jsou taktéž vypracovány metody, jednak iterační, jednak můžeme hledat vlastní vektor, známe-li, že se jedná o vlastní číslo. Tyto metody (např. LR, QR, trojúhelníkový rozklad) nemusejí vždy (všechny prvky) konvergovat, ale o vždy platí, že iterované matice mají stejná vlastní čísla jako původní matice, či hodnoty čísel na hlavní diagonále. Samozřejmě i podobně jako v následující kapitole existují iterační metody řešení úplného problému vlastních čísel (viz předchozí kapitola), my se o nich ale zde zmiňovat nebudeme.

Hledání nejmenšího vlastního čísla inverzní iterací Jelikož předmětem práce je najít nejmenší vlastní číslo, lze toto pro diagonalizovatelnou 𝑑𝑖𝑚 matici připodobnit takto: nechť je 𝐴̃(∑𝑑𝑖𝑚 𝑖=1 𝑥𝑖 ) = ∑𝑖=1 𝜆𝑖 𝑥𝑖 a nechť jsou v sumě vlastní

čísla seřazená od nejmenšího (v absolutní hodnotě) po největší, 𝑥𝑖 jsou vlastní vektory, pak lze nejmenší vlastní číslo určit inverzní iterací. Je totiž 𝐴̃−1 (𝐴̃−1 (∑𝑑𝑖𝑚 𝑖=1 𝑥𝑖 𝑎𝑖 )) = 1/𝜆2nejmenší ∑dim 𝑖=1

𝜆2nejmenší 𝜆2𝑖

. 𝑥𝑖 𝑎𝑖 a tady po 𝑘-té iteraci se koeficienty složky vektoru v bázi

u větších vlastních čísel zmenšují tak rychle, jak moc se liší od toho nejmenšího. 22

Metoda inverzní iterace tedy spočívá v opakování kroků: v k+1 = A−1 Qv k , μ¨k+1 = ‖v k ‖/‖v kˇ+1 ‖, v k+1 = v k+1 /‖v k+1 ‖. Tato metoda nalezení nejmenšího vlastního čísla tedy konverguje právě k tomu nejmenšímu vlastnímu číslu. Vektory zaručeně konvergují, je-li příslušné vlastní číslo jednonásobné. Samozřejmě, že je z hlediska numerického přístupu k řešení úlohy vhodné vektor před dosazením do nového kola iterací normovat na jedničku; zároveň je tak konvergence k vlastnímu číslu spojena s konvergencí takovéto normy výsledného vektoru iterací.

5.3 Podmíněnost matice Pro vlastní numerické řešení soustav lineárních rovnic je důležité číslo podmíněnosti matice, které odhaduje poměr relativní změny řešení ku relativní změně prvků matice a pravých stran příslušné lineární rovnice. Toto číslo také charakterizuje pracnost vyřešit soustavu. Definice

(podmíněnost

matice):

Číslo

𝜘(𝐴) = ‖𝐴‖/‖𝐴−1 ‖

nazveme

číslem

podmíněnosti matice. Zde ‖𝐴‖ je norma matice A. Nejčastěji se používá spektrální norma. Je-li A symetrická pozitivně definitní, je spektrální norma matice A rovna největšímu vlastnímu číslu matice A. Za dobře podmíněné matice se považují ty s podmíněností menší než 100. V prvním odstavci této podkapitoly se vlastně říká, jaký platí odhad a že se hodnota pro nějaké pravé strany nabývá (tedy že existuje taková kombinace vektorů).

5.4 Multigridní metoda Multigridní metody, taktéž jich existuje i v monografiích popsáno více variant, lze, jak v této práci bude prezentováno, považovat za inverzní operaci k zjemňování některých úloh, tedy místně- lokálně sloučení některých prvků. Toto lze učinit buď geometricky nebo algebraicky. Nebo se na slučování prvků můžeme podívat ještě tak, abychom úlohu na nejmenší vlastní číslo neovlivnili a aby prvky byly buď příbuzné, nebo jsme sloučením nepoškodili naši hledanou hodnotu. Principem multigridních metod je konstrukce menší úlohy, jejíž řešení je použito pro zpřesnění řešení původní úlohy. Řešit menší úlohu je méně náročné a tudíž lze očekávat, že se při vhodné volbě menší úlohy sníží celková výpočetní náročnost řešení. 23

V literatuře nalezneme multigridní metody, kde se hledaný vlastní vektor v každém cyklu opravuje aditivně, tedy nová aproximace vznikne přičtením opravy ke stávající aproximaci. V této práci uvedeme novou variantu multigridu. Pro variantu multigridní metody bylo odvozeno tvrzení o aproximaci (6.4.4) a konvergenci (6.6.12) viz [5].

6 Zvolené postupy Pro ověření vlastnosti úlohy konvergence k vlastnímu číslu stačí vzít poměr norem dvou po sobě jdoucích iterací, která musí konvergovat, podle části o lineární algebře, právě k nejmenšímu vlastnímu číslu. Z programátorského hlediska lze při vhodnosti velikosti matice 𝐴 volit takové blokové úpravy matice soustavy, že jsou pro dané uspořádání procesor–cache výhodnější, neprovádět slepé zápisy a dodržovat prostorovou a časovou lokalitu v kódu. Toto nelze ale vhodně zajistit vždy kvůli nutnosti sekvenčně-jdoucích-za-sebou kroků metody.

6.1 Řešení soustavy rovnic v inverzní iteraci pro problém vlastních čísel Samozřejmě, že z hlediska numerického přístupu k řešení úlohy je vhodné vektor před dosazením do nového kola iterací normovat na jedničku. V dalších částech zvolíme z hlediska numerické náročnosti kromě odhadu potřeby počtu operací i další přístup – budeme procesoru stopovat čas. Některé metody jsou pro typovou třídu úloh nejvhodnější a to i z hlediska konvergence, podle typu výsledku úlohy.

6.1.1 Gaussova eleminace Gaussova eliminace jako metoda řešení postupu inverzních iterací se sice řešila každý krok iterací v obecně řádově v

𝑁3 3

+

𝑁2 2

krocích. Ale pro naši volbu bazických funkcí

pro 2D úlohu se ukázalo, že to bude vzhledem k tvaru matice úlohy v asi 5,5. 𝑁 2 krocích operací (zde 𝑁 je rozměr matice). Z hlediska zaokrouhlovacích chyb bývá Gaussova eleminace vhodně doplňována částečnou (jen v jednom sloupci) nebo úplnou pivotací (podle maxima absolutních 24

hodnot vybírám i sloupec do eliminačních úprav) – pivotace je výběr v absolutní hodnotě největšího prvku a přeskupení matice soustavy výměnou řádku s tímto dominnantním řádkem. Většinou jsem používal částečnou pivotaci.

6.1.2 Choleskiho metoda Tato metoda je výhodná kvůli tomu, že na začátku můžeme naplnit matici rozkladu a dál jen iterovat výsledek, tedy přináší jistou úsporu a to už jen díky tomu, že nemusíme matici pro výpočet znovu plnit. Složitost řádově v 𝑁𝑁 krocích, pro naši matici to bude ještě méně.

6.1.3 Metoda sdružených gradientů Metoda neklade velké nároky na rezervaci paměti, za běhu potřebuje jenom jednorozměrná pole a v celku i rychle konverguje. Náročnost je dána výhradně lineárními operacemi nad vektory a minimálním využití násobení maticí. Přesné řešení poskytne teoreticky nejpozději po 𝑁 krocích, přesněji po 𝑚 krocích, kde 𝑚 je počet navzájem různých vlastních čísel matice 𝐴. Pro dobře podmíněné matice většinou mnohem dříve.

6.1.4 Použití multigridu Konstrukce hrubého prostoru Multigridní algoritmus [14], [15], [16], [17], v provedení v této práci, znovu předpokládá pozitivně definitní symetrickou matici. Algoritmus volby hrubého prostoru budeme charakterizovat jenom pro jeden krok algoritmu a to pro sloučení příbuzných indexů úlohy na nejmenší vlastní číslo. Při volbě hrubé úlohy postupujeme následovně:  Začneme s množinou všech indexů, odstraníme z ní prvky s výrazně dominantní diagonálou oproti součtu absolutních hodnot ostatních prvků.  V této množině zkonstruujeme množinu 𝑆𝑖 tvořenou těmi sloupcovými indexy z množiny indexů matice 𝐴 v 𝑖-tém řádku –mimo 𝑖, že splňují 𝑎𝑖𝑗 < −𝛽 max |𝑎𝑖𝑘 |. (Tato podmínka by šla nahradit i volnější podmínkou v absolutní 𝑎𝑖𝑘 <0

hodnotě) 25

 Vytvořme množinu 𝑚𝑗 = {𝑖|𝑗 ∈ 𝑆𝑖 }.  Položme 𝑛𝑐 = 0. Začněme iterativní algoritmus:  Vybrat 𝑖 že 𝑚𝑖 má minimální počet prvků, 𝑛𝑐 + +.  Vybrat 𝑗, že 𝑎𝑖𝑗 =

min

𝑘 je z množiny indexů

𝑎𝑖𝑘 . Eventuelně je v dalších odstavcích

vysvětleno, že by se dalo jako varianta tohoto postupu brát jednak prvky v absolutní hodnotě, nebo ještě např. ne nejmenší, ale druhý nejmenší prvek.  Jestliže je 𝑗 prvkem 𝑆𝑖 , pak 𝐺𝑛𝑐 = {𝑖, 𝑗}, jinak 𝐺𝑛𝑐 = {𝑖}.  Odebrat z množiny indexů prvky 𝐺𝑛𝑐 .  Pro 𝑘 ∈ 𝐺𝑛𝑐 pro všechny 𝑙 ∈ 𝑆𝑘 aktualizuj počet 𝑚𝑙 − −. Uzpůsobení heuristickým přístupem – konstrukce hrubého prostoru Alternativně bychom podle heuristického přístupu mohli brát nejdříve před všemi iteracemi stejný výběr 𝑚𝑖 a v něm brát záporné nediagonální elementy z indexové množiny v absolutní hodnotě malé, následně sloučit ty dvojice, které mají větší diagonálu- dojde tak s vysokou pravděpodobností k sloučení směrů pro ani velká a ani malá vlastní čísla. Tyto úvahy se zakládají na úloze pro existenci extrému formy 𝑥𝐴𝑥/‖𝑥‖2 - (První derivace i tá složka: λ𝑖 𝛼𝑖 (∑𝑗 𝛼𝑗 𝛼𝑗 ) − 2𝛼𝑖 (∑𝑘 λ𝑘 𝛼𝑘 𝛼𝑘 ), druhá derivace

i,

m

složka:

∑𝑗 𝛼𝑗 𝛼𝑗 (2λ𝑚 (∑𝑗 𝛼𝑗 𝛼𝑗 ) + 4λ𝑖 𝛼𝑖 𝛼𝑚 − 2 ∑𝑘 λ𝑘 𝛼𝑘 𝛼𝑘 −

4λ𝑖 𝛼𝑚 𝛼𝑚 ) − 2(∑𝑗 𝛼𝑗 𝛼𝑗 )𝛼𝑚 (2λ𝑖 𝛼𝑖 (∑𝑗 𝛼𝑗 𝛼𝑗 ) − 2𝛼𝑖 (∑𝑘 λ𝑘 𝛼𝑘 𝛼𝑘 )),

kde

se

vždy

přes některé indexy sčítá, rozhodující jen pro vlastní směry formy). Následovat může výše navržený multigrid. Postupovat můžeme takto ( viz [15]):  Začneme s množinou všech indexů, odstraníme z ní prvky s výrazně dominantní diagonálou oproti součtu absolutních hodnot ostatních prvků.  V této množině zkonstruujeme množinu 𝑆𝑖 tvořenou těmi sloupcovými indexy

z množiny indexů matice 𝐴 v 𝑖-tém řádku –mimo 𝑖, že 𝑎𝑖𝑗 < 0  Vytvořme množinu 𝑚𝑗 = {𝑖|𝑗 ∈ 𝑆𝑖 }. 26

 Položme 𝑛𝑐 = 0. Začněme iterativní algoritmus:  Vybrat 𝑖 že 𝑚𝑖 má minimální počet prvků, 𝑛𝑐 + +.  Vybrat 𝑗, že |𝑎𝑖𝑗 | =

min

𝑘 je z množiny indexů

|𝑎𝑖𝑘 |, v indexové množině máme jen

záporné elementy matice.  K 𝑗 najít 𝑘, že |𝑎𝑘𝑘 |není minimální, tyto seskupíme 𝐺𝑛𝑐 = {𝑖, 𝑘}, jinak 𝐺𝑛𝑐 = {𝑖}.  Odebrat z množiny indexů prvky 𝐺𝑛𝑐 .  Pro 𝑘 ∈ 𝐺𝑛𝑐 pro všechny 𝑙 ∈ 𝑆𝑘 aktualizuj počet 𝑚𝑙 − −. Tato multigridní metoda poškozuje hodnotu nejmenšího vlastního čísla méně než ta předchozí. Provedení multigridu Jako začátek iterací vezmeme libovolný počáteční vektor 𝑥0 , k němu najdeme jednou inverzní iterací vektor x, který umístíme do krajního sloupce matice 𝑃, ostatní sloupce matice tvoříme podle série množin 𝐺𝑛𝑐 . V 𝑘 + 1 ním sloupci matice 𝑃 jsou nenulové prvky v řádcích s indexy z množiny 𝐺𝑘 . Pak můžeme případně přidat ty, co nebyly obsaženy v indexové množině multigridu. Výslednou matici budeme označovat 𝑃(x). Na velikost získaného vl. čísla by nemělo mít vliv způsob jakým novou soustavu vytváříme. Přístup v této práci je víceméně algebraický a proto je třeba upozornit na to, že síla této metody by mohla spočívat v jejím řetězovém případně více-než-po-dvou seskupovacím postupu (další vlastnosti a rozšíření v [19] až [29]). Návrat do původního prostoru od hrubých multigridních iterací (s řešením 𝑔) se udělá znovu pomoci matice, dokonce platí, že nejmenší vlastní čísla hrubé úlohy 𝑃𝑇 𝐴𝑃 aproximují nejmenší vlastní číslo původní úlohy. Konkrétně 𝑃𝑔 aproximuje vlastní vektor k nejmenšímu vlastnímu číslu a to z důvodů, proč byla tato metoda heuristicky vybrána.Nejdříve

zvolíme

počáteční

aproximaci

𝑣0

𝑃𝑇 (𝑣𝑘 )𝐴𝑃(𝑣𝑘 )𝑧 =

𝑃𝑇 (𝑣𝑘 )𝑄𝑃(𝑣𝑘 )𝑧, vyřešíme a vrátíme se zpět k úloze 𝐴𝑥 = 𝑄𝑥, vztahem 𝑦𝑘 = 𝑃𝑧, provede se jeden krok inverzní iterace 𝑣𝑘+1 = 𝐴−1 𝑄𝑦𝑘 , vypočte se 𝜇𝑘 = ‖𝑣𝑘+1 ‖/‖𝑣𝑘 ‖,

27

srovnáme se zastavovacím kritériem a případně celý postup opakujeme (odhady vlastností v [14] (3.8) a v tam odkazovaných pracech). Vlastnosti varianty zvoleného multigridu s vkládáním vektoru pro dostatečný počet multigridních iterací – monotonie speciální posloupnosti aproximací nejmenšího vlastního čísla Definujme nejdříve pro matici 𝐴 tzv. Rayleighův podíl jako 𝑅(𝑢) = pro nejmenší vlastní číslo platí ekvivalent definice: 𝜆min = 𝑚𝑖𝑛𝑢≠0

𝑢𝑇 𝐴𝑢 𝑢𝑇 𝑢

𝑢𝑇 𝐴𝑢 𝑢𝑇 𝑢

. Dále

. Dále platí

vzhledem k inkluzi množin, přes které děláme minimum, že odhad nejmenšího vlastního čísla je menší než nejmenšího vlastního čísla získaného naší variantou multigridu (vzhledem k platnosti inkluze mezi prostory, přes které děláme minimum). Dále Rayleighův podíl vektrou po inverzní teraci je menší než jeho původní Rayleighův podíl (chceme:

𝑣 𝑇 𝐴𝑣 𝑣𝑇𝑣

𝑣 𝑇 𝐴−1 𝑣

≥ 𝑣𝑇 𝐴−2 𝑣, z Schwarzovy nerovnosti pro prostory se skalárním

součinem máme ale (𝑢𝑇 𝑢)2 ≤ 𝑢𝑇 𝐴−1 𝑢. 𝑢𝑇 𝐴𝑢 a (𝑢𝑇 𝐴−1 𝑢)2 ≤ 𝑢𝑇 𝑢. (𝑢𝑇 𝐴−2 𝑢), stačí aplikovat tuto na levou stranu s odmocninami (tyto existují, protože jsem volil vhodnou pozitivní matici) (𝑢𝑇 𝐴−1/2 𝐴1/2 𝑢), resp. pro skalární součin inverzního obrazu s původním vektorem). Pro přechod od řešení jedné iterace jemné úlohy k přesnému řešení hrubé úlohy platí 𝑅(jemné) ≥ 𝑅(přesně hrubé) (tj. 𝑅(z) ≥ 𝑅(P(y)z) (je tam totiž obsažen původní vektor). Důsledkem tohoto tvrzení je, že pokud se metoda přechodu k úplnému řešení hrubé úlohy doplní o jeden krok inverzní iterace „prodlouženého“ řešení na celém prostoru, můžeme opakovaném tohoto postupu dostat monotónní posloupnost odhadů vlastního čísla. Další vlastnosti multigridní metody jsou od odhadů norem [18] po další vlastnosti [10], či napojení na jiné algoritmy, či „špatně“ formulované úlohy.

28

6.2 Postup řešení úlohy 6.2.1 Příprava a ověření předpokladů Ověření předpokladů Všechny v části zvolené postupy, i za touto částí, splňují předpoklady koercitivity a spojitosti (ověřování tohoto předpokladu v [2]). Odvození koeficientů v konečně–prvkové bázi 𝜕 𝜕

𝜕 𝜕

Pro naši úlohu (1a), tj. (1b) v 2.1 bereme za základ úlohy 𝑎̃ 𝑥𝑥 𝜕𝑥 𝜕𝑥 𝑣 + 𝑎̃ 𝑦𝑦 𝜕𝑦 𝜕𝑦 𝑣 + 𝑥

𝜕 𝜕

𝜕 𝜕

𝑎̃ 𝑥𝑦 𝜕𝑥 𝜕𝑦 𝑣 + 𝑎̃ 𝑦𝑥 𝜕𝑦 𝜕𝑥 𝑣 = 𝜆𝑣 𝜕

𝜕

𝜕

𝜕

zjednodušenou 𝜕

𝜕

–

homogenní 𝜕

variační

úlohu:

𝜕

∫Ω 𝑎𝑥𝑥 𝜕𝑥 𝑢 𝜕𝑥 𝑣 + 𝑎𝑦𝑦 𝜕𝑦 𝑢 𝜕𝑦 𝑣 + 𝑎𝑥𝑦 𝜕𝑥 𝑢 𝜕𝑦 𝑣 + 𝑎𝑦𝑥 𝜕𝑦 𝑢 𝜕𝑥 𝑣 = 𝜆 ∫Ω 𝑢𝑣 . A budeme odvozovat indexované prvky v patřičných bazických funkcí pro levou a pravou stranu pro rovnici Ax = λQx.

Úvod volba bazických funkcí pro aproximaci úlohy Za bazické funkce zvolíme po částech bilineární funkce (z vhodného prostoru konečných prvků), podobné tzv. pyramidě, tj. bilineární. Uvažujme referenční obdélník s vrcholy 𝑣1 = (0,0), 𝑣2 = (ℎ𝑥 , 0), 𝑣3 = (0, ℎ𝑦 ) a 𝑣4 = (ℎ𝑥 , ℎ𝑦 ). Uvažujme čtyři tvarové funkce 𝜑𝑘 (𝑥, 𝑦), 𝑘 = 1,2,3,4, které jsou lineárními kombinacemi výrazů 1, 𝑥, 𝑦, 𝑥𝑦 a pro které platí 𝜑𝑘 (𝑣𝑗 ) = 𝛿𝑘𝑗 , kde 𝛿𝑘𝑗 , je tzv. Kroneckerovo delta, 𝛿𝑘𝑗 = 1 pro 𝑘 = 𝑗 a 𝛿𝑘𝑗 = 0 pro 𝑘 ≠ 𝑗. Uvedené funkce jsou bázovými (konečněprvkovými) funkcemi pouze na jednom elementu. Výsledné bázové funkce na celé oblasti Ω dostaneme složením těchto základních funcí tak, aby výsledné funkce byly spojité a nenulové na co nejmenší oblasti. Přesněji řečeno, zvolíme-li M x N vnitřních uzlových bodů, dostaneme M x N bázových funkcí. Každá tato funkce je nenulová pouze na čtyřech sousedních elementech (2 x 2) a její hodnota je 1 uprostřed této čtveřice elementů a v ostatních uzlech je rovna nule.

Příklady integrálů při výpočtu Uvedu příklady integrálů při výpočtu 𝐴, konst. v následujícím je konstanta, která vyjde v integrálu v závislosti na rozměrech elementu:

29

3

𝑦𝑖+1 𝑥𝑖+1 ∫𝑦 ∫𝑥 (𝑥 𝑖 𝑖

2

2

− 𝑥𝑖+1 )(𝑥𝑖+1 − 𝑥)/(ℎ𝑥 ℎ𝑦 )𝑑𝑥𝑑𝑦 =

𝑥 −𝑥 𝑑(𝑥 −𝑥 ) − 𝑖+1 𝑖+1 + 𝑖+1 𝑖 3

ℎ𝑥 2 ℎ𝑦 2

𝑦

𝑥

𝑦𝑖+1 2 −𝑦𝑖2

𝑖

𝑖

2

2 2 ∫𝑦 𝑖+1 ∫𝑥 𝑖+1 ((𝑦 − 𝑦𝑖+1 )(𝑥𝑖+1 − 𝑥))/(ℎ𝑥 ℎ𝑦 )𝑑𝑥𝑑𝑦 = (

−

𝑥𝑖+1 2 −𝑥𝑖 2 2

) /(ℎ𝑥 2 ℎ𝑦 2 ) =

konst.

3

1

= =𝑘𝑜𝑛𝑠𝑡. 3,

− 𝑑𝑦𝑖 ) (𝑑𝑥𝑖 −

− 1/4.

𝑦

𝑥

(𝑦𝑖+1 +𝑦𝑖 )(𝑦𝑖+1 2 −𝑦𝑖 2 )

𝑖

𝑖

2

2 2 ∫𝑦 𝑖+1 ∫𝑥 𝑖+1 (𝑦𝑖+1 − 𝑦)(𝑦 − 𝑦𝑖 )/(ℎ𝑥 ℎ𝑦 )𝑑𝑥𝑑𝑦 = ((−𝑦𝑖+1 𝑦𝑖 𝑑) + 1

1

−

1

− 3 ((𝑦𝑖+1 3 − 𝑦𝑖 3 )))/(ℎ𝑥 2 ℎ𝑦 2 ) =konst. (2 − 3). 𝑦𝑖+1

∫

∫

𝑦𝑖

𝑥𝑖+1 (𝑥

− 𝑥𝑖+1 )(𝑦𝑖 − 𝑦) ℎ𝑥 2 ℎ𝑦 2

𝑥𝑖

= Uvedu příklady

𝑥 2−𝑥 2 𝑦 2−𝑦 2 ( 𝑖+1 2 𝑖 − 𝑥𝑖+1 𝑑) (𝑦𝑖 𝑑 − 𝑖+1 2 𝑖 ) ℎ𝑥 2 ℎ𝑦 2

2ℎ𝑦 (2ℎ𝑦 −𝑦)

2ℎ𝑥

oblastech

se

2

0

stejným

výsledkem;

1

ℎ𝑥 ℎ𝑦 ℎ𝑥 ℎ𝑦

1

1

−2ℎ𝑥 2 )(−𝑦 2 + 3𝑦ℎ𝑦 − 2ℎ𝑦 2 ) 𝑑𝑥𝑑𝑦 = . (− 3

celkový

𝑥 (2ℎ𝑥 −𝑥)(2ℎ𝑦 −𝑦)(𝑥−ℎ𝑥 )(𝑦−ℎ𝑦 )

𝑦

ℎ𝑦 3 (8−1)

1 4

𝑑𝑦 = [𝑧 3 /3]0 ℎ𝑥 [𝑧 3 /3]ℎ𝑦 = 1/9. ℎ𝑥 ℎ𝑦 ℎ𝑥ℎ𝑥ℎ𝑦ℎ𝑦

Q potkávání rohem:∫𝑦 2 ∫𝑥 2

−2ℎ𝑥 3 )(−

= konst.

integrálů při výpočtu 𝑄 pro různé vzájemné polohy- sebe

sama:∫1ℎ𝑥 (2ℎ𝑥 − 𝑥)2 𝑑𝑥 ∫1ℎ𝑦 na 4

𝑑𝑥𝑑𝑦 =

3

ℎ𝑥 3 (8−1) 3 1

a to

výsledek 4/9. ℎ𝑥 ℎ𝑦 . 𝑦

𝑥

1

1

𝑑𝑥𝑑𝑦 = ∫𝑦 2 ∫𝑥 2(−𝑥 2 + 3𝑥ℎ𝑥 − 3

+ 2 ℎ𝑥 3 (4 − 1) − 1

+ 2 ℎ𝑦 3 (4 − 1) − 2ℎ𝑦 3 ) = (2 − 1/3)(2 − 1/3)ℎ𝑥 ℎ𝑦 ,

pro

směr

dolů doprava stejný výsledek. Q

překryv

šířkou

(2

2ℎ

2ℎ

oblasti):∫1ℎ 𝑥 (2ℎ𝑥 − 𝑥)(𝑥 − ℎ𝑥 ) /ℎ𝑥 ℎ𝑦 ℎ𝑥 ℎ𝑦 𝑑𝑥 ∫1ℎ 𝑦 (2ℎ𝑦 − 𝑥

2ℎ

𝑦

1ℎ𝑦

𝑦)(2ℎ𝑦 − 𝑦)𝑑𝑦 + ∫1ℎ 𝑥(2ℎ𝑥 − 𝑥)(𝑥 − ℎ𝑥 ) 𝑑𝑥 ∫0 𝑥

1

2)𝑑𝑥. [3 + +1/3] = (−

ℎ𝑥 3 (8−1) 3

3

𝑦 2 /ℎ𝑥 ℎ𝑦 ℎ𝑥 ℎ𝑦 𝑑𝑦 = ∫(−𝑥 2 + 3𝑥 − 1

+ 2 ℎ𝑥 3 (4 − 1) − 2ℎ𝑥 3 )2/3 = (2 − 1/3) 2/3ℎ𝑥 ℎ𝑦 .

Určení matic v problému 𝑨𝒙 = 𝝀𝑸𝒙 Zprvu bereme rozměry elementární buňky jako čtvercové rovné jedné. Pro matici 𝐴: pro překryv oblastí rohem: příklady bázových kombinací: der((2 − −𝑥)(2 − 𝑦))𝑑𝑒𝑟((𝑥 − 1)(𝑦 − 1)); 𝑑𝑒𝑟((𝑦 − 1)(2 − 𝑥))𝑑𝑒𝑟((𝑥 − 1)(2 − 𝑦)), 𝑑𝑒𝑟((3 − 𝑥)(2 − 𝑦))𝑑𝑒𝑟((𝑥 − 2)(𝑦 − 1)); s příspěvky označenými podle derivací: xx:1/3-1/2, yy:1/3-1/2, xy:-1/4, yx:-1/4; otočená situace (zprava doleva): xx:1/3-1/2, yy:1/3-1/2, xy:1/4, yx: ¼. 30

Pro překryv oblastí bokem: Příklady bázových kombinací: 𝑑𝑒𝑟((𝑥 − 1)(2 − −𝑦))𝑑𝑒𝑟((2 − 𝑦)(2 − 𝑥)),

𝑑𝑒𝑟((𝑦)(𝑥 − 1))𝑑𝑒𝑟((2 − 𝑥)(𝑦)),

𝑑𝑒𝑟((𝑥 − 2)(2 −

−𝑦))𝑑𝑒𝑟((2 − 𝑦)(3 − 𝑥)), dvě situace v jedné oblasti (prvky svisle nad sebou): 1ß) xx:-1/3, yy:-1/3+1/2, xy:1/4, yx:-1/4, 2ß) (na hranici oblasti – blíže k počátku souřadnic doleva)) xx: -1/3, yy:-1/3+1/2, xy: -¼, yx:1/4. Pro prvky vodorovně vedle sebe je xx a yy prohozeno, a pro xy a yx je směr nahoru původní vlevo. Pro

hodnotu

sama

se

sebou-

proti

směru

hodinových

ručiček:

1) xx:1/3, yy:1/3, yx: ¼, yx: ¼, 2) xx: (+)1/3,yy: (+)1/3, xy: -1/4, yx:-1/4, 3) xx: 1/3, yy: 1/3, xy: ¼, yx: ¼, 4) xx: (+)1/3, yy: 1/3, xy: -1/4, yx: -1/4. Poznámka:

-,|,x,

sebe

sama

jsou

vzájemné

geometrické

pozice

bází.

–: -a1/3+ a2(1/2-1/3)+a3/4-a4/4-a1hs/3+a2hs(1/2-1/3)-a3hs/4+a4hs/4; |: a1(1/2-1/3)-(a2)/3-a3/4+a4/4+a1hl(1/2-1/3)-(a2hl)/3+a3hl/4-a4hl/4; + směr pro druhou orientaci (v matici jdeme zpětně) prohazuje indexy 3 a 4, tedy smíšené derivace. X: a1(1/3-1/2)+a2(1/3-1/2) |-a3/4 |-a4/4(poslední dva prohodit pro prohození indexů bází)–(|-znamená

změnu

pro

pozici

dolů

doprava)

Sebe sama: a1(1/3)+a2(1/3)+a3( ¼)+ a4(¼) 2% a1(1/3)+a2(1/3)+a3(-1/4)+a4(-1/4) 3% a1(1/3)+a2(1/3)+a3(¼)+a4(¼)

4% a1(1/3)+a2(1/3)+a3(-1/4)+a4(-1/4).

Pro matici 𝐴 jsou koeficienty ještě u xx–ových derivací ℎ𝑦 /ℎ𝑥 , u yy –ových derivací ℎ𝑥 /ℎ𝑦 násobkem udaných hodnot. 4

Pro 𝑄: Sám se sebou: 9 ℎ𝑥 ℎ𝑦 . 1

1 2

X:ℎ𝑥 ℎ𝑦 (2 − 3)

1

1

+(pro dvě oblasti dohromady společný kus): ℎ𝑥 ℎ𝑦 2/3(2 − 3). Kde ℎ𝑥 je x-ový rozměr elementární oblasti a ℎ𝑦 je y-ový rozměr elementární oblasti. Při programování je výhodné zaplňovat matice 𝐴, 𝑄 tak, že postupujeme „po prvcích“ Všude při výpočtech koeficientů jsme aproximovali integrály s funkcemi 𝑎.. prostě jejich hodnotou uprostřed podoblasti, tedy jsme zvolili lokálně aproximací konstantou (předpokládáme dostatečně jemné rozdělení tak, aby aproximace byla dostatečná). Pro prohozené indexy je pro složitější úlohy vhodné koeficienty přepočítat.

31

Použití inverzní iterace Na řešení úlohy o nejmenším vlastním čísle používáme metodu inverzní iterace, kdy iterativně 𝑘 ⟶ 𝑘 + 1 řešíme znovu a znovu rovnici 𝐴𝑥𝑘+1 = 𝑄𝑥𝑘 . Vlastnosti tohoto postupu a konvergence byly probrány v části o aplikacích lineární algebry.

6.2.2 Vlastní provedení metod Známé přesné řešení pro jednu variantu zobecněného laplaciánu Na množině (0; 1) × (0; 1) má jednotkový laplacián nejmenší v absolutní hodnotě vlastní funkci sin(𝜋𝑥)sin(𝜋𝑦) k vlastnímu číslu −2𝜋 2 .

Konvergence multigridu Pro vhodnou volbu hrubého prostoru by měla tato metoda být rychle konvergentní. Konvergence výrazu 𝑥 𝑇 𝐴𝑥 − 𝑥 𝑇 𝑥 k nule je úměrná mocninám výrazu ℎ(2𝛼−1) , kde 𝐻 = ℎ𝛼 , jsou po řadě rozměry hrubé (𝐻) a jemné sítě (ℎ). Příkladem práce zabývající se konvergencí je např. [29].

6.2.3 Vlastní aplikace metod na některé konkrétní případy V praxi vznikají u tohoto typu úloh obrovské matice, je třeba proto použít vhodnou výpočetní techniku, metody uložení hodnot a vhodný numerický postup. Z těchto důvodů se tyto úlohy řeší na různých typech paralelních či distribuovaných počítačů, clusterů, dnes rozvíjejících se gridů, open projektů, linuxové free či placené knihovny. Ukazuje se také, že architektura některých počítačů a jejich procesorů je pro ty či ony úlohy vhodnější (např. matice na dřívějších Cray a pod.). S konkrétní problematikou svázané výpočty, např. v chemii [20] či subjaderné fyzice se dají mezi top 500 počítači najít konkrétní příklady (i jejich operační systémy a architektury).

6.2.4 Numerické výsledky experimentů a postupy při programování a ladění Iterace byly činěny většinou na síti 30 × 30 a pro řádově 1000 (pro špatně konvergující postupy) (až 10) iterací, nebo dosažení přesnosti nalezení vlastního čísla, tedy pro řádově jednotky iterací. 32

Postupy v ladění, rozdíly v distribucích, nastavení project properties ve Visual Studio a některé potíže a rysy procesorů a operačního systému důležité pro úlohu a další možnosti implementace algoritmu úlohy Podařilo se program spustit i pro větší rozměry, jenže zde nelze spustit najednou na jednom počítači vícero různých programů, tedy nejen doba omezuje experimentování s nastaveními ve Visual Studiu a různých jeho platforem. Z optimalizace používáme nastavení pro maximální rychlost (tuto volbu je třeba zadat ve Visual Studiu na dvou místech), dále pro přesnost v plovoucí čárce. Některé nastavení Visual Studia v project properties je nekompatibilní se zvolenou optimalizační volbou, Visual Studio dále umožňuje činit další volbou přímo z „commandline“- např. vypisování reportů o průběhu překladu- a v jedné verzi jsou instrukce navíc, které nejsou v pozdější, která má navíc zase jiné, např. AVX2 související s rozšířenými instrukcemi procesoru a jsou jen v rodině Haswell (dávají hlášku: „Source stopped working.“). Defaultní instrukce pro překlad jsou SSE2, maximum pro pole je 2GB, 64- bitové Visual Studio má větší nároky na jednu buňku typu long double než 32- bitové. Pro počítání s rozměrnějšími poli je třeba upravit velikost Stacku a Heap a to dle specifikací jednotlivých datových typů, které se do okénka v Linker nastavení v podzáložce Systém udávají v Bytech. Dále je potřeba zapnout podporu „Enable Large Address“. Poznámka: hodnota reserve musí být větší než commit, protože tam se jedná na prvním místě o virtuální a na druhém fyzickou paměť. V podzáložce C/C++ jsem dále musel upravit v záložce „General“ formát zpráv při ladění, tzv. „Debug information format“ na volbu /Zi; v podzáložce „Optimization“ jsem nastavil /O2 a /Ot pro největší rychlost běhu, volba „whole program optimization“, pro vícemodulární projekty, se nepodařila uspokojivě odzkoušet; v podzáložce „Code Generation“ bylo nutné změnit „Basic runtime checks“, zde lze trochu experimentovat multithreaded debugem, instrukční sady bylo na používaných počítačích pro rychlost běhu nejlépe nechat na defaultní hodnotě, zatímco pro větší přesnost jsem používal volbu „Floating Point Model“ /fp:precise. Nastavení pro provedení optimalizace apod. jsou však nekompatibilní s nastavením pro ladění programů ( rozdíl mezi release a debug verzí). Na msdn se nepodařilo najít vhodné volby pro kompilaci pro x64 kompilaci v novém prosředí VS 2013, či to školní verze nepodporovala. V různých počítačových učebnách byly dostupné různé možnosti funkčnosti překladače dané stářím strojů. 33

Zaznamenány byly dále případy, kdy na počítači s 8GB paměti nešly spustit s získáním správného výsledku na konci správně běžící verze s dynamickou alokací paměti, kdežto verze používající stack- statickou alokaci paměti- běžely správněji. Tyto potíže s dynamickým polem se opakovaly i po přenesení kódu a znovuzaložení projektu v VS 2013 i VS 2010 a byly vyřešeny kompletním přetypováním výrazů a zohledněním variant přetížených funkcí a precedence operátorů, protože překladač interpretoval některé reálné výrazy jako celočíselné. Podobně při rutinní práci se neobjevující chyba při ukládání zdrojového kódu je hůře sledovatelná a to jak díky přenášení verzí programu mezi 2010 a 2013, tak díky pravděpodobné možnosti chyby při posledním ukládání před ukončením VS – jako sice nepravděpodobná, ale chyba VS. Zkoušeny byly i verze vytvořeného programu zohledňující strukturou kódu velikost okénka cache, tzv. cache line či cache block, kdy se za sebe řadí instrukce s blízkými indexy, protože jsou do vyrovnávací paměti cache tedy přemísťovány po takto velkých blocích. Cache může být víceroasociativní (např. „2way-associative“), s přívlastkem „set“; tento parametr lze v aktuální verzi Visual Studia zjistit stažením z msdn programu, překompilováním a spuštěním, který po spuštění načte parametry procesorové cache, paměti a další údaje. Nejrychlejší byly programy pro pole o velikosti 10x10, protože takové matice, s jakými se pracuje, se pak vejdou celé do L2 cache. Pro větší musí procesor přetahovat data do vícestupňové vyrovnávací paměti a tyto následně zase přepisovat novými, což program zpomaluje. Rozdíl v rychlosti byl určen

přepočtením

na

jednu

operaci

řádově

zlepšení

o

asi

3

%.

Dostupné stroje měly cache line velikou 64 B (zjištěno spuštěním programu z msdn pro aktuální VS) a program zobrazoval několik typů cache pro každý jednotlivý využívaný procesor. Dále mají novější operační systémy v sobě zahrnuty vlastní optimalizační nástroje zajišťující rychlejší obsluhu konkrétních operací. Programy byly napsány tak, aby po ukončení výpočtu programem byl otevřen soubor, do něho zapsána doba běhu, velikost aproximace vlastního čísla, případně počet iterací do dosažení přesnoti pro jednotlivé typy iterací používané v programu a dále eventuálně po složkách ověřený podíl charakterizující vlastní čísla po složkách; zde je nutno poznamenat a zopakovat poznatky z přednášky, které jsou mj. dostupné na msdn pro systém windows, že pro tyto systémy nemůže název souboru obsahovat následující znaky, vždy za pomlčkou: -? –„ -\ -/ -< -> -* -| -: (jmenovitě otazník, uvozovky, lomítko, zpětné lomítko, znaménko větší, znaménko menší, hvězdičku, | a dvojtečku).

34

Pro tvorbu 64 bitové verze programu je potřeba změnit toto nastavení v nastaveních projektu, dále v configuration manageru, kde je potřeba dát pozor na to, aby se nevytvářela nová konfigurace, ale aby se tato překopírovala z té původní. Rovněž tak je potřeba v nastaveních projektu změnit „include“ a „library“ volbu na proměnnou, kterou VS používá pro 64 bitové knihovny. Konkrétní název knihovny *.lib nebo dll je potom

potřeba

uvést

v položce

Input

v nastaveních

linkeru.

VS se mezi sebou ve verzích 2010 a 2013 lišilo i nově zahrnutými parametry pro překlad a generování kódu „Code Generation“. Při ladění bylo použito #pragma specifická použitých překladačů a volby kompilátorů (gcc, VS icl kompiler, od výrobce procesorů). Po celou dobu byl používán interní VS překladač icl.exe- dostupný v instalačních adresářích VS; toto prostředí ale umožňuje v nastaveních projektu, či při doinstalaci používat jiný překladač jako další optimalizační možnosti, eventuálně překlad pomocí nově zahrnutých instrukčních sad. Podobně je tomu i u prostředí linux s grafickým IDE eclipse, které umožňuje taktéž použít pestrou škálu překladačů. Pro řešení úlohy se v obou prostředích dají použít rozšiřující knihovny pro paralelní vykonávání instrukcí, rovněž prostředí eclipse má svoji paralelní verzi, která umožňuje ladění paralelních programů. Konkrétní volby koeficientů jednoduché inverzní iterace Rozdílnost hodnot se dá vysvětlit tím, že byly aplikovány různé verze překladačů a postupy při kompilování programů. Co je určující pro vlastnosti metod, jsou počty iterací.

Volba 𝒂𝒙𝒙 = 𝒂𝒚𝒚 = 𝟏; 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = 𝟎; Pro Gaussovu eliminaci pro naši oblast 30x30 byl výsledek nalezen za 83962 ms při spouštění exe, a za 90573 ms při spouštění z VS pomocí volby debug start new instance (což je výhodné protože VS vypíše pro ladění potřebné chybové hlášení); hodnota při iterování do dosažení přesnosti byla určena zde přibližně jako 1/0,05~-2ππ (tj. podíl tam je kvůli užití inverzní iterace). Pro choleskiho metodu v pořádku 80ms v tabulce výsledků. Pro metodu konjugovaných gradientů nezapsalo VS na disk výsledek, i když býl kód v pořádku (výsledky pro přesnost 1e-10 652ms).

35

Obrázek 1 Výsledek pro koeficienty 1,1,0,0.

Volba koeficientů 𝒂𝒙𝒙 = 𝟏𝟎𝟎; 𝒂𝒚𝒚 = 𝟏; 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = 𝟎; Gaussova metoda pro 30x30 skončila za 88038 ms pro spuštění z adresáře debug, za 89038 ms při spuštění z VS a 100 iterací (1247820ms). Choleskiho metoda 2707ms, metoda konjugovaných gradientů 80 (2206/213) iterací za 5360 viz tabulky dále. Vlastní číslo bylo určeno iterováním do dosažení přesnosti přibližně jako 1/0,001 (tj. podíl tam je kvůli užití inverzní iterace) s tím, že pro tyto koeficienty je už konvergence horší a tudíž odhad není tak přesný.

Obrázek 2: Výsledky pro koeficienty 100, 1,0,0

Volba koeficientů 𝒂𝒙𝒙 = 𝟎. 𝟎𝟏𝒄𝒐𝒔𝟐 (𝟎. 𝟎𝟏) + 𝒔𝒊𝒏𝟐 (𝟎. 𝟎𝟏); 𝒂𝒚𝒚 = 𝒄𝒐𝒔𝟐 (𝟎. 𝟎𝟏) + 𝟎. 𝟎𝟏𝒔𝒊𝒏𝟐 (𝟎. 𝟎𝟏); 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = (𝟏 − 𝟎. 𝟎𝟏)cos(𝟎. 𝟎𝟏)sin(𝟎. 𝟎𝟏); Gaussova eleminace na 30x30 za 89175 ms při spuštění z VS a 84974 ms při spuštění exe souboru z debugu (177 iterací za 3507770ms). Užity též choleski 3105ms a konjugované gradienty 7000ms (87(2274) iterací). Vlastní číslo bylo určeno přibližně jako 0.015, nebo-li pro úlohu ze zadání práce 1/0.0995049 (tj. podíl tam je kvůli užití inverzní iterace).

36

Obrázek 3: Výsledná funkce pro koeficienty 𝟎. 𝟎𝟏𝐜𝐨𝐬 𝟐 (𝟎. 𝟎𝟏) + 𝐬𝐢𝐧𝟐 (𝟎. 𝟎𝟏); (𝟎. 𝟎𝟏) + +𝟎. 𝟎𝟏𝐬𝐢𝐧𝟐 (𝟎. 𝟎𝟏); (𝟏 − 𝟎. 𝟎𝟏)cos(𝟎. 𝟎𝟏)sin(𝟎. 𝟎𝟏); (𝟏 − 𝟎. 𝟎𝟏); cos(𝟎. 𝟎𝟏)sin(𝟎. 𝟎𝟏).

Srovnání základních variant metod Pro matici 10x10 byly všechny metody stejně náročné na čas i paměť i co do počtu operací (z hlediska úkonů pro procesor). Více se lišily potom u úloh s většími maticemi, kde pro úlohy s pozitivně definitní symetrickou maticí i co do nároků na paměť v celém rozsahu byly nejlepší metody konjugovaných gradientů (které ale potřebují, aby matice úlohy byla pozitivně definitní symetrická). Zaznamenány byly rozdíly mezi jednotlivými úlohami mezi koeficienty v první a druhé variantě, tedy 𝒂𝒙𝒙 = 𝒂𝒚𝒚 = 𝟏; 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = 𝟎; a 𝒂𝒙𝒙 = 𝟏𝟎𝟎; 𝒂𝒚𝒚 = 𝟏; 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = 𝟎; , kdy druhá varianta potřebovala více iterací do dosažení stejné přesnosti, tedy nekonvergovala tak rychle. V konvergenci byla na tom ještě hůře exotická volba koeficientů, která potřebovala dvaavícekrát více iterací. Konkrétní volby koeficientů pro konkrétní variantu multigridní metody Zpočátku byly metody zkoušeny nejdříve na síti 10 × 10 a pro malý počet iterací. V konečném provedení se ukázala na první volbu koeficientů lepší varianta metody z [15]. Vlastní multigrid byl lepší na druhou volbu koeficientů (tam kde se u multigridu z literatury bere nejzápornější pro slučování, my bereme největší záporný). Výběr druhého nejmenšího v algoritmu (tedy upravený multigrid z literatury, čili druhý vlastní, kde v kroku výběru nejmenšího pro slučování se bere druhý nejmenší) dle [15] měl méně jemných iterací u všech variant koeficientů. Varianty s vlivem na počty ve skupinách mohly mít několikero vyhazovacích postupů během slučování.

Multigridní metoda pro volbu 𝒂𝒙𝒙 = 𝒂𝒚𝒚 = 𝟏; 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = 𝟎; Multigridní metoda konvergovala 3386 s, tj. pro přepočet na větší síť z menší. Pro větší síť je rychlost konvergence metody v krát s zachycena v následující tabulce- při zapnutí dostupných optimalizací v Visual Studiu 2010 a pro malý počet iterací pro metodu 37

pracující se stackem,

hodnotově

jsou

výsledky

dány

použitými

nástroji

a nemusí být dosažitelné. Pro síť 20x20 se nejlepší ukázala vlastní multigridní metoda s výběrem při slučování od v absolutní hodnotě nejmenšího záporného čísla (5(8) (361, 182, 92 skupin), ostatní: multigrid z literatury 5(12) (361, 186 skupin), jeho verze s druhým nejmenším číslem v pořadí 5(11)) (361, 183, 94). Případně z hlediska ladění je vhodné kontrolovat počty vnitřní a vnější iterace metody i s např. výpisy úprav matic. Výsledky časů závisejí na buildování programu. Rozměr pole oblasti 50x50 vlastní multigrid 60x60 vlastní multigrid 70x70 vlastní multigrid 150x150

Doba řešení úlohy do dosažení přes. v krát [s]: 42400 47360 1d4690 3d4306

Multigridní metoda pro volbu koeficinetů 𝒂𝒙𝒙 = 𝟏𝟎𝟎; 𝒂𝒚𝒚 = 𝟏; 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = 𝟎; Pro tuto volbu konverguje iterační metoda pomaleji než u jednotkových koeficientů, pro upravený multigrid z literatury 5 hrubé a 212 jemných iterací (361, 199, 106 skupin), , pro vlastní multigrid 5 hrubých a 208 jemných (361, 182, 92 skupin), upravený multigrid z literatury 5 hrubé 218 jemných (361, 191 skupin); oproti 70 iteracím

Gaussovou

eliminací

na

úroveň

dosažení

stejné

přesnosti.

Totéž na síti 20x20 multigrid z literatury 5(276), vlastní multigrid 5(210) a druhý vlastní multigrid, tedy upravený z literatury 5(356).

Multigridní metoda pro exotickou volbu koeficientů 𝒂𝒙𝒙 = 𝟎. 𝟎𝟏𝒄𝒐𝒔𝟐 (𝟎. 𝟎𝟏) + 𝒔𝒊𝒏𝟐 (𝟎. 𝟎𝟏); 𝒂𝒚𝒚 = 𝒄𝒐𝒔𝟐 (𝟎. 𝟎𝟏) + 𝟎. 𝟎𝟏𝒔𝒊𝒏𝟐 (𝟎. 𝟎𝟏); 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = (𝟏 − 𝟎. 𝟎𝟏)cos(𝟎. 𝟎𝟏)sin(𝟎. 𝟎𝟏); Multigridní metoda s Gaussovou metodou na 30x30 7045 s a 2307 ms na malé síti 10x10, viz tabulky. Opět zkoušeno pro pevný počet iterací a s pevným rozdílem a fixní přesností pro porovnání metod. Na síti 20x20 multigrid z literatury 5(72) (při vyházení více indexů z dalšího slučování 5(137)) (asi 360 a 191 skupin; jinak 361 a 264 skupin), vlastní multigrid 5(173) (při vyházení více indexů z dalšího slučování 5(75)) (361,182,92 skupin; asi 360 a asi 190 skupin), druhý vlastní multigrid upravený z literatury 5(178) (při vyházení více indexů z dalšího slučování 5(155)) (361, 215, 129; jinak 361, 191, 101 skupin). Pro všechny varianty multigridní metody tato úloha konvergovala ze všech koeficientů nejpomaleji a hodnoty odhadu nejmenšího vlastního čísla se i nejvíce mezi metodami 38

lišily, i když rozdíly konvergovaly stejným způsobem (tj. se stejnou ukončovací konstantou). Konvergence je zachycena v počtu iterací, protože tato konstanta byla dána všem stejná.

Obrázek 4: Multigridní metoda pro exotickou volbu laplaciánu

Obrázek 5: Výsledek vlastního multigridu pro exotickou volbu laplaciánu pro malý počet iterací.

6.2.5 Doby běhů pro pevně stanovený počet iterací Používáme zkratky kvůli velikosti buněk tabulky vl. grid, tj. multigrid braný od v absolutní hodnotě nejmenšího záporného čísla s upřednostněním těch s větší diagonálou (vcelku stejný počet skupin jako následující); mgrid, tj. multigrid z literatury (tj z [15]) (braný od nejzápornějšího - vcelku stejný počet skupin jako předešlý); vl. druhý grid, tj. varianta multigridu z literatury, když se vezme ne nejzápornější číslo, ale druhé nejzápornější číslo v příslušném řádku matice v algoritmu multigridu (ten byl nejlepší z hlediska velikosti skupin). Ještě bylo učiněno několik ojedinělých pokusů s iterováním na každé úrovni multigridu, či tvorbou vícero slučovacích zón v matici a ovlivňováním následujícího kroku algoritmu jen tím, že prvek vymažu jen z 𝑚. Tedy ne z 𝑈. (viz pasáž o multigridu), což umožňuje zvětšit počet skupin, s dopadem na vlastnosti algoritmu. U počtu iterací je dodržen formát hrubé (pak jemné). Visual Studio 2013 Výsledky jsou na síti 30x30, místy i na 10x10. Verze paměti stack stack stack

dle Počet iterací

Koef.

Volby VS

1000 1000

1,1,0,0 1,1,0,0

100/1000

1,1,0,0

bez Intrin. optimal Bez

metoda Gauss. El. Fce, Gauss. El. Choleski

Doba běhu [ms] 483900 476859 15866

39

stack

1000

1,1,0,0

Heap+stack

1000

1,1,0,0

70x70, stack Stack Heap+stack stack

1000 100/1000 1000 1000

1,1,0,0 1,1,0,0 100,1,0,0 100,1,0,0

Heap+stack

1000

100,1,0,0

Stack+heap Stack+heap Stack+heap

10000/1000 10000/1000 10000/1000

1,1,0,0 1,1,0,0 1,1,0,0

Heap+stack

1000

100,1,0,0

Stack

1000

exotické

Stack+heap

10000/1000

1,1,0,0

stack heap

1000 6 (1)…

1,1,0,0 1,1,0,0

Sse2, intrin., optimal Sse2, intrin, optimal. Bez Avx, optimal Avx, optimal Intr., avx, optimal Intr., avx, vectorizér, optimal Bez Intr., optimal Intr., sse2, optimal Intr., avx, vectorizér, optimal Intr., avx, vectorizér, optimal

Gauss. El.

480068

Gauss. El

448552

Choleski Choleski Gauss. El. Gauss. El

1d6737210 16006 444699 575761

Gauss. El

444824

Vl. multigrid Vl. multigrid Vl. multigrid

13670300 13537000 13532500

Gauss. El.

446412

Gauss. El.

581977

Intr., sse2, Vl. multigrid optimal Intr, avx, p Konj. Grad. Optimal., intr, Vl.grid G. bez; s par. eliminace

13534800 4479070 13626500; 13512400 (c)

Hodnoty v tabulece jsou doby běhu v ms, Heap označuje alokaci proměnných na haldě, Stack v zásobníku, přičemž kombinace obou může znamenat to, že se metoda nastartovala počátečnímu hodnotami ze stacku a potom už se v cyklech programu jenom pracovalo s proměnnými alokovanými na haldě, tedy v C++ s dynamicky alokovanými proměnnými pomocí konstruktu s operátorem new. Obvyklé na 30x30. Visual Studio 2010 Výsledky jsou na síti 30x30, místy i na 10x10. Verze paměti heap Heap Heap Heap Stack Stack Heap stack heap heap

dle Počet iterací 1000 10000 1000 1000 100 100 100 100 100 100

koeficienty

Volby VS

metoda

100,1,0,0 Exot. Exot. 100,1,0,0 100,1,0,0 100,1,0,0 Exot. Exot.

optimal Optimal, intr Optimal, intr Optimal, intr Optimal, intr Optimal, intr Optimal, intr Optimal, intr Optimal, intr Optim., intr.

Gauss. El. Gauss. El. Gauss. El. Gauss. El. Choleski Gauss. El. Gauss. El. Choleski Gauss. El. G. eliminace

1,1,0,0

Doba běhu [ms] 1247820 12488100 1120270 1118020 2707 122821 100006 2677 123405 134387

40

heap heap stack heap stack Heap+stat heap Heap+stat Heap+stat heap+stat

6 (1)… 1000 100 1000 1000 70 100

1,1,0,0 1,1,0,0 exot 1,1,0,0 exot 100,1,0,0 exot 1,1,0,0 100,1,0,0 exot

Optim., intr.

Vl.grid G. El G. eliminace Choleski G. eliminace Choleski G. eliminace G. eliminace G. eliminace G. eliminace G. eliminace

22315800 7922970 87 930142 3105 31233 448994 462451 499566

6.2.6 Parametry a výsledky při iterování do dosažení přesnosti 0.00001 Zastavovacím kritériem pro tuto kapitolu v programech pro většinu metod bylo |𝜆(𝑗) − 𝜆(𝑗−1) | < 0.00001. |𝜆(1) − 𝜆(2) |, kde zde 𝜆 znamená nejmenší vlastní číslo v iteraci dané indexem, tedy samozřejmě bráno vždy jako aproximace pomocí inverzní iterace Galerkinovy aproximace tohoto nejmenšího vlastního čísla (0,00001=1e-5). Visual Studio 2013 Některé výsledky z následující tabulky mohou být zvýšeny díky tomu, že kompilace a build byl udělány na starších procesorech a exe bylu spuštěno na novějším s rozšířenými instrukcemi. V. dle paměti 10x10 heap

Počet iterací 8 (2)

Koef. 1,1,0,0

30x30 heap

1 (88)…

1,1,0,0

30x30 heap 30x30 heap

56 1(15)…

100,1,0,0 1,1,0,0

30x30 heap

1(15)…

1,1,0,0

100x100 heap

8

1,1,0,0

100x100 heap 50x50 heap

8 93

1,1,0,0 100,1,0,0

Volby VS Intrin.opt avx, p Intrin.opt avx, p Intr.o. avx p Intrin. avx, p

metoda Vl. grid

Doba [ms] 2836

mgrid

2955420

G. eliminace Mgrid G. eliminace Intrin. avx, p Mgrid G. eliminace Intrin., avx, G. eliminace p. Intrin., p. G. eliminace Intrin., p G. eliminace

25833 6015720 5818170 8637330 8961310 1474760

41

Visual Studio 2010 Obvyklé na síti 30x30. Číselné parametry ve sloupci Volby VS znamenají pro koeficienty pevně stanovený rozdíl pro všechny metody bez rozdílu- jako zastavující kritérium. V. dle pam./ přes. 30x30 h. 30x30 h. H. 30x30 H 30x30 50x50 H přes. 10−9 přes. 10−9 100x100h. H.100x100 100x100h. 10x10 h. 10x10 h. 10x10 h. 30x30h. 30x30h. 10x10 30x30h. 30x30h.

30x30 h. 30x30 h. 30x30 h. heap Heap heap heap Heap Heap Heap

Počet iterací 10 (2) 1 (14) 1 (14) 6 (1) 177

Koef.

169 7 11(2) 11(2) 13(10) 7(3) 7(11) 1(10) 11(72) 13(10) 1(7) 7(147) 11(152) 11(130) 1(10) 3(58) 10(152) 4(147) 3(71) 11(72) 13(10) 11(152) 13(10) 5(147) 3(58) 11(130) 33(10) 11(72) 1(7) 33(10) 13(10) 4 72 71 10(15139) 7(187) 7(67) 10(118)

100,1,0,0 1,1,0,0 1,1,0,0

1,1,0,0 1,1,0,0 100,1,0,0 1,1,0,0 exot

1,1,0,0 1,1,0,0 100,1,0,0 1,1,0,0 1,1,0,0 exot 1,1,0,0

Volby VS/rozdil intrin intrin Optim., intr. Optim., intr.

0,0679044 0,0026 0,00002 0,0367324

0,0035 0,0000243 0,01 826 0,00000237 0,01

1,1,0,0 exot

826 841 838

100,1,0,0 1,1,0,0 1,1,0,0 1,1,0,0 1,1,0,0 100,1,0,0 exot 1,1,0,0 1,1,0,0 1,1,0,0 1,1,0,0

0,01 0,01 0,01

1e-10(2e-15) 1e-5(2e-10) 1e-5(2e-3) 1e-10(2e-5)

Metoda

Doba b. [ms]

Vl. grid Mgrid G. E. Mgrid G. E.. Vl.grid G. E.. G. eliminace

12088700 32125900 80672900 22284700

Choleski G. eliminace Vl.grid G. E. Vl.grid G. E. Mgrid G. E. Vl.grid G. E. Vl.druhý gr. G. E.

200 2d10522900 1d10200500 1d10535255 2007 33500000

Vl.grid G. E. Vl.grid G. E.. Vl.druhý gr. G.E. Vl.grid G. E. Vl.grid G. E. Vl.grid G. E. Vl.druhý gr. G. E. . Mgrid G. E. Multigrid Vl.druhý gr. G. E.

Vl.grid G. E. Vl.druhý gr. G. E. Vl.druhý gr. G. E. Vl.grid G. El. Mgrid G. E. Vl.druhý gr. G. E. Mgrid G. e. Vl.grid G. E. G. eliminace G. eliminace G. eliminace Konj. Grad. Konj. Grad. Konj. Grad. Konj. Grad.

11000000 33000000 1404 10080000 10290000 37000000 12300000 13900000 2891 33626000 10360000 14370000 10600000 14167000 27860000 10290000 37843000 14800000 40043800 34121000 14167000 39040 541230 1085780 137870 559 235 374

42

heap 20x20 hi 20x20 hi 20x20 hi 20x20 hi 20x20 hi 20x20 hi 20x20 hi 20x20 hi 20x20 hi

10(205) 5(12) 5(8) 5(11) 5(356) 5(210) 5(276) 5(72)/ (137) 5(173)/ (75) 5(178)/ (155)

1,1,0,0 1,1,0,0 1,1,0,0 1,1,0,0 100,1,0,0 100,1,0,0 100,1,0,0 Exot

1e-10(2e-10) 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5

Konj. Grad. Mgrid g.e. Vl.gr.g.e. Vl.druhý gr. G.e. Vl.druhý gr. G.e. Vl.gr.g.e. Mgrid g.e. Mgrid g.e.

652 12578 11984 11984 34471 23266 28251 41080/ 75410

Exot

1e-5

Vl.gr.g.e.

84660/ 25376

Exot

1e-5

Vl.druhý gr. G.e.

27267/ 25720

Koeficienty 𝒂𝒙𝒙 = 𝟏; 𝒂𝒚𝒚 = 𝟏; 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = 𝟎;.parametry běhů programů rozměr 10 20 40

g.e.čas [ms] 15

Počet it.

k.g. čas [ms] 0 0 1313

6 6 6

62831

Počet it. 7 7 7

Vnitřní it.

It.

Skup.

5(8) 5(8) 5(8)

22 92 382

Choleski čas [ms] 0

33 60 114

Počet it. 5 5…

Multigridní metody: rozměr 10 20 40

Gr.lit [ms] 93 23156 5.37e6

Iter.

Skup.

5(9) 5(12) 5(10)

43 186 771

Vl.gr. [ms] 78 18360 5.558e6

Druhý g. [ms] 78 20110 8.589e6

Iter.

Skup.

5(8) 5(11) 5(10)

24 94 384

Koeficienty 𝒂𝒙𝒙 = 𝟏𝟎𝟎; 𝒂𝒚𝒚 = 𝟏; 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = 𝟎;.parametry běhů programů Základní verze metod: rozměr 10 20 40

g.e.čas [ms] 46 4187 333955

Počet it.

k.g. čas [ms] 25173 656 25173

44 76 77

Počet it. 77 72 77

Vnitřní it. 2700 1465 787

Choleski čas [ms]

Počet it.

Multigridní metody: rozměr 10 20 40

Gr.lit [ms] 296 48643 6.673e6

Iter.

Skup.

5(246) 5(276) 5(269)

46 191 781

Vl.gr. [ms] 250 30064 9.6752e6

It.

Skup.

5(208) 5(210) 5(210)

382 92 22

Druhý g. [ms] 296 42314 1.056e7

Iter.

Skup.

5(222) 5(356) 5(279)

30 106 414

Koeficienty exotické 𝒂𝒙𝒙 = 𝟎. 𝟎𝟏𝒄𝒐𝒔𝟐 (𝟎. 𝟎𝟏) + 𝒔𝒊𝒏𝟐 (𝟎. 𝟎𝟏); 𝒂𝒚𝒚 = 𝒄𝒐𝒔𝟐 (𝟎. 𝟎𝟏) + 𝟎. 𝟎𝟏𝒔𝒊𝒏𝟐 (𝟎. 𝟎𝟏); 𝒂𝒙𝒚 = 𝒂𝒚𝒙 = (𝟏 − 𝟎. 𝟎𝟏)cos(𝟎. 𝟎𝟏)sin(𝟎. 𝟎𝟏);parametry běhů programů Základní verze metod: rozměr 10 20 40

g.e.čas [ms] 46 4234 445538

Počet it.

k.g. čas [ms]

Počet it.

72 67 77

15(328) 82() 687(1500) 80() 898475(328) 78()

Vnitřní it. 3084 1480 5694

Choleski čas [ms]

Počet it.

43

Multigridní metody: rozm ěr 10 20 40

Gr.lit [ms]

Iter.

125 20610 4.9494e6

6(99) 5(72) 5(72)

Skup . 46 191 781

Vl.gr. [ms] 171 29798 9.0444e6

It. 6(193) 5(173) 5(153)

Sku p. 92 382

Druhý g. [ms] 250 29814 9.08311e 6

Iter. 5(199) 5(178) 5(81)

Skup . 28 101 404

Pro tyto koeficienty byl adaptován program konjugovaných gradientů s ohledem na tvar matic na ukládání a ukládání hodnot a operace jen s nenulovými prvky matice. Pro síť 320x320 prvků dávala tato metoda dobu běhu 292927 s odchylkou na 4. platném desetinném místě.

Obrázek 6: Multigrid pro 3 hrubé iterace z literatury [15], vlastní multigrid, upravený multigrid z literatury; všechny pro případ axx = 1; ayy = 1; axy = ayx = 0;.

Obrázek 7: Multigrid pro 3 hrubé iterace z literatury [15], vlastní multigrid, upravený multigrid z literatury; všechny pro případ 𝐚𝐱𝐱 = 𝟏𝟎𝟎; 𝐚𝐲𝐲 = 𝟏; 𝐚𝐱𝐲 = 𝐚𝐲𝐱 = 𝟎;.

Obrázek 8: Pozměněný multigrid pro 3 hrubé iterace z literatury [15], vlastní multigrid, upravený multigrid z literatury; všechny pro případ 𝐚𝐱𝐱 = 𝟏𝟎𝟎; 𝐚𝐲𝐲 = 𝟏; 𝐚𝐱𝐲 = 𝐚𝐲𝐱 = 𝟎;.

44

Nejlepší výsledky vybraných verzí programů Nejlepší výsledky byly zaznamenány pro Gaussovu eliminaci pro 30x30 (tedy pole v programu 294) 100 iterací (70 iterací už má přesnost asi 1e-5) pro druhou volbu koeficientů na 12 až 18 sekund použitím pokročilých optimalizačních postupů programovacího prostředí a kompilátoru a 40 ms pro konjugované gradienty pro tutéž úlohu (méně spolehlivých pokusů pro stejnou přesnost); při tomto byl učiněn pokus využít intrinsické a další rozšířené instrukce z palety defaultních knihoven VS a dalších, bylo to však neúspěšné – vše na procesorech i5 třetí generace s 8GB pamětí (i na počítačích a programátorských prostředích s jinými parametry, ne však takové dobré hodnoty). Metoda konjugovaných gradientů se ukázala nejlepší pro pevně stanovenou vhodnou přesnost; pro zmenšující se přesnost nebo hodnotu zastavujícího kritéria začíná počet kroků metody prudce narůstat. Celé experimentování s multigridními metodami bylo zajímavé, protože se dá najít ještě více obměn multigridu z literatury; snad jen u vícestupňového multigridu iterování vždy na každé úrovni nepřineslo úsporu jemných iterací v závěru. Zkoušeny byly ojediněle i obměny podle výběru koeficientů, které se slučovaly, zavedením vícero zón slučování (v matici i v oblasti), dále změnou podmínky, kdy po sloučení se koeficient odstraní jen z příslušné indexové množiny namísto z množiny vhodných indexů. Zkoušeny byly i některé další algoritmické obměny programů a reprezentace dat. Závěry byly činěny jak pro techniku kompilace a programování, tak pro chování jednotlivých metod podle koeficientů, velikosti sítě, metody apod.

45

7 Závěr Práce byla zajímavým seznámením s počítačovými a matematickými metodami nejenom numerického řešení některých úloh z praxe. Při studiu literatury jsem se seznámil i s teorií obecnějších úloh než homogenní a i tzv. velmi slabých řešení, což bylo velmi zajímavé. Vidět např. propojení vlastností bilineární formy a jejího operátoru, případně jejich a odvozených funkcí a jejich diferencovatelnosti. Při přípravě programátorského řešení jsem se seznámil se širokou paletou manuálů využitelných knihoven a to od těch, co bylo třeba instalovat na propojovacím serveru do active directory domény, tak po nové normy a doporučené postupy v C++/C programování. Bylo potřeba překlenout potíže typu, kdy nefungoval pro účely ladění v jazyce C zápis do souboru, docházelo k „Stack overflow“. Tyto se podařilo vyřešit překopírováním kódu do programu v nově založeném projektu, kde bez potíží fungoval, jak měl, alternativně nastavením ve Visual Studiu. Při tomto ladění skákal debugger do dissasembly, což se stalo několikrát i při pokusu využít větší počet indexů a tím i virtuální paměť. Při spouštění programů multigridními metodami se někdy nepodařilo spustit úlohu na stejně velkých maticích jako obyčejnou Gaussovu eliminaci, protože se tam konstruuje řada dalších objektů (kvůli využití paměti). Podobně při „únavě“ počítač nevykonal na konci výpočtu uložení do souboru- při přechodu na jiný typ stroje se toto později přestalo objevovat, dále bylo někdy třeba zvětšit Stack v nastavení projektu „project properties“ a zapnout podporu rozšířeného adresování, někdy se projekt ve VS 2013 zablokoval do té míry, že nebylo možné otevřít vlastnosti projektu, které se zobrazovaly jako prázdné, a bylo nutné restartovat procházení vlastností vývojového prostředí. Pro zobrazení výsledků byl vhodným nástrojem program scilab, jeho schopnost načítat textový soubor do matic a následně je zobrazovat. Podobně bylo zajímavé srovnat výsledky pokusů z počítačů s různými parametry. Při studiu literatury jsem nalezl propojení symetrických forem a dobře řešitelných obecně pojatých problémů variačního počtu. Vše je poměrně obšírně pojednané v literatuře. Pro vlastnosti objektů lineární algebry stejně jako diferenciálních rovnic jsem čerpal z přednáškových materiálů VŠPJ. 46

Volba diferenciálního operátoru byla dostatečně obecná ze spektra vhodných podobně řešitelných úloh na jednoduchá vlastní čísla. V průběhu práce byly citovány časopisecké zdroje skýtající další možnosti použití a zobecnění metod. Multigridní metoda a jí podobné a další obdoby operací s oblastmi jsou možnostmi další orientace. Stejně jako bylo při rešerších objeveny kombinace těchto metod, např. multigrid doplněný o zjemňování u nespojitosti nebo singularity; úlohy se opravdu mohou řešit tak, že se toto místo vyjme a obejde, nebo se zvolí obecnější formulace problému. Nabyl jsem praktické zkušenosti při optimalizačních a jiných volbách kompilace ve Visual Studiu i např. při pokusech spouštět kódy kompilované na starších typech procesorů přímo z operačního systému novějšího data a širší instrukční sady. Stalo se i, že programy se nedaly spustit, jindy spustit šly, pravděpodobně to bylo způsobeno nejenom přetížením počítačů, na kterých jsem pracoval, ale i tím, že někde nebylo nainstalované to či ono potřebné dll, jednalo se o starší počítače bez některých nověji zařazovaných instrukcí apod. Vhodnou oporou při programování funkcí metod a korekci

programů

byly

dále

přednáškové

materiály

VŠPJ

z celé

palety

programovacích předmětů a konzultování nefunkčností s pracovníky VŠPJ katedry KEI. Dále mohly některé nepřesnoti vzniknout při přepisování dokumentu, protože se např. změnily parametry metod, přechodem na jiný způsob zápisu parametrů metody apod. Zajímavé bylo využívat teoretické poznatky pro psaní rychlejších programů se zohledněním

chování

procesorů pro

různě velká pole

a přesvědčit

se,

že např. pro rozměr o velikosti 10x10 jsou běhy programů nejrychlejší, protože se matice celé vejdou do L2 cache a procesor tudíž pravděpodobně nemusí přetahovat a vracet z a do paměti RAM hodnoty počítaných veličin. Přece jenom došlo i k chybám při přenášení a tvorbě programů, VS znehodnotilo jeden flashdisk, další odešel za rutinní činnosti, tudíž některé výsledky bylo nutno přehodnotit a znovu vytvořit nebo znovu provést úpravy programů, protože některé výsledky byly ztraceny.

47

Práce poskytla zajímavý pohled na informatiku a další obory včetně jejich aplikací. Multigridní metody šetří čas redukcí úlohy na menší a využitím těchto výsledků dále; nejlepší z hlediska velikosti sítě se ukázala upravená varianta z literatury s výběrem druhého nejmenšího záporného prvku, případně s omezením vyhazovaní prvků při slučovacím postupu; z hlediska hůře konvergujících úloh se metoda z literatury dostala až na druhé místo za vlastní multigrid s primárním výběrem největšího záporného prvku z příslušné množiny (ten navíc upřednostňuje prvky, které mají větší diagonální prvek), i když vcelku poskytla taky za malý počet iterací výsledek a v počtu iterací byla dokonce jedna její varianta i lepší. Stěžejními partiemi této práce jsou pasáže zabývající se provedením a náročností různých variant metod a algoritmů; možnosti další orientace jsou v odkazované literatuře.

48

Seznam použité literatury 1. Siegfried, C.; Le, V. K.; Motrenu, D. Nonsmooth variational problems and their inequalities: comparison principles and applications. Springer, 2007. ISBN-13: 978-0-387-30653-7. 2. RNDr. Rektorys, K., DrSc. Variační metody v inženýrských problémech a v problémech matematické fyziky. 6. vyd., Praha: Academia, 1999. ISBN 80200-07148-8. 3. Brenner, S. C.; Scott, L. Ridgway. Direct methods in the theory of elliptic equations. 3. vyd. Springer, 2012. ISBN 978-0-387-75933-3. 4. Nečas, Jindřich. Direct methods in the theory of elliptic equations. Springer, 2012. ISBN 978-3-642-10454-1. 5. Dr. Sauter, S. Finite elements for elliptic eigenvalue problems. Lecture notes for the Zürich Summer School [online]. 2008. [cit.: 30.7.2014 ] dostupné na: http://www.math.uzh.ch/conferences/uploads/media/Sauter1.pdf 6. Auchmuty, G. Variational principles of eigenvalues of compact operators. Siam Journal of mathematical analysis. 1989, 20, č. 6, s. 1321-1335. DOI: 10.1137/0520087. 7. Auchmuty, G. Variational principles for eigenvalues of nonsymmetric matrices. Siam Journal of matrix analysis appl. 1989, 10, č.1, s. 105-117. DOI: 10.1137/0610008. 8. Sangalli, G. A uniform analysis of nonsymmetric and coercive linear operators. Siam Journal of mathematical analysis. 2005, 36, č.6, s. 2033-2048. DOI: 10.1137/S0036141008434996. 9. Chatelin, F. The spectral approximation for linear operators with applications to the computation of eigenelements of differential and integral operators. Siam review. 1981. 23, č. 4. DOI: 10.1016/0378-4754(86)90031-5. 10. Arnold, D. N.; Falk, R. S.; Winther, Ragnar. Finite element exterior calculus, homological techniques, and applications. Acta Numerica. 2006, s. 1-55. DOI: doi: 10.1017/S0962492906210018. 49

11. Saberi Najafi, H; Khaleghi, E. A new restarting method in the Arnoldi algorithm for computing the eigenvalues of a nonsymmetric matrix. Applied mathematics and computation. 2004, I, č. 156, s. 59-71. DOI: doi 10.1016/j.amc.2003.07.018. 12. Weiss, R. Properties of generalized conjugate gradient methods. Numerical linear

algebra

with

applications.

1994,

I,

č.

1,

s.

45-63.

DOI:

10.1002/nla.1680010106. 13. Eisenstatt, S.C.; Elman, H. C.; Schultz, M. H. Variational iterative methods for nonsymmetric systems of linear equations. Siam Journal of numerical analysis. 1983, 20, č.2, s. 345-357. DOI: 10.1137/0720023. 14. Fraňková, P.;Hanuš, M.; Kopicová, H.; Kužel, R.;Marek, I.; Pultarová, I.; Vaněk, P.; Vastl, Z. Convergence theory for the exact interpolation scheme with approximation vector used as the first column of the prolongator: The partial eigenvalue problém. [Odesláno k publikování]. 15. Notay, Y. An aggregation based algebraic multigrid method. Electronic Transactions on Numerical Analysis. 2010, č. 37, s. 123-146. 2010, Kent State University. ISSN: 1068-9613. 16. Brandt, A.; McCormick, S.; Ruge, J. Multigrid methods for differential eigenproblems. Siam journal of scientific computing. 1983, II, č. 4, s. 244-260. id. číslo: 0196-5204/83/0402-0009 $01.25/0. DOI: 10.1137/0904019. 17. Cai, Z.; Mandel, J.; McCormick, S. Multigrid methods for nearly singular linear equations and eigenvalue problems. Siam journal of numerical analysis. 1997, I, č. 34, s. 178-200. id. číslo: PII. S1064827594261139. DOI: 10.1137/S1064827. 18. Bornemann, F.; Yserentant, H. A basic norm equivalence for the theory of multilevel methods. Numerical mathematics. 1993, 64, s. 455-476. DOI: 10.1007/BF01388699. 19. Zhang, S.; Yu, D. Multigrid algorithm for the coupling system of natural boundary element method and finite element method for unbounded domain problems. Journal of computational mathematics. 2007, 25, č. 1, s. 13-26.

50

20. Ed.: Grotendorst, J.; Attig, N.; Blügel, S.; Marx, D. Multiscale simulation methods in molecular science. Jülich Supercomputing Centre, 2009. ISBN 9783-9810843-8-2. 21. Kraus, J. K. Algebraic multilevel preconditioning of finite element matrices using local Schur complements. Numerical linear algebra with applications. 2006, 13, s. 49-70. DOI: 10.1002/nla.462. 22. Rivara, M.-C. Local modification of meshes for adaptive and/or multigrid finiteelement methods. Journal of computational and applied mathematics. 1991, 36, s. 79-89. DOI: 10.1016/0377-0427(91)90227-B. 23. Douglas, C. C.; Douglas, J. Jr. A unified convergence theory for abstract multigrid or multilevel algorithms, serial and parallel. Siam Journal of numerical analysis. 1993, 30, č. 1, s. 136-158. DOI: 10.1137/0730007. 24. Aricò, A.; Donatelli, M. A V-cycle multigrid for multilevel matrix algebras: proof of optimality. Numerical mathematics. 2007, č. 105, s. 511-547. DOI: 10.1007/s00211-006-0049-7. 25. Kraus, J. K. Algebraic multilevel preconditioning of finite element matrices using local Schur complements. Numerical linear algebra with applications. 2006, č. 13, s. 49- 70. DOI: 10.1002/nla.462. 26. Dendy, J. E. Jr.; Moulton, J. D. Black box multigrid with coarsening by a factor of three. Numerical linear algebra with applications. 2010, č. 17, s.577-598. DOI: 10.1002/nla.705 27. Xu, X.; Chen, W. Standard and economical cascadic multigrid methods for the mortar finite element methods. Numerical mathematics: Theory, methods and applications. 2009, 2, č. 2, s. 180-201. 28. Zhang, S. Multigrid algorithm for the coupling system of natural boundary element method and finite element method for unbounded domain problems. Journal of computational mathematics. 2007, 25, č. 1, s. 13-26. 29. Mayer, P.; Marek, I.; Pultarová, I. Convergence issues in the theory and practice of iterative aggregation/disaggregation methods. Electronic transactions on numerical analysis [online]. 2009, 35, s. 185-200. [cit. 4.8.2014]. ISSN: 51

1068-9613. dostupné na: etna.mcs.kent.edu/vol.35.2009/pp185-200.dir/pp185200.pdf,

nebo

na:

www.emis.de/journals/ETNA/vol.35.2009/pp185-

200.dir/pp185-200.pdf

Seznam obrázků Obrázek 1 Výsledek pro koeficienty 1,1,0,0. ................................................................. 36 Obrázek 2 Výsledek pro koeficienty 100,1,0,0. ............................................................. 36 Obrázek 3 Výsledek pro koeficienty exot....................................................................... 37 Obrázek 4 Výsledek pro multigrid po zpětněm vrácení do soustavy proměnných. ....... 39 Obrázek 5 Výsledek pro malý počet iterací multigridu. ................................................. 39 Obrázek 6 Výsledek různých multigridů pro koeficienty 1,1,0,0 pro 3 iterace. ............. 44 Obrázek 7 Výsledek různých multigridů pro koeficienty 100,1,0,0 pro 3 iterace. ......... 44 Obrázek 8 Výsledek pozměněných různých multigridů pro koeficienty 100,1,0,0 pro 3 iterace. ................................................................................................................... 44

52

Seznam použitých zkratek – reálná (komplexní) funkce reálné (komplexní) proměnné – integrovatelné;

𝑎𝑥𝑥

analogicky s vlnovkou 𝑣 – funkce z konkrétního prostoru funkcí 𝜆 – reálné (komplexní) číslo ∫Ω. – Integrace přes oblast A, Q – matice 𝑥 – vektor 𝑉 – vektorový lineární prostor 𝐴 – lineární zobrazení, jeho matice, operátor, konstanta v pasáži o variacích ‖. ‖𝑊𝑝𝑘(Ω) – norma na prostoru pod znaménkem indexu \ – množinové mínus 𝑇1 – lineární (kompaktní) operátor ℋ – Hilbertův prostor, tj. úplný prostor se skalárním součinem. 𝜀 – malé reálné číslo – blízko nule 𝜕

𝐷𝑖 , 𝜕𝑥 – parciální derivace (. , . ) – skalární součin 𝑊𝑝𝑘 (Ω) – Sobolevův prostor Γ1 , Γ – hranice oblasti 𝐺, Ω – oblast 𝐵1𝜇 , 𝐵 – diferenciální operátory na hranici oblasti 𝑑𝑆 – diferenciál přes hranici 𝑎(. , . ) – bilineární forma 𝑋 ′ – duální prostor 𝑅 + \{0} – kladná reálná čísla bez nuly (𝑛)

𝜆1

– vlastní číslo v pořadí aproximace dané horním indexem a jeho pořadí dle

velikosti dolním indexem L.-M– Lax – Milgramova věta 𝜉 – konstanta odhadující formu odspoda pro definice koercitivity 𝐾 – oblast pro definici konečného prvku 𝑃 – konečně dimenzionální prostor funkcí 53

ℐ𝐾 – interpolant funkce vzhledem ke konečnému prvku 𝐶0∞ (Ω), 𝐶(Ω) – prostory spojitých (dostatečně diferencovatelných) funkcí (s kompaktním supportem) ℰ𝑒 – lokální indikátor chyby 𝛼 – lokálně konstanta – reálné číslo, či prvek tělesa 𝑢𝑗 – výsledná funkce v pořadí aproximace 𝜘 – podmíněnost matice Diam – charakteristický rozměr oblasti/elementu 𝛾 – robustnost (chunkiness) 𝑙𝑜𝑔 – logaritmus 𝑁 – počet bazických funkcí konečného prvku 𝐴𝑘𝑟 , 𝐿𝑘𝑠 – indexované prvky matice 𝑥𝑘 , 𝑏, 𝑟𝑘 , 𝑝𝑘 – vektory (posloupnost vektorů) α𝑘 , 𝛽𝑘 – posloupnost čísel Det – determinant GMRES –generalized minimal rezidual method BiCG – metoda bikonjugovaných gradientů FOM – full orthogonalization method 𝐸 – jednotková matice ∑.𝑖=.. – suma přes index v mezích od dolní po horní ‖𝐴‖ – norma matice 𝑚𝑗 , 𝑆𝑖 , 𝐺𝑛𝑐 – indexované množiny 𝑑𝑒𝑟 – označuje možnost doplnit za tento symbol libovolnou parciální derivaci prvního stupně (dáno typem úlohy) 𝑥 𝑇 – transpozice vektoru nebo matice j – jiný postup nebo podmínky počet iterací – hrubší(hrubé(jemné)) iterace vícestupňových metod … – jiný pokus (v tabulce výsledků běhů programů, tyto i dále) h.– heap (tj. halda) s.– stack st. – static hi – heap kombinováno s v kódu dalšími úpravami pam. – paměti 54

Přílohy 1 Obsah přiloženého DVD Na přiloženém DVD se v kořenovém adresáři nachází tato bakalářská práce ve formátu bakalarska_prace.pdf s jednoduchým návodem navod.txt pro obsluhu programů. Struktura projektu je tvořena: vlc_ge_heapinl, rkoh_kongposym, rkoh_mgridnot, rkoh_notgr_sec, rkoh_owngrid, rkoh_real_choleski, rkoh_ridke_psymkonjgr, a program pro sledování časovače během běhu programu.

55