Nelineární systémy
Fázové portréty Hezké příklady nelineárních systémů
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 2
Numerická konstrukce fázových portrétů
Pro numerické řešení obyčejných diferenciálních rovnic existuje mnoho programů Můžeme je použít ke kreslení fázových portrétů Několik praktických rad, jak nakreslit hezké portréty (další rady v literatuře)
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 3
Výpočet trajektorie
Trajektorii (řešení) rovnice , která prochází bodem najdeme tak, že 1. řešíme rovnici z bodu kupředu,
tj. v rostoucím (kladném) čase 2. řešíme rovnici z bodu dozadu,
tj. v klesajícím (záporném) čase a to tak, že v kladném (normálním) čase řešíme rovnici
Platí totiž
( )x f x=
( )x f x= 0x
0(0)x x=
0( ), (0)x f x x x−= =
( )x f x= 0(0)x x=
( ) ( ) ( )( ( )) ( ( )) ( ( ))( ) ( )
dx t dx t dx tx f x t f x t f x tdt d t d t
− −= = → = − → = − −
−
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 4
fwd=dsolve('Dx=exp(x)','x(0)=-2'), ezplot(fwd), hold on fwd = -log(-t+exp(2)) bwd=dsolve('Dx=-exp(x)','x(0)=-2'),ezplot(bwd) bwd = -log(t+exp(2))
Příklad
Rovnice řešíme dopředu a dozadu
2( ) ln, (0) 2 ( )x xx e x t t e= − −− → += =2( ) ln( ), (0) 2x xx te tx e= − = − → = − +
, (0) 2xx e x= = −
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 5
Použití pro fázový portrét
Fázový portrét systému
kreslíme takto: 1. vyberme jeden bod trajektorie 2. část trajektorie, vycházející z tohoto bodu („dopředu“)
vypočteme integrací rovnic s počátečními podmínkami
3. část trajektorie, končící v tomto bodě („dozadu“) vypočteme integrací rovnic s počátečními podmínkami
1 1 1 2
2 2 1 2
( , )( , )
x f x xx f x x==
1 2[ , ]P Px x
1 1 1 2
2 2 1 2
( , )( , )
x f x xx f x x==
1 1 2 2(0) , (0)P Px x x x= =
1 1 1 2
2 2 1 2
( , )( , )
x f x xx f x x= −= −
1 1 2 2(0) , (0)P Px x x x= =
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 6
Další důležité věci
Uvedený postup je jedinou možností, jak dostat hezký portrét v okolí nestabilních uzlů, nestabilních ohnisek nestabilních cyklů
Důležitou věcí je také správně zvolit rozsahy: aby portrét obsahoval všechny zajímavé jevy všechna ekvilibria a aby v něm všude integrační metoda fungovala když toho dopředu moc nevíme, postupujeme iteračně
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 7
Přehled
Tlumené kyvadlo – podle Newtona Oscilátor – podle Van der Pola Nosník – podle Eulera Klarinet a housle – podle Reyleigho Dravci a oběti – podle Volterry Tunelová dioda Vítr DW oscilátor Adaptivně řízený systém další příklady
Příklady nelineárních systémů
Tlumené kyvadlo – podle Newtona
Isaac Newton Rovnice Simulink
Fázový portrét
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 9
Sir Isaac Newton 1643 – 1727
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 10
Isaac Newton
Newton [njútn], sir Isaac, * 4. 1. 1643, † 31. 3. 1727, anglický fyzik, matematik a astronom; člen, v letech 1703 – 27 předseda Královské spol. v Londýně. Roku 1705 povýšen do šlechtického stavu. Zakladatel klasické mechaniky. Objevil gravitační zákon, podílel se (souběžně s G. W. Leibnizem) na vytvoření infinitezimálního počtu. Prováděl optické výzkumy, objasnil rozklad světla, zkonstruoval zrcadlový dalekohled. Zabýval se též alchymií.
V díle Philosophiae naturalis principia mathematica (Matematické základy přírodovědy) formuloval tři základní zákony dynamiky (dnes nazývány Newtonovy).
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 11
Tlumené kyvadlo - rovnice
Pohybová rovnice v tečném směru Pro stavové proměnné dostaneme model
sinml mg klϕ ϕ ϕ= − −
1 2,x xϕ ϕ= =
1 2
2 2 1sin( )
x xk gx x xm l
=
= − −
Podobné rovnice:
• synchronní generátor připojený na nekonečné vedení
• Josephsonovy obvody
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 12
Tlumené kyvadlo - Simulink
1 2, 0π= =x x
cokoli jiného
Model v Simulinku pend_damp
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 13
Tlumené kyvadlo - fázový portrét
fázový portrét tlumeného kyvadla
demoph2
Oscilátor – podle van der Pola
Balthasar van der Pol Oscilátor
Fázový portrét
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 15
Balthasar van der Pol
Balthasar van der Pol, 1889-195918 Holandský elektroinženýr. Jako jeden z prvních studoval experimentálně dynamiku v laboratoři v dvacátých a třicátých letech. Zkoumal elektrické obvody s vakuovými lampami a objevil, že mohou mít stabilní oscilace, dnes nazývané limitní cykly. V rovce 1927 publikoval (s kolegou van der Markem) v Nature článek o „nepravidelném šumu“, který pozoroval pro některé budící frekvence. Z rekonstrukce jeho pokusí dnes víme, že objevil deterministický chaos. Sestrojil mnoho elektronkových modelů lidského srdce a na nich studoval jeho dynamiku. Zkoumal vliv buzení analogický vlivu kardiostimulátoru. Snažil se najít způsob, jak stabilizovat srdeční arytmie.
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 16
Oscilátor
RLC obvod s lineárními L, C a nelineárním rezistorem s kubickou charakteristikou
Pro stavy jsou stavové rovnice
2( 1)i v vα= −
1 2,L Cx i x v= =
( )1 2
22 1 2 2
1
1 ( 1)
x xL
x x x xC
α
=
= − + −
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 17
Oscilátor - fázový portrét
Fázový portrét pro triviální řešení periodické řešení
1L C= =
nestabilní ekvilibrium
limitní cyklus
demophVanDerPol
Nosník – podle Eulera
Leonard Euler Rovnice
Netlumený nosník Malé zatížení Větší zatížení
Ekvilbria Bifurkace
Tlumený nosník
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 19
Leonard Euler 1707 – 1783
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 20
Leonard Euler
Euler [ojlr] Leonhard * 15. 4. 1707, † 18. 9. 1783, švýcarský matematik, fyzik a astronom. Působil v Petrohradě a v Berlíně. Napsal řadu učebnic matematiky. Pracoval v teorii čísel a integrálů; zabýval se diferenciálními rovnicemi, variačním počtem, exponenciálními a goniometrickými funkcemi i úlohami geometrickou tematikou.
Známá je např. Eulerova hypotéza (později přesně dokázaná), že tzv. úloha o 36 důstojnících nemá řešení. Pro geofyziku má význam zejména jeho studie rotace tuhého tělesa. Viz též Eulerova perioda.
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 21
Ohnutý nosník – rovnice
Pohybová rovnice ve směru x Pro stavové proměnné dostaneme model
3 0mx dx x x xµ λ+ − + + =
1 2,x x x x= =
1 2
32 2 1 1
1x x
dx x x xm m m
µ λ=
−= − + −
tlumení zatížení pružná síla v nosníku
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 22
Nosník netlumený - malé zatížení
Pro malé zatížení se nosník stlačí, ale neprohne
µ λ<demophBuckledBeam1eq
Hamiltonovský systém
0d =
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 23
Nosník netlumený - větší zatížení
Pro velké zatížení se nosník prohne
µ λ>demophBuckledBeam2eq
0d =
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 24
Ohnutý nosník - ekvilibria
23
2 1 1
0
0 ( )µ λ
=
= − + − −
xdx x x
2 12
2 1
0, 0
0,
x xx x µ λ
= =
= = −
1 2
32 2 1 1
1x x
dx x x xm m m
µ λ=
−= − + −
1 2
1 2
: 0
: 0, ; 0
x x
x x
µ λ
µ λ µ λ
≤ = =
> = ± − =
1 ekvilibrium
3 ekvilibria
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 25
Nosník - Bifurkace
µ λ>
µ λ<
µ λ=
Přechod mezi neohnutým a ohnutými stavy
Bifurkace
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 26
Nosník - tlumený
není Hamiltonovský
0.310.8
dλµ
===
demophBuckledBeam1eqDumped
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 27
Nosník – tlumený, větší zatížení
0.313
dλµ
===
demophBuckledBeam2eqDumped
Klarinet a housle – podle Rayleigho
Baron Rayleigh Zvuk hudebních nástrojů
Foukání na klarinet Vliv parametrů
Smyčec a housle Fázový portrét
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 29 John William Rayleigh, 1842 – 1919
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 30
Baron Rayleigh
Rayleigh [rejli] John William, baron, 12. 11. 1842, † 30. 6. 1919, britský fyzik; profesor na univerzitě v Cambridgi. 1879 – 84 ředitel Cavendishovy laboratoře, člen a 1905 – 08 ředitel Královské společnosti v Londýně. Zabýval se optikou, akustikou a elektromagnetismem, objevil (spolu s W. Ramsayem ) argon a zákon o záření absolutně černého tělesa. Rayleighovy výzkumy přispěly k rozšíření fyzikálních znalostí o stavu hmoty. Nobelova cena v roce 1904 za výzkumy o hustotě nejdůležitějších
plynů a za objev argonu.
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 31
Teorie zvuku hudebních nástrojů
J. W. Rayleigh: The Theory of Sound: Vols. I and II. Dover (1945 edition), 1887.
Rayleigh rozlišuje dva druhy hudebních nástrojů:
„perkusní (bicí)“: bubny, kytary, piána modeluje tlumenými oscilacemi, jednoduchá dynamika vlastně jen přechod do ustáleného stavu „udržované“: smyčcové a dechové moduluje „udržovanými oscilacemi“ = uzavřenými orbity
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 32
Foukání na klarinet – podle Rayleigho
Rayleigh popisuje jazýček klarinetu jako lineární oscilátor a efekt klarinetistova foukání členem na pravé straně (záporné tlumení pro malé a kladné pro velké). Tedy celkem nebo
0x kx+ =
3( ) , , 0x xα β α β− >
x
3( ) 0x x x kxα β− + + =
1 23
2 2 2 1
x xx x x kxα β
=
= − −
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 33
Klarinet - fázový portrét
Fázový portrét demophClarinet
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 34
Klarinet – Vliv parametrů na zvuk
1kα β= = =
tužší jazýček (tužší pružina, větší k)
bohatší tón , 1, 3kα β = =
tvrdší foukání (širší charakteristika tření, menší )
hlasitější zvuk , 1, 0.5kα β= =
β
demophClarinet,demophClarinet2,demophClarinet3
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 35
Rayleigh navrhl model analogický systému hmotnost-pružina položenému na pásovém dopravníku s konstantní rychlostí Tedy celkem
nebo stavově
Smyčec a houslová struna
( ) 0mx kx f x b+ + − =
1 2
2 1 21 ( )
x xkx x f x bm m
=
= − − −
struna bez smyčce ~ hmotnost-pružina
použití smyčce ~ tření mezi pásem a tělesem
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 36
Smyčec a struna - fázový portrét
limitní cyklus - asi špatně kreslí ?
demophBowingString
Dravci a oběti – podle Volterry
Vito Volterra Dravec a oběť Omezený růst
Soutěžící populace
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 38
Vito Volterra 1860-1940
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 39
Vito Volterra
Vito Volterra, 1860 - 1940
Volterra publikoval články o parciálních diferenciál-ních rovnicích, zejména o rovnicích cylindrických vln. Nejslavnější je jeho práce o integrálních rovnicích, zejména o té, které se dnes říká „Volterrova integrální rovnice“. Další info na:
www-history.mcs.st-andrews.ac.uk/history/ Mathematicians/Volterra.html
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 40
Dravec a oběť
Hrabě Vito Volterra: Model vývoje dvou populací v Jaderském moři: oběť x dravec y Model je velmi zjednodušený. Možno přidat např.
omezený růst (o. bez d. a d. s o.) kvůli sociálním faktorům soutěžící druhy (každý druh jí ten druhý)
x ax bxyy cxy dy= −= −
, , , , , 0a b c d x y >
( )( )
x a by x xy cx d y y
λµ
= − −= − −
, , , , , 0a b c d λ µ >
x ax bxyy dy cxy= −= −
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 41
Dravec a oběť
další kvadranty nemají fyzikální smysl
(0,0);(1,1);ex =
oběti
drav
ci
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 42
Dravec a oběť s omezeným růstem
(0,0);(1,0); (0, 1)ex = −
drav
ci
oběti
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 43
Soutěžící populace
(0,0);(1,1);ex =
druh1
druh
2
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 44
Obvod s tunelovou diodou
RLC obvod s tunelovou diodou
cC
LL
dvi Cdtdiv Ldt
=
=
( )R Ri h v=
( )
( )
1 1 2
2 1 2
1 ( )
1
x h x xC
x x Rx uL
= − +
= − + +
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 45
Charakteristika tunelové diody
Polynomiální charakteristika ( )R Ri h v=
2 3 4 5( ) 17.76 103.79 229.62 226.31 83.72R R R R R Rh v v v v v v= − + − +
až 3 rovnovážné stavy podle zátěže
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 46
Fázový portrét
x=0:.01:1;plot(x,tunneldiode(x))
demophTunnel
[ ]Rv V
[ ]Ri mA
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 47
Počet rovnovážných stavů
Vítr – podle
Oscilace vynucené větrem
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 49
Systém s oscilacemi vynucenými větrem
Sastry, s. 342
tlumení rozladění
1 1 1 2 2 1 2
2 22 2 1 1 2 1 2
1 ( )2
x x x x x
x x x x x
µ µ
µ µ
= − − +
= − + +
1
2
01
µµ
==
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 50
Vítr
případ limitní cykly (heteroklinické orbity)
(-2,0),(0, 0)ex =
1
2
01
µµ
==
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 51
Vítr
případ
(-2.383,0.704),(0, 0)ex =
1 2 1µ µ= =
Double Well oscilátor
Rovnice Vliv buzení
Chaos Citlivost na počáteční podmínky
Poincarého řez Bifurkace
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 53
Double Well oscilátor
magneto-elastický mechanický systém v klidu má 2 ustálené stavy periodická budící síla popsán Duffingovou rovnicí tlumení amplituda budící síly budící frekvence
3 cos( )x bx x x F tω+ − + =
bFω
1 23
2 2 1 1 cos( )
x xx bx x x F tω=
= − + − +
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 54
DW Oscilátor - bez buzení
Bez buzeni 0.05
01
bFω
===
0, 1ex = ±
double well = dvojitá jáma
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 55
DW Oscilátor
Bez buzeni
1 20.05, 0, 1 ( (0) 1.5, (0) 1)b F x xω= = = = = 0, 1ex = ±
DWoscilatorMS
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 56
DW Oscilátor - malé buzení
Malé buzeni
2 stabilní limitní cykly
0.250.221
bFω
===
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 57
DW Oscilátor
Malé buzeni
1 20.25, 0.22, 1 ( (0) 1, (0) 0)b F x xω= = = = =
DWoscilatorMS
Stabilní limitní cyklus
je tam ještě jeden
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 58
DW Oscilátor - větší buzení
Větší buzeni – přechodný chaos
jako dříve 2 stabilní limitní cykly ale přechodový jev je chaotický při malé změně počát. podmínek může skončit na druhé straně hranice oblastí počátečních podmínek vedoucích k jednotlivým cyklům jsou složité: mají fraktální charakter
0.250.2451
bFω
===
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 59
DW Oscilátor
Větší buzeni
1 20.25, 0.245, 1 ( (0) 1, (0) 0)b F x xω= = = ≅ =
1(0) 1x =
1(0) 1.01x =
1(0) 0.985x =
1(0) 1.02x =
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 60
DW Oscilátor
Ještě větší buzeni - chaos
chaotické chování
0.250.41
bFω
===
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 61
DW Oscilátor
Ještě větší buzeni
1 20.25, 0.4, 1 ( (0) 1, (0) 0)b F x xω= = = = =
Chaos
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 62
DW
Při dalším zvětšování buzení se chování zase mění Např. pro F = 0.5 se zdá být zase limitní cyklus obdobně při změnách dalších parametrů zkuste sami 0.15, 0.3, 1b F ω= = =
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 63
DW Citlivost na počáteční podmínky
Co to znamená velká citlivost na počáteční podmínky?
0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 a 4.0
Jak se rozpliznou blízké počáteční stavy po
cyklech délky 2*pi
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 64
DW Citlivost na počáteční podmínky
A takhle to vypadá po 25 cyklech: Poznáte, kde to začalo?
Poznámka: Výpočet těchto 25
cyklů při 200 krocích integrace v cyklu pro 40000
počátečních podmínek trval přes 7 hodin na rychlém počítači
s 9000 mips! www.apmaths.uwo.ca/ ~bfraser/
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 65
DW
další zajímavý obrázek: extrém vs. extrém
extrém v kroku i
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 66
DW – Poincarého řezy
Přidáme-li jako třetí rozměr velikost budící síly, dostaneme třírozměrný fázový portrét. Jeho dvourozměrný řez (třeba pro sílu ) odhalí, že to je (3-D) podivný atraktor. je překvapivé že po 2pi skáče po tak hezké křivce (srovnej obr
cos(2 )f F n Fπ= =
Zimní semestr 2006 SRI | M. Šebek | FEL ČVUT 67
I zde jsou bifurkace
1, [1.0,1.1], 1b F ω= ∈ =
animace bifurkací