Lorenzův atraktor. MM semestrální práce. Jméno a příjmení: Pavel Martínek Osobní číslo: A08N0203P. Datum odevzdání: 12.2

Fakulta aplikovaných věd

Pavel Martínek MM – semestrální práce

Lorenzův atraktor

Jméno a příjmení: Osobní číslo: Obor: E-mail: Datum odevzdání:

Pavel Martínek A08N0203P MA [email protected] 12.2.2009

Strana 1 (celkem 25)



Obsah Lorenzův atraktor........................................................................................................................1 Úvod............................................................................................................................................3 Dynamický systém.................................................................................................................3 Lineární systém......................................................................................................................3 Stavový prostor.......................................................................................................................3 Atraktor dynamického systému..............................................................................................4 Chaos......................................................................................................................................4 Model konvektivního proudění atmosféry .................................................................................8 Numerické experimenty v programu Mathematica...................................................................12 Stacionární body...................................................................................................................16 Přesnost při výpočtu v programu Mathematica....................................................................22 Závěr.........................................................................................................................................26 Zdroje........................................................................................................................................27




Úvod

Dynamický systém Dynamický systém sestává ze stavového prostoru, jehož souřadnice popisují stav systému v daném čase a z dynamických podmínek, které popisují změnu tohoto systému v čase. Stav systému je potom popsán vektorem, který celý leží ve stavovém prostoru. Dynamické podmínky jsou většinou zadány soustavou diferenciálních rovnic, které popisují změnu stavového vektoru v čase. Změna stavu dynamického systému se děje provedením těchto diferenciálních rovnic a nahrazením starého stavového vektoru vektorem novým. Dynamický systém může být deterministický nebo stochastický (náhodný). Deterministický dynamický systém lze poměrně přesně popsat, zatímco u systému stochastického jsme odkázáni pouze na statistické vlastnosti takového systému (například střední hodnota, disperze, směrodatná odchylka, centrální moment a jiné).

Lineární systém Lineární systém je takový systém, v němž lze uplatnit princip superpozice. Superpozice využíváme při řešení velkého množství problémů, například při řešení průtoku elektrického proudu v elektronických obvodech nebo ve fyzice při skládání působení sil na hmotný bod. Obecně platí, že je-li systém lineární a lze využít superpozice, je řešení takového systému často velmi jednoduché a jednoznačné. Chování takových systémů lze předpovědět i do budoucnosti. Při práci se systémy, které nejsou lineární, častou používáme postupu zvaného linearizace. Nelineární závislost se nahradí závislostí lineární.

Stavový prostor Stavový prostor určuje, jakých hodnot může nabývat stavový vektor dynamického systému. Stavový vektor je tvořen množinou proměnných, které mohou nabývat hodnot z určitého intervalu. Interval všech těchto hodnot potom určuje celý stavový prostor. Stavový prostor může být několika typů: • konečný • spočitatelný • nekonečný




Atraktor dynamického systému Atraktor (anglicky attractor) dynamického systému je stav, do kterého systém směřuje. Je to tedy množina, ve které je stavový vektor, když je systém v nekonečném čase. Atraktory rozdělujeme do několika tříd: • • • • •

atraktorem jsou pevné body atraktorem jsou periodické body atraktorem jsou kvaziperiodické body atraktor je chaotický podivný atraktor

Jsou-li atraktorem dynamického systému pevné body, jde o nejjednodušší případ. Systém se tedy v nekonečném čase ustálil v nějakém stabilním stavu a v podstatě už nejde o dynamický systém. Příkladem může být kyvadlo, které se vlivem odporu vzduchu a odporu ložisek zastaví v nejnižším bodě své dráhy. Jsou-li atraktorem periodické (resp. kvaziperiodické) body, jde také o jednoduchý případ. Systém se ustálil tak, že osciluje mezi několika stavy. Příkladem je těleso, které se na své cestě vesmírem dostane do blízkosti velmi hmotného tělesa. Po určitém čase se pohyb tohoto tělesa ustálí na eliptické dráze. Je-li atraktor chaotický, znamená to, že výsledný atraktor nelze v podstatě nijak dopředu předpovědět. To je způsobeno tím, že je systém velmi citlivý na počáteční podmínky. Chaotičnost v tomto případě neznamená náhodnost, protože se bavíme o deterministických systémech. Příkladem může být koule postavená na vrcholku jehlanu. Jakýkoliv vnější podnět způsobí, že koule tento stav opustí a dostane se do některého atraktoru (místo pod jehlanem). Tento atraktor nelze předpovědět, protože nemůžeme bez zásahu do měření zjistit počáteční podmínky. V kvantové fyzice existuje takzvaný princip neurčitosti, který má obdobný význam pro kvantové jevy. Podivný atraktor (anglicky strange attractor) je nejzajímavějším případem atraktoru. Tento typ atraktoru vzniká, je-li systém popsán minimálně třemi diferenciálními rovnicemi. Takový systém může mít velmi komplikovaný atraktor, který sice bude chaotický, ale přesto bude vykazovat určité pravidelnosti. Termín podivný atraktor není ještě přesně matematicky definován, ale považujeme za něj takový atraktor, který vykazuje stejné vlastnosti, jaké mají fraktály (podivný atraktor je tedy fraktálem).

