´ NI´ VE FYZICE MATEMATICKE´ MODELOVA
Prˇ´ırodoveˇdecka´ fakulta UP Olomouc za´ˇr´ı 2003
Matematicke´ modelova´nı´ ve fyzice
OBSAH
Obsah 1
Za´kladnı´ pojmy
2
2
Matematicky´ model 2.1 Stavova´ rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3 3
2.2 2.3
8 9
3
4
Rovnova´zˇny´ stav dynamicke´ho syste´mu
Rovnice statiky dynamicke´ho syste´mu . . . . . . . . . . . . . . . . . . . . . . . . . 9 Rovnova´zˇny´ stav linea´rnı´ho syste´mu . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3 3.4 3.5
Astatismus dynamicke´ho syste´mu . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Singula´rnı´ body statvove´ rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Singula´rnı´ bod typu uzel a ohnisko . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Numericke´ modelova´nı´ 12 4.1 Diskretizace spojite´ho syste´mu numerickou metodou ˇresˇenı´ . . . . . . . . . . . . . . 12 Chyby numericky´ch vy´pocˇtu˚ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Oveˇrˇenı´ funkcˇnosti modelu 14 5.1 Optimalizace parametru˚ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 5.2
6
9
3.1 3.2
4.2 5
Deduktivnı´ a induktivnı´ formulace modelu . . . . . . . . . . . . . . . . . . . . . . . Linearita ve stavove´ formulaci . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Validace modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Seznam pouzˇite´ literatury
16
1
1 ZA´KLADNI´ POJMY
Matematicke´ modelova´nı´ ve fyzice
1 Za´kladnı´ pojmy Syste´m Pojmem syste´m rozumı´me jake´koli u´cˇelove´ usporˇa´da´nı´ jednodusˇsˇ´ıch objektu˚, resp. komponentu˚ ve slozˇiteˇjsˇ´ı celek charakterizovany´ soucˇinnostı´, resp. interakcı´ teˇchto objektu˚ a vykazujı´cı´ urcˇite´ vy´sledne´ vlastnosti odpovı´dajı´cı´ u´cˇelu, k neˇmuzˇ bylo toto usporˇa´da´nı´ vytvorˇeno. Syste´m jako model V pu˚vodnı´m slova smyslu jsou jednodusˇsˇ´ımi, dı´lcˇ´ımi objekty fyzicky existujı´cı´ a oddeˇlitelne´ komponenty a jejich soucˇinnost spocˇ´ıva´ v jejich vza´jemne´m pu˚sobenı´ naprˇ. silove´m, tepelne´m, elektricke´m apod. Sta´le cˇasteˇji se vsˇak pojem syste´m rozsˇirˇuje na proble´my abstraktnı´. Komponenty syste´mu a jejich interakce jsou potom naprˇ. promeˇnne´ velicˇiny a rovnice, jimzˇ tyto promeˇnne´ vyhovujı´, anebo jinak vyja´drˇene´ prˇ´ıznaky a relace mezi nimi. Zpravidla tento abstraktnı´ syste´m formulujeme tak, aby popsal urcˇitou trˇ´ıdu jevu˚ pozorovany´ch v rea´lne´m sveˇteˇ a podstatneˇ se s nimi shodoval. Takto vytvorˇeny´ abstraktnı´ syste´m pak nahrazuje v nasˇem modelova´nı´ pu˚vodnı´ rea´lne´ jevy a proble´my a sta´va´ se jejich modelem. Kauzalita, vstup, vy´stup Pojem syste´m byl vytvorˇen prˇedevsˇ´ım pro studium slozˇiteˇjsˇ´ıch kauza´lnı´ch vztahu˚, tedy vztahu˚ mezi souborem urcˇity´ch prˇ´ıcˇin a souborem jimi zpu˚sobeny´ch du˚sledku˚. Velicˇiny, ktere´ charakterizujı´ vliv okolı´ na syste´m, jsou vstupy syste´mu. Vneˇjsˇ´ı projevy syste´mu, ktere´ charakterizujı´ vliv syste´mu na jeho okolı´, jsou jeho vy´stupy.
Obra´zek 1: Syste´m a jeho okolı´
2
Matematicke´ modelova´nı´ ve fyzice
2
MATEMATICKY´ MODEL
Dynamicky´ syste´m U veˇtsˇiny kauza´lnı´ch vztahu˚ rea´lne´ho sveˇta pozorujeme, zˇe okamzˇita´ hodnota vy´stupu˚ nenı´ da´na jen soucˇasnou hodnotou vstupu˚, ale i urcˇitou historiı´, tedy stavem vstupu˚ v minulosti. Takovy´ syste´m se nazy´va´ dynamicky´. (Naprˇ. snazˇ´ıme-li se rozjet teˇzˇky´ vozı´k, neza´lezˇ´ı jeho okamzˇita´ rychlost jen na soucˇasneˇ pu˚sobı´cı´ sı´le, ale na integra´lu te´to sı´ly na cele´m intervalu rozjezdu a na pocˇa´tecˇnı´ rychlosti, s nı´zˇ tento rozjezd zacˇ´ınal.) Stav dynamicke´ho syste´mu Stav syste´mu lze charakterizovat pomocı´ tzv. stavovy´ch velicˇin. Stavovou velicˇinou mu˚zˇe by´t naprˇ. teplota syste´mu atp. Vstupy syste´mu spolu se stavovy´mi velicˇinami urcˇujı´ vy´stup syste´mu. Vy´voj syste´mu (tedy vy´voj jeho stavovy´ch velicˇin) lze vyja´drˇit v tzv. stavove´m prostoru. Stavovy´ prostor je zobrazenı´ stavovy´ch velicˇin v za´vislosti na sobeˇ samy´ch. V kazˇde´m cˇasove´m okamzˇiku dostaneme ve stavove´m prostoru jeden bod. Spojenı´m teˇchto bodu˚ dostaneme stavovou trajektorii. Na obra´zku 2 je prˇ´ıklad za´vislosti dvou stavovy´ch velicˇin X1 a X2 na cˇase. Stavova´ trajektorie teˇchto velicˇin je na obra´zku 3.
Obra´zek 2: Stavove´ velicˇiny
Obra´zek 3: Stavova´ trajektorie
2 Matematicky´ model 2.1 Stavova´ rovnice Za´kony dynamicky´ch jevu˚ jsou da´ny rovnicemi typu dxi = fi (x1 (t), ..., xn (t), u1 (t), ..., um (t)) i = 1, 2, ..., n dt
3
(2.1)
Matematicke´ modelova´nı´ ve fyzice
MATEMATICKY´ MODEL
2
kde xi je i-ta´ stavova´ velicˇina a ui je i-ta´ vstupnı´ velicˇina. Pomocı´ teˇchto rovnic mu˚zˇeme modelovat vy´voj stavovy´ch velicˇin v cˇase. Prˇ´ıklady • pro pohyb mechanicke´ soustavy pohybujı´cı´ se vlivem pu˚sobı´cı´ch sil platı´
d hybnost X sı´ly pu˚sobı´cı´ = dt soustavy na soustavu • pro rotujı´cı´ teˇlesa platı´ obdobna´ podmı´nka pro moment hybnosti
(2.2)
momenty sil d moment X = dt hybnosti pu˚sobı´cı´ch na soustavu • pro teplo akumulujı´cı´ se v urcˇite´m prostoru nebo za´sobnı´ku platı´
(2.3)
tepelne´ toky do d akumulovane´ X = dt a z prostoru akumulace teplo
(2.4)
Rovnici 2.1 lze vyja´drˇit i v integra´lnı´m tvaru xi (t) = xi (0) +
Zt
fi (x1 (τ ), ..., xn (τ ), u1 (τ ), ..., um (τ )) dτ i = 1, 2, ...n.
(2.5)
0
Pokud soubor stavovy´ch velicˇin xi (t) zapı´sˇeme jako vektor x = [x1 ... xn ]T
(2.6)
a da´le oznacˇ´ıme vektor vstupu, resp. vektor derivacı´ u = (u1 , ..., um ) ,
(2.7)
f = (f1 , ..., fn ),
(2.8)
resp. mu˚zˇeme rovnici 2.5 zapsat ve tvaru x(t) = x(0) +
Zt 0
4
f(τ )dτ
(2.9)
Matematicke´ modelova´nı´ ve fyzice
2
MATEMATICKY´ MODEL
Pomocı´ teˇchto symbolu˚ bude i soustava 2.1 zapsa´na jedinou vektorovou, tzv. stavovou rovnicı´ dx(t) = f (x(t), u(t)) dt
(2.10)
Vy´stupnı´ rovnice dynamicke´ho syste´mu Stavove´ velicˇiny xi volı´me tak, aby rovnice pro derivace dxi dt
(2.11)
byly tvaru (2.1). Proto se tyto velicˇiny obecneˇ nemusejı´ shodovat s vy´stupy y 1 , ..., yp zkoumane´ho syste´mu. Oznacˇ´ıme-li i vy´stup jako vektor y = [y1 , ..., yp ]T ,
(2.12)
mu˚zˇeme vy´stupnı´ rovnici obecneˇ vyja´drˇit jako funkci y = g(x, u)
(2.13)
jejı´mizˇ argumenty jsou pouze stavove´ a vstupnı´ promeˇnne´. Prˇ´ıklad Na obra´zku 4 je trˇ´ıramenna´ hadicova´ vodova´ha slouzˇ´ıcı´ k urcˇenı´ roviny. Takova´ vodova´ha ma´ vy´razne´ dynamicke´ vlastnosti s ma´lo tlumeny´mi vlastnı´mi kmity hladin vody v koncovka´ch. Popisˇme tento pohyb ve vhodny´ch stavovy´ch promeˇnny´ch. Za´kladnı´ vlastnost stavovy´ch promeˇnny´ch majı´ zrˇejmeˇ vy´sˇky hladin v koncovka´ch h1 ,h2 ,h3 , nebot’ rychlosti jejich zmeˇn souvisı´ s pru˚toky (vyja´drˇeny´mi v jednotka´ch hmotnosti na jednotku cˇasu) ve spojovacı´ch hadicı´ch q1 ,q2 ,q3 rovnicemi
Obra´zek 4: Trˇ´ıramenna´ hadicova´ vodova´ha
dh1 = q2 − q 1 dt dh2 = q3 − q 2 Sρ dt
Sρ
5
(2.14) (2.15)
Matematicke´ modelova´nı´ ve fyzice
2
MATEMATICKY´ MODEL
kde ρ je hustota pouzˇite´ kapaliny. Trˇetı´ stejneˇ sestavena´ rovnice by vsˇak jizˇ nemeˇla by´t do soustavy 2.10 zahrnova´na, nebot’o vy´sˇka´ch hladin platı´, zˇe h1 + h2 + h3 = 3hS
(2.16)
protozˇe objem na´plneˇ vodova´hy zu˚sta´va´ konstantnı´ (hS je urcˇena rovnova´zˇny´m stavem dosazˇeny´m prˇi stejny´ch vy´sˇka´ch hladin). Z toho plyne, zˇe pouze dveˇ (libovolne´) ze trˇ´ı vy´sˇek hladin patrˇ´ı do vektoru x. Trˇetı´ je pak urcˇena jejich konstantnı´m soucˇtem. Ky´va´nı´ na´plneˇ vodova´hy je vy´sledkem soucˇasne´ho vlivu tı´hovy´ch sil, setrvacˇnosti a pasivnı´ch odporu˚. Oznacˇme L de´lku hadicovy´ch spoju˚ a F jejich pru˚tocˇny´ pru˚ˇrez a da´le prˇedpokla´dejme, zˇe L >> hS . Sı´la pu˚sobı´cı´ na na´plnˇ hadicove´ho spoje mezi koncovkami h1 a h2 je FN = pS = (h1 − h3 )ρgS,
(2.17)
kde g je gravitacˇnı´ konstanta. Zmeˇna pru˚toku s cˇasem je tedy FN (h1 − h3 )ρgS dq1 = = , dt m F Lρ
(2.18)
kde m je hmotnost kapaliny v hadicove´m spoji. Dynamiku pru˚toku˚ lze tedy vyja´drˇit rovnicemi dq1 = γF (h1 − h3 ) − Rq1 dt
(2.19)
dq2 = γF (h2 − h1 ) − Rq2 , dt
(2.20)
L L kde
γ=
gS FN
(2.21)
R prˇedstavuje konstantu lamina´rnı´ho odporu v hadicove´m spoji. Ovsˇem prˇi uzavrˇenosti a nemeˇnnosti na´plneˇ vodova´hy nejsou ani pru˚toky q1 , q2 , q3 vza´jemneˇ neza´visle´. Pouze vyrovna´vajı´ rozdı´ly mezi h1 , h2 , h3 a protozˇe neexistuje vneˇjsˇ´ı sı´la, ktera´ by nutila na´plnˇ obı´hat hadicemi, musı´ soucˇet vsˇech trˇ´ı zrychlenı´ by´t nulovy´ dq1 dq2 dq3 + + =0 (2.22) dt dt dt takzˇe trˇetı´ rovnici pro q3 uzˇ nema´ smysl psa´t, byla by jen soucˇtem prvnı´ch dvou. Ovsˇem nejen zrychlenı´, ale i pru˚toky qi musejı´ mı´t nulovy´ soucˇet, protozˇe jsou vy´sledkem integracı´ zrychlenı´ prˇi shodneˇ nulovy´ch pocˇa´tecˇnı´ch podmı´nka´ch, neboli q1 + q2 + q3 = 0
6
(2.23)
Matematicke´ modelova´nı´ ve fyzice
2
MATEMATICKY´ MODEL
a tak stavovy´ vektor tohoto syste´mu budou tvorˇit pouze dveˇ vy´sˇky a dva pru˚toky x = [h1 , h2 , q1 , q2 ]T Po oznacˇenı´ koeficientu˚
1 =a Sρ
γF =K L
(2.24)
R =r L
(2.25)
je stavova´ rovnice 2.10 da´na soustavou cˇtyrˇ rovnic dh1 = a (q2 − q1 ) dt
(2.26)
dh2 = a (−q1 − 2q2 ) dt
(2.27)
dq1 = K (2h1 − 3hS + h2 ) − rq1 (2.28) dt dq2 = K (h2 − h1 ) − rq2 (2.29) dt Na obra´zku 6 je pru˚meˇt stavove´ trajektorie1 do roviny h1 , q2 a za´vislost kmitu˚ h1 a h2 na cˇase. Dany´ model je pocˇ´ıta´n pro parametry a = 2, K = 0, 1 a r = 0, 32.
Obra´zek 5: Pru˚meˇt stavove´ trajektorie h1 − q2 a za´znam kmitu˚ h1 , h2 trˇ´ıramenne´ vodova´hy 1
Vzhledem k tomu, zˇe zde ma´me 4 stavove´ velicˇiny, stavovou trajektorii bychom museli zobrazit ve cˇtyrˇrozmeˇrne´m prostoru. Mu˚zˇeme si ovsˇem zobrazit jejı´ pru˚meˇt do roviny libovolne´ dvojice stavovy´ch velicˇin.
7
Matematicke´ modelova´nı´ ve fyzice
2
MATEMATICKY´ MODEL
2.2 Deduktivnı´ a induktivnı´ formulace modelu Prˇi sestavova´nı´ matematicke´ho modelu stojı´me prˇed u´kolem vystihnout matematicky´mi prostrˇedky urcˇite´ projevy zvolene´ho syste´mu, neboli identifikovat prˇ´ıslusˇny´ syste´m. Rozlisˇujeme dva za´sadneˇ odlisˇne´ principy hleda´nı´ vhodne´ho matematicke´ho modelu, postup deduktivnı´ a postup induktivnı´. Deduktivnı´ postup Dedukcı´ sestavujeme matematicky´ model na za´kladeˇ obecneˇ platny´ch za´konu˚ prˇ´ıslusˇne´ veˇdnı´ disciplı´ny. Tento prˇ´ıstup lze pouzˇ´ıt tam, kde podstata vysˇetrˇovany´ch jevu˚ je dobrˇe zna´ma a kde je propracova´na a oveˇˇrena jejı´ teorie. Jako prˇ´ıklad budeme vysˇetrˇovat pohyb teˇlesa vrzˇene´ho vertika´lneˇ z povrchu Zemeˇ s pocˇa´tecˇnı´ rychlostı´ v0 v gravitacˇnı´m poli Zemeˇ. Vliv ostatnı´ch nebesky´ch teˇles budeme zanedba´vat. Podle Newtonova gravitacˇnı´ho za´kona se dveˇ teˇlesa prˇitahujı´ silou F =κ
Mm , r2
(2.30)
kde κ je gravitacˇnı´ konstanta, M hmotnost Zemeˇ, m hmotnost teˇlesa a r vzda´lenost teˇlesa od strˇedu Zemeˇ. Na za´kladeˇ druhe´ho Newtonova za´kona F =m
dv(t) dt
(2.31)
sestavı´me pohybovou rovnici
M dv = −κ 2 , (2.32) dr r kde v je rychlost a cha´peme ji jako funkci vzda´lenosti od strˇedu Zemeˇ. Z odvozene´ho matematicke´ho modelu lze naprˇ´ıklad stanovit nejmensˇ´ı pocˇa´tecˇnı´ rychlost, prˇi ktere´ se teˇleso jizˇ nevra´tı´ na Zemi atd. v
Induktivnı´ postup Induktivnı´ postup jsme nuceni pouzˇ´ıt tam, kde nejsou splneˇny podmı´nky pro pouzˇitı´ deduktivnı´ho postupu, tj. neexistujı´ nebo nejsou zna´my exaktnı´ za´kony. Vzˇdy je potrˇeba vycha´zet z konkre´tnı´ch jevu˚, na za´kladeˇ ktery´ch (veˇtsˇinou podle zkusˇenosti) sestavı´me model, ktery´ se co nejle´pe shoduje s pozorovany´mi skutecˇnostmi. Tento model sice popisujeme pomocı´ matematicky´ch vztahu˚, ale tyto vztahy nemajı´ rea´lnou fyzika´lnı´ interpretaci. Takto formulovany´ model ma´ zpravidla me´neˇ obecnou platnost, nezˇ model sestaveny´ deduktivnı´m postupem, na cozˇ je prˇi ˇresˇenı´ fyzika´lnı´ch proble´mu˚ nutno bra´t zrˇetel.
8
Matematicke´ modelova´nı´ ve fyzice
3
ROVNOVA´ZˇNY´ STAV DYNAMICKE´HO SYSTE´MU
2.3 Linearita ve stavove´ formulaci Vzhledem k vy´hodnosti linea´rnı´ch vztahu˚ vyuzˇ´ıva´me i u modelu˚ dynamicky´ch jevu˚ co nejvı´ce ty modely, jejichzˇ vztah vstup-vy´stup je linea´rnı´ relacı´. Uva´zˇ´ıme-li, zˇe simulta´nnı´ integrace jsou linea´rnı´ relace, tudı´zˇ o lineariteˇ cˇi nelineariteˇ stavove´ formule (2.10),(2.13) rozhoduje jen charakter vektorovy´ch funkcı´ f a g. Z algebry plyne, zˇe kazˇda´ linea´rnı´ funkce f(x, u) je vzˇdy vyja´drˇitelna´ jako soucˇin urcˇite´ matice tzv. parametru˚ P a vektoru vsˇech argumentu˚ te´to funkce, prˇicˇemzˇ nejen x a u, ale i prvky matice P mohou by´t dany´m zpu˚sobem promeˇnne´ v cˇase. Protozˇe vsˇak x a u uvazˇujeme oddeˇleneˇ, rozdeˇlujeme i matici parametru˚ na dveˇ matice A a B, z nichzˇ matice A typu (n, n) je tzv. matice dynamiky a B (n, m) je tzv. matice buzenı´. Linea´rnı´ stavova´ rovnice dynamicke´ho syste´mu bez zpozˇdeˇnı´ je tedy obecneˇ tvaru dx = A (t) x (t) + B (t) u (t) . (2.33) dt Vy´stupnı´ rovnice 2.13 je rovneˇzˇ da´na maticovy´mi soucˇiny y(t) = g(x, u) = C(t)x(t) + D(t)u(t)
(2.34)
a linea´rnı´ stavovou formulacı´ dynamicke´ho syste´mu ve spojite´m cˇase (bez zpozˇdeˇnı´) pak rozumı´me takovou, jejı´zˇ stavova´ rovnice ma´ tvar (2.33) a vy´stupnı´ rovnice tvar (2.34).
3 Rovnova´zˇny´ stav dynamicke´ho syste´mu Modelem dynamicke´ho syste´mu popisujeme prˇedevsˇ´ım pohyb modelovane´ho objektu, avsˇak prakticky ve vsˇech proble´mech dynamicke´ho chova´nı´ objektu˚ ma´ za´sadnı´ vy´znam take´ ota´zka, kdy a za jaky´ch podmı´nek se vsˇechny stavove´ promeˇnne´ usta´lı´ na konstantnı´ch hodnota´ch tzv. rovnova´zˇne´ho stavu. Jiny´mi slovy, rovnova´zˇny´ stav znamena´ takovou situaci, kdy vsˇechny prˇ´ıcˇiny pohybu objektu se dostanou do vza´jemne´ rovnova´hy a pohyb ve vsˇech stavovy´ch promeˇnny´ch ustane.
3.1 Rovnice statiky dynamicke´ho syste´mu Budeme studovat podmı´nku rovnova´zˇne´ho stavu z hlediska stavove´ rovnice (2.10). Protozˇe vesˇkery´ pohyb syste´mu je vyja´drˇen n derivacemi, v rovnova´zˇne´m stavu je nutne´ dosa´hnout jejich anulova´nı´: dxi = 0, i = 1, 2, ..., n. dt
(3.1)
K tomu, aby splneˇnı´ tohoto pozˇadavku znamenalo rovnova´zˇny´ stav, je nutne´, aby derivace se anulovaly trvalejsˇ´ım zpu˚sobem, a proto tento pozˇadavek ma´ smysl jen za konstantnı´ch hodnot
9
Matematicke´ modelova´nı´ ve fyzice
3
ROVNOVA´ZˇNY´ STAV DYNAMICKE´HO SYSTE´MU
pu˚sobı´cı´ch vstupu˚: ui = uis , i = 1, 2, ..., m.
(3.2)
Souhrnneˇ tyto podmı´nky rovnova´zˇne´ho stavu znamenajı´ soustavu rovnic fi (x, us ) = 0, i = 1, 2, ..., n
(3.3)
f (x, us ) = 0
(3.4)
jejı´zˇ vektorovy´ tvar nazy´va´me rovnice statiky dynamicke´ho syste´mu. Tato rovnice mu˚zˇe by´t velmi rozmanite´ho tvaru a pokud ma´ ˇresˇenı´ (mu˚zˇe jich by´t vı´ce), jsou to sourˇadnice mozˇny´ch rovnova´zˇny´ch stavu˚ syste´mu. Tato ˇresˇenı´ budou vsˇak znamenat mozˇnost rovnova´zˇne´ho stavu jen v prˇ´ıpadeˇ rea´lny´ch korˇenu˚ x = x s .
3.2 Rovnova´zˇny´ stav linea´rnı´ho syste´mu Pro linea´rnı´ dynamicky´ syste´m je i rovnice statiky linea´rnı´. Uplatneˇnı´m podmı´nek (3.3) na linea´rnı´ stavovou rovnici (2.33) dostaneme soustavu linea´rnı´ch algebraicky´ch rovnic Ax + Bus = 0
(3.5)
Te´to rovnici statiky vyhovuje • bud’ jedine´ ˇresˇenı´ xs , pokud matice A je regula´rnı´ • nebo vu˚bec neexistuje ˇresˇenı´, pokud matice A je singula´rnı´, tj. det A = 0 V prvnı´m prˇ´ıpadeˇ ma´ jedine´ ˇresˇenı´ rovnice statiky tvar xs = −A−1 Bus
(3.6)
a pra´veˇ inverze matice A vyzˇaduje splneˇnı´ podmı´nky jejı´ regula´rnosti. Toto ˇresˇenı´ ovsˇem nijak nesignalizuje, zda tento rovnova´zˇny´ stav je cˇi nenı´ stabilnı´.
3.3 Astatismus dynamicke´ho syste´mu V prˇ´ıpadeˇ, zˇe matice A je singula´rnı´, neˇktere´ jejı´ ˇra´dky jsou linea´rneˇ za´visle´ na ostatnı´ch, cozˇ mu˚zˇe by´t zpu˚sobeno dveˇma prˇ´ıcˇinami • bud’ je za´vislost ˇra´dku˚ da´na nadbytecˇny´mi stavovy´mi promeˇnny´mi • nebo dojde k za´vislosti ˇra´dku˚ A i tam, kde mezi stavovy´mi velicˇinami algebraicka´ za´vislost nenı´ - pak jde o astaticky´ dynamicky´ syste´m.
10
Matematicke´ modelova´nı´ ve fyzice
3
ROVNOVA´ZˇNY´ STAV DYNAMICKE´HO SYSTE´MU
V prˇ´ıpadeˇ astaticke´ho syste´mu lze najı´t rea´lna´ cˇ´ısla C1 , C2 , ..., Cn takova´, zˇe pro konstantnı´ vstupy u = us rychlosti stavovy´ch promeˇnny´ch trvale splnˇujı´ podmı´nku n X i=1
Ci
dxi (t) =0 dt
(3.7)
kde alesponˇ jedna z derivacı´ ma´ nenulovou limitu dxi (t) = konst t→∞ dt lim
(3.8)
Kromeˇ trivia´lnı´ho prˇ´ıpadu us = 0 takovy´ syste´m nemu˚zˇe vyhoveˇt podmı´nka´m rovnova´zˇne´ho stavu, nema´ tedy statickou charakteristiku a nazy´va´me jej proto astaticky´.
3.4 Singula´rnı´ body statvove´ rovnice Funkcı´ f(x, u) prave´ strany stavove´ rovnice (2.10) je v kazˇde´m okamzˇiku pohybu syste´mu definova´n vektor rychlosti zmeˇny stavu dx/dt, urcˇujı´cı´ tecˇnu, neboli okamzˇity´ smeˇr stavove´ trajektorie ve stavove´m prostoru Rn . V pru˚meˇtu stavove´ trajektorie do ktere´koli roviny xi , xj je tecˇna tohoto pru˚meˇtu v bodeˇ x da´na derivacı´ ∂xj fj (x, u) = (3.9) ∂xi fi (x, u) Smeˇr stavove´ trajektorie je takto da´n jednoznacˇneˇ v kazˇde´m bodeˇ x (pro soucˇasnou okamzˇitou hodnotu u), nenı´ vsˇak definova´n v bodeˇ, ktery´ splnˇuje podmı´nky rovnova´zˇne´ho stavu (3.3), nebot’ zde derivace (3.9) vedou na neurcˇite´ vy´razy typu „0/0“. Z hlediska pohybu syste´mu ma´ tedy bod xs vyjı´mecˇne´ vlastnosti a je proto v teorii diferencia´lnı´ch rovnic tzv. singula´rnı´m bodem. V pohybu syste´mu to znamena´, zˇe tı´mto bodem procha´zı´ nekonecˇneˇ mnoho stavovy´ch trajektoriı´, zˇe tedy bod xs je hromadny´m bodem vsˇech stavovy´ch trajektoriı´ smeˇˇrujı´cı´ch bud’ k te´muzˇ rovnova´zˇne´mu stavu (stabilnı´ xs ) anebo naopak z neˇj vycha´zejı´cı´ch (nestabilnı´ xs ).
3.5 Singula´rnı´ bod typu uzel a ohnisko Singula´rnı´ body, tj. rovnova´zˇne´ stavy mohou by´t stabilnı´ a nestabilnı´. Kromeˇ toho se rovnova´zˇne´ stavy rozlisˇujı´ take´ podle kmitavosti cˇi nekmitavosti v okolı´ xs . Neˇktere´ rovnova´zˇne´ stavy jsou charakteristicke´ tı´m, zˇe syste´m se v nich ustaluje nekmitavy´m pohybem (stabilnı´ prˇ´ıpad) nebo se nekmitavy´m pohybem od nich vzdaluje (nestabilnı´ prˇ´ıpad). Takovy´ singula´rnı´ bod se nazy´va´ stabilnı´, resp. nestabilnı´ uzel. Dalsˇ´ım typem xs je stabilnı´, resp. nestabilnı´ ohnisko vyznacˇujı´cı´ se tı´m, zˇe ustalova´nı´ ve stabilnı´m xs , resp. vzdalova´nı´ se od nestabilnı´ho xs se deˇje pohybem kmitavy´m, prˇi neˇmzˇ se opakovaneˇ meˇnı´ zname´nka derivacı´ dxi /dt, i = 1, ..., n.
11
4 NUMERICKE´ MODELOVA´NI´
Matematicke´ modelova´nı´ ve fyzice
Obra´zek 6: Stavove´ trajektorie pohybu v okolı´ singula´rnı´ho bodu
4 Numericke´ modelova´nı´ Ve veˇtsˇineˇ prakticky´ch u´loh je ˇresˇenı´ stavovy´ch (2.10) a vy´stupnı´ch (2.13) rovnic nemozˇne´ nebo velmi obtı´zˇne´. Z tohoto du˚vodu jsou tyto rovnice ˇresˇeny numericky pomocı´ vy´pocˇetnı´ techniky.
4.1 Diskretizace spojite´ho syste´mu numerickou metodou rˇesˇenı´ Prˇi numericke´m ˇresˇenı´ stavovy´ch a vy´stupnı´ch rovnic je nutne´ tyto rovnice diskretizovat, tj. zvolit prˇimeˇˇreny´ interval ∆t a tı´m prˇejı´t od spojite´ho cˇasu t k diskre´tnı´mu cˇasu k. Mı´sto spojiteˇ promeˇnny´ch pu˚vodnı´ho modelu na pocˇ´ıtacˇove´m modelu pak pracujeme jen s odpovı´dajı´cı´mi posloupnostmi jejich vzorku˚. Aplikacı´ numericke´ metody se takto k pu˚vodneˇ spojite´mu syste´mu prˇirˇazuje novy´ diskre´tnı´ 12
4 NUMERICKE´ MODELOVA´NI´
Matematicke´ modelova´nı´ ve fyzice
syste´m tak, aby posloupnosti vzorku˚ jeho stavovy´ch velicˇin se co nejle´pe shodovaly s odezvou stavu syste´mu pu˚vodnı´ho. Podstatu diskre´tnı´ho nahrazenı´ numerickou metodou ˇresˇenı´ si uka´zˇeme na nejjednodusˇsˇ´ı Euleroveˇ metodeˇ. Rˇesˇenı´ pocˇa´tecˇnı´ u´lohy stavove´ rovnice (2.10) pro pocˇa´tecˇnı´ podmı´nku x(0) = x 0 a pro dany´ pru˚beˇh vstupu u = u(k) vycha´zı´ u te´to metody z nejjednodusˇsˇ´ı prˇiblizˇne´ na´hrady derivacı´ doprˇedny´mi diferencˇnı´mi podı´ly
dxi (t) . xi (k + 1) − xi (k) = fi (x(k), u(k)) = dt t=tk ∆t
(4.1)
„Nove´“ sourˇadnice stavu v na´sledujı´cı´m okamzˇiku diskre´tnı´ho cˇasu pocˇ´ıta´me tedy z „dosavadnı´ch“ x(t) a u(k) diferencˇnı´mi rovnicemi xi (k + 1) = xi (k) + ∆t fi (x(k), u(k)) , i = 1, 2, ..., n
(4.2)
ktere´ jsou stavovy´mi rovnicemi diskre´tnı´ho dynamicke´ho syste´mu nahrazujı´cı´ho pu˚vodnı´ spojity´ (2.10).
4.2 Chyby numericky´ch vy´pocˇtu˚ Prˇi numericke´m ˇresˇenı´ matematicky´ch modelu˚ vznikajı´ dva druhy chyb: 1. Aproximacˇnı´ chyba - tato chyba souvisı´ s neprˇesnostı´, s jakou zvolena´ numericka´ metoda aproximuje prˇesne´ ˇresˇenı´. Meˇˇr´ıtkem toho, jak prˇesneˇ aproximujeme skutecˇnou funkci na dane´m intervalu je globa´lnı´ (celkova´) chyba, ktera´ je u Eulerovy metody u´meˇrna´ zvolene´mu ∆t. Cˇ´ım kratsˇ´ı je tento cˇasovy´ interval, tı´m prˇesneˇji aproximujeme danou funkci. Z hlediska aproximacˇnı´ chyby je tedy vy´hodne´ volit integracˇnı´ krok ∆t co nejmensˇ´ı. 2. Zaokrouhlovacı´ chyba - tato chyba je da´na tı´m, zˇe pocˇ´ıtacˇ pracuje vzˇdy s urcˇity´m pocˇtem desetinny´ch mı´st a tedy i s urcˇitou prˇesnostı´. Velikost zaokrouhlovacı´ chyby je neprˇ´ımo u´meˇrna´ velikosti integracˇnı´ho kroku ∆t, protozˇe cˇ´ım mensˇ´ı je jeho hodnota, tı´m na dane´m intervalu provedeme vı´ce matematicky´ch operacı´ a tı´m se na´m zaokrouhlovacı´ chyba akumuluje. 3. U´hrnna´ chyba - je to celkova´ chyba, ktere´ se dopousˇtı´me prˇi vy´pocˇtu matematicke´ho modelu pomocı´ vy´pocˇetnı´ techniky. Je da´na soucˇtem chyby globa´lnı´ a zaokrouhlovacı´. Na obra´zku 7 je zna´zorneˇny´ pru˚beˇh aproximacˇnı´ (e1), zaokrouhlovacı´ (e2) a celkove´ (e) chyby. Existuje jista´ optima´lnı´ hodnota ∆t, prˇi ktere´ je celkova´ chyba nejmensˇ´ı. Prˇi urcˇova´nı´ kroku, s jaky´m bude vy´pocˇet probı´hat se snazˇ´ıme te´to hodnoteˇ co nejvı´ce prˇiblı´zˇit.
13
5 OVEˇRˇENI´ FUNKCˇNOSTI MODELU
Matematicke´ modelova´nı´ ve fyzice
Obra´zek 7: Pru˚beˇh aproximacˇnı´ (e1), zaokrouhlovacı´ (e2) a celkove´ (e) chyby
5 Oveˇrˇenı´ funkcˇnosti modelu Po sestavenı´ matematicke´ho modelu je potrˇeba oveˇˇrit jeho funkcˇnost, tj. oveˇˇrit, jestli dany´ model vystihuje rea´lne´ chova´nı´ syste´mu, jestli jsme do modelu zahrnuli vsˇechny relevantnı´ velicˇiny ovlivnˇujı´cı´ chova´nı´ syste´mu atd. Da´le je potrˇeba prove´st vhodne´ nastavenı´ vy´pocˇetnı´ch konstant tak, aby se vy´stup modelu co nejvı´ce prˇiblı´zˇil vy´stupu rea´lne´ho syste´mu.
5.1 Optimalizace parametru˚ Jestlizˇe budeme vycha´zet z prˇedem zna´me´ nebo prˇedpokla´dane´ struktury matematicke´ho modelu, potom musı´me nale´zt takove´ hodnoty parametru˚ pocˇ´ıtacˇove´ho modelu, aby shoda rea´lne´ho syste´mu a pocˇ´ıtacˇove´ho modelu byla maxima´lnı´. Takova´ u´loha se oznacˇuje jako optimalizace parametru˚ dynamicke´ho syste´mu. Prˇedpokla´dejme, zˇe stejny´ vstupnı´ signa´l pu˚sobı´ na vstup rea´lne´ho syste´mu a take´ na vstup pocˇ´ıtacˇove´ho modelu. Vy´stupnı´ hodnoty, ktere´ budou zmeˇˇreny na raa´lne´m syste´mu oznacˇ´ıme yexp . Tyto hodnoty se se dı´ky urcˇity´m chyba´m neshodujı´ s vy´sledky pocˇ´ıtacˇove´ho modelu ymodel . Na za´kladeˇ urcˇenı´ teˇchto velicˇin je nutno vypocˇ´ıtat, jak se lisˇ´ı chova´nı´ rea´lne´ho syste´mu od chova´nı´ pocˇ´ıtacˇove´ho modelu. Definujme velicˇinu e, ktera´ bude prˇedstavovat odchylku mezi modelem a rea´lny´m procesem. ei = yexp (ti ) − ymodel (ti ) , i = 1, ..., n
14
(5.1)
5 OVEˇRˇENI´ FUNKCˇNOSTI MODELU
Matematicke´ modelova´nı´ ve fyzice
Na za´kladeˇ odchylek ei je nutno matematicky´ model upravit tak, aby se jeho vy´stupy co nejvı´ce shodovaly s vy´stupy rea´lne´ho fyzika´lnı´ho syste´mu. Odchylku mezi modelem a rea´lny´m syste´mem definujeme pomocı´ tzv. kriteria´lnı´ funkce S (yexp , ymodel ), ktera´ musı´ mı´t takovou vlastnost, zˇe pro nejlepsˇ´ı matematicky´ model ma´ minima´lnı´ hodnotu. V praxi se nejcˇasteˇji pouzˇ´ıva´ krite´rium strˇednı´ kvadraticke´ odchylky, ktere´ je definovane´ na´sledujı´cı´m vztahem: S (yexp , ymodel ) =
n X i=1
(yexp − ymodel )2
(5.2)
Cı´lem je najı´t takove´ hodnoty parametru˚ modelu, prˇi ktery´ch kriteria´lnı´ funkce dosahuje co nejmensˇ´ıch hodnot.
5.2 Validace modelu Nelze povazˇovat za dostatecˇne´, je-li matematicky´ model pouze verifikova´n, tj. nalezeny takove´ hodnoty konstant, prˇi ktery´ch se jeho chova´nı´ shoduje s chova´nı´m rea´lne´ho syste´mu, na ktere´m oveˇˇrujeme jeho funkcˇnost. Te´meˇˇr kazˇdy´ slozˇity´ matematicky´ model dovoluje nale´zt takove´ kombinace parametru˚, ktere´ umozˇnˇujı´ dobrou shodu s urcˇity´m experimentem. Tato shoda vsˇak zdaleka nezarucˇuje matematicke´mu modelu schopnost prˇedpovı´dat skutecˇnost i za jiny´ch podmı´nek. Proto musı´me s matematicky´m modelem prova´deˇt tzv. validaci. Validace je testova´nı´, zda se matematicky´ model chova´ za ru˚zny´ch vstupnı´ch podmı´nek podobneˇ jako rea´lny´ syste´m. Mu˚zˇeme rozlisˇovat dva zpu˚soby validace: 1. Vy´sledky zı´skane´ ze simulacˇnı´ch experimentu˚ se porovna´vajı´ s experimenta´lnı´mi daty, ktere´ jsme zı´skali pozorova´nı´m rea´lne´ho syste´mu za urcˇite´ obdobı´ v minulosti a za podmı´nek, ktere´ byly pouzˇite´ prˇi simulacˇnı´ch experimentech. Data, ktera´ pouzˇ´ıva´me k validaci, vsˇak nesmı´ by´t pouzˇita stavbeˇ a optimalizaci matematicke´ho modelu. 2. Na za´kladeˇ vy´sledku˚ simulace se prˇedpovı´da´ chova´nı´ rea´lne´ho syste´mu v budoucnosti. Tyto vy´sledky se porovna´vajı´ s u´daji, ktere´ se zı´ska´vajı´ z na´sledne´ho pozorova´nı´ rea´lne´ho syste´mu za urcˇite´ obdobı´.
15
6 SEZNAM POUZˇITE´ LITERATURY
Matematicke´ modelova´nı´ ve fyzice
6 Seznam pouzˇite´ literatury [1]
Pazourek J.: Simulace biologicky´ch syste´mu˚, Praha, Grada 1992
[2]
Zı´tek P.,Petrova´ R.: Matematicke´ a Simulacˇnı´ modely, Praha Vydavatelstvı´ CˇVUT 1996 (skriptum)
16