Chaos Chaos je z časového hlediska budoucí stav deterministického dynamického systému, který není předpovídatelný v důsledku velké citlivosti na počáteční podmínky. Chaos může nastat v systému, který má více než dvě stavové proměnné, tedy například v trojrozměrném prostoru. Pro diskrétní procesy existuje i jednorozměrný chaos. Jako příklad můžeme uvést takzvanou Logistickou mapu, která je definována předpisem: f(x) = r*x(1-x) Tato funkce je pro r > 4 chaotická.




Obr. 1: Celkový pohled na Mandelbrotovu množinu

Obr. 2: Detail z Mandelbrotovy množiny




Obr. 3: Dynamický systém - Eulerova substituce

Obr. 4: Dynamický systém - simulace difúze




Obr. 5: Dynamický systém - simulace difúze




Model konvektivního proudění atmosféry Proudící tekutina je poměrně složitý spojitý dynamický systém, s velkou variabilitou okrajových podmínek, lze jej charakterizovat velmi vysokým počtem stupňů volnosti. Studium vlastností takového systému je obecně technicky velmi obtížné. Proto budeme demonstrovat chaotické chování dynamických systémů na případě daleko jednoduššího systému s nízkým počtem stupňů volnosti. Základní mechanismus vzniku chaosu je společný všem dynamickým systémům bez ohledu na jejich složitost. Jako příklad jsme zvolili Lorenzův systém, který vlastně odstartoval éru systematického studia chaosu.

I když by se na první pohled mohlo zdát, že x, y, z jsou prostorové souřadnice, ale není tomu tak. Jejich fyzikální význam je poněkud abstraktní. Proměnná x představuje rychlost rotace pohybu částice, kladná hodnota je ve směru hodinových ručiček. Proměnná y je rozdíl teplot stoupající a klesající tekutiny. Proměnná z charakterizuje odchylku svislého profilu teploty od lineárního průběhu. Parametr r je Rayleigho číslo (v normovaném tvaru), sigma Prandtlovo číslo a parametr b představuje štíhlost válce tekutiny při konvekci, tedy poměr jeho délky a průměru. Tento matematický model zachycuje základní vlastnosti konvektivního proudění atmosféry, která je zahřívána povrchem ze spodu a ochlazována svrchu. Vzniká tak rotační pohyb částic vzduchu, kdy ohřátá částice stoupá, tím se ochlazuje a začne klesat, aby se opět zahřála a stoupala. Tento jev je známý jako Rayleigh-Bénárdova nestabilita. Okrajové podmínky jsou poněkud idealizovány: proudění v horní oblasti je považováno bez smykového napětí místo realističtější podmínky stejných rychlostí, v příčném směru je uvažována periodická okrajová podmínka místo omezení stěnami a celý případ je modelován jako rovinný místo prostorového. Lorenzovy rovnice mají následující vlastnosti: • Jsou autonomní, to znamená, že jejich pravá strana explicitně neobsahuje čas, koeficienty jsou konstantní; • Obsahují pouze první časové derivace. Důsledkem tohoto spolu s uvážením autonomie systému je zřejmé, že jeho vývoj závisí pouze na okamžitých vlastnostech proměnných (x, y, z) a nikoli na jejich historii; • rovnice jsou nelineární, viz členy xz a xy ve druhé a třetí rovnici; • řešení soustavy rovnic je omezené v prostoru proměnných.




Atraktor příslušející Lorenzovu systému je prvním z tzv. "podivných atraktorů" charakterizujících chaotické chování dynamického systému a má některé vskutku podivné vlastnosti: • je tvořen spojitou křivkou v prostoru, která obecně začíná v jistém počátečním bodě, může však mít nekonečně velkou délku. Přitom vyplňuje jistý přesně vymezený podprostor ve fázovém prostoru, ze kterého nikdy nevybíhá; • nikdy neprotíná sám sebe, nekříží se ani se neopakuje; • má vlastnost fraktálů, tj. jeho struktura se opakuje na různých měřítkách; • jeho průběh v prostoru je náhodný, chaotický, nepředpověditelný. Ukazuje se, že kritická hodnota parametru r při výše uvedených hodnotách parametrů s a b je rovna asi 24,74, pro hodnoty nižší směřuje vývoj systému do jediného bodu ve fázovém prostoru, pro hodnoty vyšší dostáváme nekonečný pohyb s prvky chaosu.

Obr. 7: Rozložení teploty pro nízké r (nižší než kritická hodnota)

Obr. 6: Rozložení teploty pro vysoké r (vyšší než kritická hodnota)




Obr. 8: Schéma Lorenzova systému – Reygleigh-Bénárdova buňka




Obr. 9: Grafické znázornění proudění pro různé r (vyšší než kritická hodnota)




Numerické experimenty v programu Mathematica Při numerických experimentech jsme vycházeli z následujícího tvaru Lorenzových rovnic. dx =−b x t y t  z t dt dy =−s y ts z t dt dz =− y t x tr y t−z t dt Při zkoumání tohoto modelu se nejčastěji používá následujících hodnot pro parametry b = 8/3, s = 10, r = 28. yHt L 20

10

10

20

30

40

-10

Obr. 10: Fázový diagram - {x, y} pro y0 = [0,3,10]


50

x Ht L



zHt L 30 20 10 10

20

30

40

50

xHt L

-10 -20 Obr. 11: Fázový diagram - {x, z} pro y0 = [0,3,10]

zHt L 30 20 10 -10

10 -10 -20

Obr. 12: Fázový diagram - {y, z} pro y0 = [0,3,10]


20

yHtL



8xHt L, yHt L, zHtL<

40

30

20

10

t 10

20

30

40

50

60

-10 -20

Obr. 13: Všechny složky řešení pro y0 = [0,3,10]


70



Stacionární body a) S1 = [0, 0, 0] det(I – J(S1)) = - (b + ) * (s – r * s +  + s *  + 2 ) 1 = -b, 2,3 = ½ (-1 – s  (1 – 2s + 4r * s + s2)½). 2 = -½((1 + s) - ((1 + s)2 + 4s * (r - 1))½) > 0 => S1 je nestabilní. b) Z druhé rovnice vyplývá, že z = y. Nyní třetí rovnici vydělíme y (pro y = 0 bychom dostaly S1) a vyjde nám že x = r - 1. Po dosazení za x do první rovnice nám vyjde, že y = z = n => S2,3 = [r – 1, n, n], kde n = [b(r-1)]½. det(I – J(S2)) = 3 + (b + s + 1) * 2 + (b + b * s + n2) *  + 2n2 * s Nyní dosadíme za b, r, s a vyjdou nám komplexně sdružené kořeny s kladnou reálnou částí => S2,3 je nestabilní.




Pro počáteční hodnotu rovnu jednomu ze stacionárních bodů S2,3 je výsledkem opět tento bod. Pro obvyklou volbu b, r, s je y0 = [27,  6 * 2½,  6 * 2½]. a) Vychýlení počáteční hodnoty y0 pro x o 5 * 10 -5 Výsledné řešení se pro nízká t příliš neliší od konstantního: y HtL 8.488 8.487 8.486 8.485 8.484 8.483 26.994

26.996

26.998

27.002

27.004

27.006

x Ht L

Obr. 14: Fázový diagram – {x, y}

zHtL 26.994

26.996

26.998

27.002 8.488

8.486

8.484

8.482

Obr. 15: Fázový diagram - {x, z}


27.004

27.006

xHtL



zHtL 8.483

8.484

8.485

8.486

8.487

yHtL

8.488

8.488

8.486

8.484

8.482

Obr. 16: Fázový diagram - {y, z}

8 yHtL - zHt L<

0.003 0.002

0.001 t 10

20

30

40

50

-0.001 -0.002 -0.003

Obr. 17: Rozdíl složek řešení {y, z}


60

70



b) Chování systému pro r = 24, y0 a t  (0, 200). y Ht L

9

8

7

22

24

26

xHtL

Obr. 18: Fázový diagram – {x, y}

zHtL

10 9 8 7 6

22

24

Obr. 19: Fázový diagram – {x, z}


26

xHtL



zHt L

10 9 8 7 6

7

8

9

yHtL

Obr. 20: Fázový diagram - {y, z}

8 yHt L -zHt L<

2

1

t 50

100

-1


150

200



Přesnost při výpočtu v programu Mathematica Pro výpočet relativní chyby jsem využil vzorce: er(t) = Max[xok(t) – xn(t), yok(t) – yn(t), zok(t) – zn(t)] * (Max[xok(t), yok(t), zok(t)])-1 xok – x složka řešení, které bylo numericky spočítáno na min 100 platných cifer xn – x složka řešení, které bylo spočtená danou metodou (y0 = [0,3,10]). t  {1, 2, …, 20} NDSolve – funkce programu Mathematica pro řešení dif. rovnic. V tomto případě bez jakýchkoliv nastavení.

Relativní c hyba 3.5  10 6 3  10 6 2.5  10 6 2  10 6 1.5  10 6 1  10 6 5  10 7 5

10

15

20

Čas

Pro klasický NDSolve

Relativní c hyba

Parametr WorkingPrecision -> 32 znamená, že v průběhu výpočtu bude program pracovat s čísly o 32 číslicích. Dále již jen 32.

8  10 8 6  10 8 4  10 8 2  10 8

5

10

15

20

NDSolve s WorkingPrecision -> 32


Čas



Relativní c hyba

NDSolve[LorenzEquations, {x,y,z},{t,0,20}, WorkingPrecision->32,Metod-> StiffnessSwitching ]. Pro vyplněný parametr metod počítá mathematica pouze pomocí této (těchto) numerických metod. Dále vždy uvádím jen jméno použité metody (v anglickém originálu, jméno je shodné s probíranými metodami).

0.001 0.0008 0.0006 0.0004 0.0002

5

10

15

20

Čas

Mé nastavení

Relativní c hyba 8  10 8 6  10 8 4  10 8 2  10 8

5

10

15

20

15

20

Čas

Adams + 32 Relativní c hyba 3.5  10 6 3  10 6 2.5  10 6 2  10 6 1.5  10 6 1  10 6 5  10 7 5

10

Čas

Adams




Relativní c hyba 0.000025 0.00002 0.000015 0.00001 5  10 6 5

10

15

20

Čas

BDF Relativní c hyba 8  10 8 6  10 8 4  10 8 2  10 8

5

10

15

20

5

10

15

20

Čas

BDF + 32

Relativní c hyba 8  10 8 6  10 8 4  10 8 2  10 8 Čas

ExplicitRungeKutta + 32




Relativní c hyba

4  10 6 3  10 6 2  10 6 1  10 6

5

10

15

20

5

10

15

20

5

10

15

20

Čas

ExplicitRungeKutta Relativní c hyba 8  10 7 6  10 7 4  10 7 2  10 7 Čas

ImplicitRungeKutta Relativní c hyba 7  10 7 6  10 7 5  10 7 4  10 7 3  10 7 2  10 7 1  10 7 Čas

ImplicitRungeKutta + 32




Závěr V této práci jsme popsali Lorenzův systém, který popisuje cirkulaci tepla v zemské atmosféře. A i když jde o model značně zjednodušující skutečnost, podařilo se nám ilustrovat jeho chování pro různé volby parametru r a různé počáteční podmínky. Pokud bychom chtěli předpovídat počasí na základě tohoto jednoduchého modelu, museli bychom se vypořádat s jeho nepříjemným chaotickým chování pro r vetší než je jeho kritická hodnota. Přestože Lorenz formuloval tento model už v roce 1963 a dnešní meteorologie již jistě používá mnohem sofistikovanější modely, zcela určitě budou vykazovat podobné chování jako rovnice námi zkoumané. Proto věřím, že po shlédnutí příští relace o počasí se spokojíte s přáním slunce v duši a nevydá se všanc stochastickému fenoménu jakým je počasí.




Zdroje 1. Pavel Martínek, seminární práce ze SNM1 2. Ing. Václav Uruba, Csc, Náhoda v exaktní vědě, http://www.essentia.cz/ 3. Pavel Tišnovský, Fraktály, http://www.vood.mysteria.cz/fraktaly


Lorenzův atraktor. MM semestrální práce. Jméno a příjmení: Pavel Martínek Osobní číslo: A08N0203P. Datum odevzdání: 12.2

Recommend Documents