MOKSLASplius.lt

Briuseliatorius

Dalyvaujančių cheminėje reakcijoje medžiagų tankiai reakcijos metu paprastai didėja arba mažėja, todėl praėjus pakankamai ilgam laiko tarpui reakcija užgęsta. Tačiau yra ir tokių cheminių reakcijų, kuriose reaguojančių medžiagų tankiai osciliuoja, t.y. reakcija periodiškai vyksta tai į vieną, tai į kitą pusę [Zaikin70]. Čia panagrinėsime osciliuojančią cheminę reakciją, kurią atrado Briuselio mokslininkai I.Prigogine, G.Nicolis ir R.Lefever [Prigogine67,Prigogine68,Lefever67]. Briuselio mokyklos bei miesto garbei ji buvo pavadinta briuseliatoriumi. Briuseliatorių aprašysime netiesinių spartuminių lygčių sistema, kuri nusako reakcijos produktų tankio kitimą laikui bėgant. Panašias, tiesa, paprastesnio elgesio spartumines lygtis jau esame aptarę Monte Carlo metodo. eksperimente . Skaitiškai rasime briuseliatoriaus lygčių osciliuojančius ir relaksuojančius sprendinius, o pasitelkę mažų nuokrypių nuo pusiausvyros analizės metodą įvertinsime šios cheminės reakcijos stabilumo sąlygas.

Nustatymai

Pirmiausia nustatome atsakymų šrifto dydį šiame puslapyje taip, kad gerai matytume atsakymus
Šrifto dydis:

  

Briuseliatoriaus lygtys

Tegul raidės $ X $, $ Y $, $ A $, $ B $, $ D $ ir $ E $ žymi cheminėje reakcijoje dalyvaujančias medžiagas, o kai tai nesukels dviprasmiškumo, — ir jų tankius. Briuseliatorių [Tyson73a,Tyson73] nusako tokios keturios elementarios reakcijos ir jų spartos koeficientai $ k_i $:

  1. $ A\rightarrow X $ reakcijos sparta $ k_1 $;
  2. $ B+X\rightarrow Y+D $ reakcijos sparta $ k_2 $;
  3. $ 2  X+Y\rightarrow 3X $ reakcijos sparta $ k_3 $;
  4. $ X\rightarrow E $ reakcijos sparta $ k_4 $.

Laikysime, kad koeficientai $ k_i $ yra iš anksto žinomi dydžiai. Jie nepriklauso nuo reakcijoje dalyvaujančių medžiagų kiekio ar jų santykio mišinyje. $ A $ ir $ B $ yra dviejų pradinių medžiagų tankiai. Juos visos reakcijos metu palaikysime pastoviais. $ D $ ir $ E $ yra galinių produktų, o $ X $ ir $ Y $ — tarpinių produktų tankiai. Sudarysime spartumines lygtis, nusakančias cheminių medžiagų tankio kitimą vykstant reakcijai. Konkretaus reakcijos produkto kitimo sparta, t.y. jo tankio išvestinė pagal laiką $ T $, yra proporcinga sąveikojančių medžiagų tankiams. Pavyzdžiui, $ Y $ produkto kiekis keičiasi tik 2 ir 3 reakcijose. Į pastarąsias įeina $ B $ ir $ X $ medžiagų tankiai. Todėl $ Y $ tankio kitimą laike, t.y. dydį $ \mathrm{d} Y/\mathrm{d} T $, nusako dviejų narių suma:

\[ \mathrm{d} Y/\mathrm{d}  T = k_2 B  X- k_3 X^2 Y.\tag{1} \]

Kadangi antroji reakcija didina $ Y $ medžiagos tankį, o trečioji mažina, spartuminėje lygtyje pirmąjį narį parašėme su pliuso, o antrąjį — su minuso ženklu. Pasinaudoję keturiomis reakcijų išraiškomis bei panašiai samprotaudami užrašome spartumines likusių medžiagų $ X $, $ D $ ir $ E $ lygtis:

\[ \begin{split}\mathrm{d} X/\mathrm{d}  T&= k_1 A - k_2 B  X+(3-2) k_3 X^2 Y- k_4 X,\\\mathrm{d} D/\mathrm{d}  T&= k_2 B X,\\\mathrm{d} E/\mathrm{d}  T&= k_4 X .\end{split}\tag{2} \]

Medžiagos $ X $ kitimo spartą nusako net keturi nariai, nes ši medžiaga dalyvauja visose keturiose reakcijose. Supaprastinsime gautas lygtis. Kadangi tankiai $ A $ ir $ B $ nuo laiko nepriklauso, reakcijų pastoviąsias $ k_1 $, $ k_2 $ ir $ k_3 $ galima įtraukti į tankius $ A $ ir $ B $. Naujuosius tankius, kuriuos dar papildomai padalysime iš $ k_3 $, toliau taip pat žymėsime tomis pačiomis raidėmis: $ k_1 A/k_3 \rightarrow A $, $ k_2 B/k_3 \rightarrow B $. Be to, įvesime bedimensį laiką $ T k_3 \rightarrow t $ ir naują koeficientą $ k_4/ k_3\rightarrow k $. Dvi paskutiniosios diferencialinės lygtys nusako galutinių produktų susidarymą. Jas lengvai išspręsime po to, kai rasime $ X $ medžiagos tankio kitimą laike. Tokiu būdu briuseliatorių iš esmės aprašo pirmosios dvi redukuotos diferencialinės lygtys. Užrašytos naujaisiais kintamaisiais jos atrodo taip:

\[ \begin{split}\mathrm{d} X/\mathrm{d} t&=A-B  X+ X^2 Y-k  X\\\mathrm{d} Y/\mathrm{d} t&=B  X-X^2 Y.\end{split}\tag{3} \]

Tai netiesinė diferencialinių lygčių sistema. Netiesiškumą įneša tankių kvadratai ir sandaugos. Panašaus pavidalo nariai pasirodo daugelyje nestabilaus ir chaotinio elgesio sistemų. Užrašytas lygtis spręsime skaitiškai. Tam reikia nurodyti pradinių medžiagų $ A $, $ B $ tankius, tarpinių produktų $ X $, $ Y $ tankius pradiniu laiko momentu, $ X(0)={X_0} $ ir $ Y(0)={{Y}_0} $, bei koeficientą $ k=k_4/k_3 $. Chemines osciliacijas stebėsime, jei pradinių medžiagų tankiams ir koeficientui $ k $ suteiksime tokias vertes:

\boldmath$oscParametrai=\bigl\{$
\boldmath$A$\boldmath$\to$pradinės medžiagos $ A $ tankis
\boldmath$B$\boldmath$\to$pradinės medžiagos $ B $ tankis
\boldmath$k$\boldmath$\to$spartos koeficientų santykis
\boldmath$\bigr\}$

  

Lygtis spręsime \boldmath$NDSolve[~]$ nuo pradinio $ t=0 $ iki galinio $ t=\mathbf{tend} $ laiko momento.

\boldmath\begin{eqnarray*}&&sprendinys=\mathrm{NDSolve}[\{X'[t]==Y[t] X[t]^2-B X[t]-k X[t]+A,Y'[t]==B X[t]-X[t]^2 Y[t],\\&&\quad X[0]==\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{X_0}$},Y[0]==\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{Y_0}$}\}/.oscParametrai,\{X[t],Y[t]\},\{t,0,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{tend}$}\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys}$}]\end{eqnarray*}


parinktys =
$ X_0 $=
$ Y_0 $=
tend=

  

Sprendinius pavaizduojame komanda \boldmath$\mathrm{Plot}[~]$. Komandos \boldmath$\mathrm{Show}[~]$ ir \boldmath$\mathrm{GraphicsArray}[~]$ padeda grafikus surikiuoti vieną šalia kito.

\boldmath\begin{eqnarray*}&&\mathrm{Show}[\mathrm{GraphicsRow}[\{\\&&\qquad \mathrm{Plot}[\mathrm{Evaluate}[X[t]/.sprendinys], \{t, 0, tend\},      \mathrm{AxesLabel}\to \{t, X\}, \mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys1}$}],\\&&\qquad \mathrm{Plot}[\mathrm{Evaluate}[Y[t]/.sprendinys], \{t, 0, tend\},      \mathrm{AxesLabel}\to \{t, Y\}, \mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys2}$}]\\&&\}],\mathrm{DisplayFunction}\to \$\mathrm{DisplayFunction},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys3}$}];\end{eqnarray*}

Pirmojo grafiko piešimo parinktys:
Antrojo grafiko piešimo parinktys:
Grafikų grupės piešimo parinktys:

  

Matome, kad tarpinių produktų $ X $ ir $ Y $ tankiai reakcijos metu osciliuoja. Osciliacijų pobūdis nepriklauso nuo pradinių sąlygų. Nuo jų priklauso tik pradinio pereinamojo proceso eiga. Pastarasis mūsų atveju yra gana trumpas. Tai geriau matysime, jei atvaizduosime $ Y $ medžiagos tankio priklausomybę nuo $ X $.

\boldmath\begin{eqnarray*}&&\mathrm{ParametricPlot}[\\&&\ \ \mathrm{Evaluate}[\{X[t],Y[t]\}/.sprendinys],\{t,0,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{t_{max}}$}\},\mathrm{AxesLabel}\to \{X, Y\},\\&&\qquad \mathrm{Epilog} \to \{\mathrm{PointSize}[0.03], \mathrm{Point}[\{X_0,Y_0\}]\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys}$}]\end{eqnarray*}

Piešti iki \boldmath$t_{max}$=
Piešimo parinktys:

  

$ X-Y $ plokštumoje parametriškai nuo laiko priklausanti trajektorija prasideda centriniame taške (jį pavaizdavome papildoma parinktimi \boldmath$\mathrm{Epilog}\to\{\mathrm{PointSize}[0.03], mathrm{Point}[\{X_0,\,Y_0\}]\}$) ir laikui bėgant greitai pereina į ribinį, netaisyklingo trikampio formos ciklą. Įdomu tai, kad nepriklausomai nuo pradinio taško vietos $ X-Y $ plokštumoje, t.y. nuo pradinių sąlygų, kreivės trajektorija visada susisuka į tą patį ribinį ciklą.

Stabilumo analizė

Kadangi netiesinės diferencialinės lygtys analiziškai išsprendžiamos tik išimtinais atvejais, apie sistemos evoliuciją dažnai tenka spręsti iš skaitinių sprendinių. Deja, turėdami skaitinius sprendinius, niekada negalime būti tikri, ar išnagrinėjome visus įdomius atvejus. Todėl bet kokia informacija apie diferencialinės lygties analizines savybes visada išlieka svarbi. Pasirodo, kad apie reakcijos dinamiką kai ką galima sužinoti vien iš spartuminių lygčių sandaros, pačių lygčių visai nesprendžiant.

Tiesinių ir netiesinių sistemų dinamikoje svarbu žinoti sistemos savybes arti sistemos pusiausvyros taškų. Mažus nukrypimus nuo pusiausvyros tada jau galima aprašyti tiesinėmis diferencialinėmis lygtimis. Pastarųjų sprendiniai yra matematikų išnagrinėti ir suklasifikuoti. Kaip matysime iš briuseliatoriaus lygčių, netiesinės diferencialinės lygties ištiesinimas, t.y. jos pakeitimas tiesine lygtimi arti pusiausvyros taškų, leidžia kai ką pasakyti apie sistemos stabilumą šių taškų aplinkoje.

Įsivaizduokime, kad parinkome tokį parametrų $ A $, $ B $ ir $ k $ rinkinį, kurį įstatę į lygtis gauname gęstančius briuseliatoriaus sprendinius. Tai reikštų, kad praėjus pakankamai ilgam laikui tarpinių medžiagų $ X $ ir $ Y $ kiekis nustoja keitęsis. Matematiškai šį teiginį suformuluotume dešiniąsias mūsų diferencialinių lygčių puses prilyginę nuliams: $ \mathrm{d} X/\mathrm{d} t=0 $ ir $ \mathrm{d} Y/\mathrm{d} t=0 $.

\boldmath\begin{eqnarray*}&&pusiausvyra=\{0==Y[t] X[t]^2-B X[t]-k X[t]+A,\\&&\hphantom{pusiausvyra=\{} 0==B X[t]-X[t]^2 Y[t]\}/.\{X[t]\to X_e,Y[t]\to Y_e\}\end{eqnarray*}


Komanda-programa \boldmath$\mathrm{Solve}[~]$ išsprendę algebrinę lygčių sistemą randame pusiausvirąjį sprendinį.

\boldmath\begin{eqnarray*}sprendinysPusiausvyroje=\mathrm{Flatten}[\mathrm{Solve}[pusiausvyra,\{X_e,Y_e\}]]\end{eqnarray*}


Tokiu būdu, kai briuseliatorius yra pusiausvyroje, tarpinių produktų tankiai išsireiškia per pradinių medžiagų tankius $ A $ ir $ B $ bei spartos koeficientą $ k $: $ X_e=A/k $ ir $ Y_e=B\, k/A $. Imkime nedidelius tankių nukrypimus $ \Delta X(t) $ ir $ \Delta Y(t) $ nuo pusiausvyros padėties, kurie tenkina nelygybes $ \Delta X\ll X_e $ ir $ \Delta Y \ll Y_e $. Nuokrypio $ \Delta  $ ženklą (didžioji graikiška delta) gausite surinkę klaviatūroje \boldmath$\boldsymbol{\backslash}[\mathrm{CapitalDelta}]$, arba, trumpiau, Esc D Esc.,Dešinėse diferencialinių lygčių pusėse vietoje $ X(t) $ ir $ Y(t) $ įstatykime artimas pusiausvirajam sprendiniui išraiškas $ X_e+\Delta X(t) $ ir $ Y_e+\Delta Y(t) $. Kadangi nuokrypius imame mažus, didžiausią įnašą duos tik nariai, proporcingi $ \Delta X $ ir $ \Delta Y $. Šiems dviems nariams išskirti pasinaudosime nedidele gudrybe: $ \Delta X(t) $ ir $ \Delta Y(t) $ papildomai padauginsime iš formalaus ,,mažo'' dydžio $ h $. Jo atžvilgiu komanda \boldmath$\mathrm{Series}[~]$ skleiskime eilute iki norimo mažumo narių (šiuo atveju iki $ \mathrm{O}(h)^2 $). Formali ,,skleidimo'' procedūra čia reikalinga tik aukštesnės eilės nariams numesti. Tada komanda \boldmath$\mathrm{Normal}[~]$ skleidinį vėl paversime įprastine išraiška, o formalų parametrą $ h $ pašalinsime, jį tiesiog prilyginę vienetui. Gautose išraiškose narius prie $ \Delta X(t) $ ir $ \Delta Y(t) $ sugrupuosime. Pradžioj įstatykime $ X_e+\Delta X(t) $ ir $ Y_e+\Delta Y(t) $ į briuseliatoriaus diferencialinių lygčių dešiniąsias puses.

\boldmath\begin{eqnarray*}tarpinis=\mathrm{Expand}[\{Y[t] X[t]^2-B X[t]-k X[t]+A, B X[t]-X[t]^2 Y[t]\}/.\{X[t]\to X_e+h \Delta X[t],          Y[t]\to Y_e+h \Delta Y[t]\}]\end{eqnarray*}


Išskyrę narius, proporcingus $ \Delta X $ ir $ \Delta Y $, turime

\boldmath\begin{eqnarray*}\mathrm{Collect}[\mathrm{Normal}[\mathrm{Series}[tarpinis,\{h,0,1\}]]/.h\to 1,\Delta X[t]]\end{eqnarray*}


Iš pirmosios išraiškos matome, kad pirmos eilės nariai yra $ (-B-k+2  X_e Y_e) \Delta X $ ir $ X_e^2 \Delta Y $. Pastebėsime, kad gautieji nariai (koeficientai) sutampa su diferencialinių lygčių dešiniųjų pusių pirmosiomis dalinėmis išvestinėmis, apskaičiuotomis pusiausviriesiems sprendiniams $ X_e $ ir $ Y_e $. Todėl reikalingus narius taip pat galima ,,išrinkti'' ir diferencialiniu operatoriumi \boldmath$\mathrm{D}[~]$. Tam pakanka apskaičiuoti atitinkamų lygčių dalines išvestines $ X(t)=X_e $, $ Y(t)=Y_e $ taškuose.

\boldmath\begin{eqnarray*}m11=\frac{\partial [Y[t] X[t]^2-B X[t]-k X[t]+A]}{\partial X[t]}/.\{X[t]\to X_e,Y[t]\to Y_e\}\end{eqnarray*}


\boldmath\begin{eqnarray*}m12=\frac{\partial [Y[t] X[t]^2-B X[t]-k X[t]+A]}{\partial Y[t]}/.\{X[t]\to X_e,Y[t]\to Y_e\}\end{eqnarray*}


\boldmath\begin{eqnarray*}m21=\frac{\partial [B X[t]-X[t]^2 Y[t]]}{\partial X[t]}/.\{X[t]\to X_e,Y[t]\to Y_e\}\end{eqnarray*}


\boldmath\begin{eqnarray*}m22=\frac{\partial [B X[t]-X[t]^2 Y[t]]}{\partial Y[t]}/.\{X[t]\to X_e,Y[t]\to Y_e\}\end{eqnarray*}


Dabar iš mažų pokyčių $ \Delta X $ ir $ \Delta Y $ sudarysime vektorių $ \vec{v} $, o iš koeficientų $ m_{ij} $ — matricą\boldmath$m$.

\boldmath\begin{eqnarray*}v=\{\Delta X[t],\Delta Y[t]\};\quad m=\{\{m11, m12\},\{m21, m22\}\};\end{eqnarray*}

Linearizuotą diferencialinių lygčių sistemą užrašysime tokiu matriciniu pavidalu: $ \mathrm{d} \vec{v}/\mathrm{d} t = m \vec{v} $.

\boldmath\begin{eqnarray*}(linearizuotosLygtys = \mathrm{Thread}[\mathrm{Equal}[\mathrm{D}[v,t] ,m . v]])//\mathrm{TableForm}\end{eqnarray*}


Matrica $ m $ yra vadinama Jacobi matrica arba jakobianu. Ji vaidina nepaprastai svarbų vaidmenį netiesinių sistemų stabilumo analizėje [Moon87]. Linearizuotų diferencialinių lygčių sprendinių ieškosime kompleksinių eksponenčių pavidalu: $ \Delta  X(t)=\Delta X_C \mathrm{e}^{s t} $ ir $ \Delta Y(t)=\Delta Y_C \mathrm{e}^{s t} $, kur $ \Delta X_C $, $ \Delta Y_C $ ir $ s $ yra nuo laiko nepriklausantys parametrai. Tiesą sakant, mūsų tikslas yra ne rasti ,,ištiesintų'' diferencialinių lygčių sprendinius (\boldmath$\mathrm{DSolve}[~]$ jas nesunkiai išspręs), o tik nustatyti parametrą $ s $, kuris vadinamas Liapunovo rodikliu. Būtent šio parametro vertė, kaip tuojau matysime, ir charakterizuoja sprendinių stabilumą.

\boldmath\begin{eqnarray*}sprendiniai=\{\Delta X[t]->\Delta X_C*\mathrm{Exp}[s* t],\Delta Y[t]->\Delta Y_C*\mathrm{Exp}[s*t]\}\end{eqnarray*}


Sprendinius, kurių koeficientus ir parametrą $ s $ dar turime rasti, įstatykime į linearizuotas lygtis. Abi lygybės puses padauginkime iš $ \mathrm{e}^{-s t} $, kad suprastintume bendrą eksponentinį daugiklį $ \mathrm{e}^{s t} $.

\boldmath\begin{eqnarray*}(linearizuotosLygtysKoeficientams=\mathrm{Thread}[\mathrm{Equal}[\mathrm{Exp}[-s*t]*\mathrm{D}[v/.sprendiniai,t],\\\hphantom{linearizuotosLygtysKoeficientams=\mathrm{Thread}[\mathrm{Equal}}\mathrm{Expand}[\mathrm{Exp}[-s*t]*m .v /.sprendiniai]]])//\mathrm{TableForm}\end{eqnarray*}


Pasinaudojus matrica \boldmath$m$ ir vektoriumi $ \vec{v} $, šią tiesinių algebrinių lygčių sistemą galima perrašyti matriciniu pavidalu:

\[     \begin{pmatrix}    m_{11} - s &m_{12} \\    m_{21}&m_{22} -  s  \end{pmatrix}\begin{pmatrix}    \Delta X_C \\    \Delta Y_C\end{pmatrix}=\begin{pmatrix}    0 \\    0\end{pmatrix}.\tag{4} \]
\boldmath\begin{eqnarray*}(\{\{m11, m12\}, \{m21, m22\}\} - s*\mathrm{IdentityMatrix}[2])\, .\, (v /. \{\Delta X[t] \to \Delta X_C, \Delta Y[t] \to \Delta Y_C\}) // \mathrm{Expand}) //\mathrm{TableForm}\end{eqnarray*}


Kadangi parametrai $ \Delta  X_C, \Delta Y_C $ yra nepriklausomi, sprendžiamoji lygtis turės netrivialius (t.y. nelygius nuliui) sprendinius tik tuo atveju, kai matricos ($\mathbf{m} - s \mathbf{1}$) determinantas

\boldmath\begin{eqnarray*}determinantas=\mathrm{Det}(\{\{m11, m12\}, \{m21, m22\}\}]\end{eqnarray*}


bus lygus nuliui (\boldmath$1$ žymi vienetinę matricą). O tai gali atsitikti tik esant šioms parametro $ s $ vertėms:

\boldmath\begin{eqnarray*}\mathrm{Flatten}[\mathrm{Solve}[determinantas==0,s]]\end{eqnarray*}


Nagrinėjamą algebrinę tiesinių lygčių sistemą (4) buvo galima suformuluoti ir kaip tikrinių verčių uždavinį: $ m \vec{v}=s\vec{v} $. Tada $ s $ interpretuotume kaip tikrinę vertę, o $ \vec{v} $ — kaip tikrinį vektorių. Tikrinių verčių uždaviniams spręsti Mathematica turi tris komandas: \boldmath$\mathrm{Eigenvalues}[~]$, \boldmath$\mathrm{Eigenvectors}[~]$ ir \boldmath$\mathrm{Eigensystem}[~]$. Pavyzdžiui, komanda \boldmath$\mathrm{Eigenvalues}[~]$ nesunkiai randa tas pačias tikrines reikšmes $ s=s_1 $ ir $ s=s_2 $, sutampančias su anksčiau gautomis:

\boldmath\begin{eqnarray*}\mathrm{Eigenvalues}[m]\end{eqnarray*}


Linearizuoto sprendinio stabilumą lemia eksponentės rodiklio realiosios dalies ženklas $ \mathrm{Re}(s) $. Jei bent vieno iš sprendinių realioji dalis, $ \mathrm{Re}(s_1) $ arba $ \mathrm{Re}(s_2) $, yra teigiamas skaičius, tada sistemos evoliucija taško $ (X_e,Y_e) $ aplinkoje bus nestabili (mažo trikdžio paveiktos artimos trajektorijos tolsta viena nuo kitos eksponentiškai greitai). Linearizuota sistema osciliuos tik tuo atveju, kai $ s $ bus grynai menamas dydis. Aišku, norėdami spręsti apie globalią sistemos evoliuciją turime rasti pirminės netiesinės diferencialinės lygties sprendinius.

Grižkime prie mūsų uždavinio. Kadangi visi dydžiai, įeinantys į rastą $ s $ išraišką, yra teigiami, grynai menamo rodiklio, $ \mathrm{Re}(s)=0 $, galime tikėtis tik tuo atveju, kai išraiškos prieš kvadratines šaknis pavirs nuliais: $ B+k+X_e^2-2X_e Y_e\rightarrow 0 $. Tuo pačiu pošaknio reiškiniai įgis neigiamas vertes, todėl Liapunovo rodikliai $ s=s_1 $ ir $ s=s_2 $ bus grynai menami (teigiamo ir neigiamo ženklo). Įstatykime į minėtą išraišką pusiausvirojo sprendinio \boldmath$sprendinysPusiausvyroje$ tankių $ X_e $ ir $ Y_e $ vertes ir gautą reiškinį prilyginkime nuliui.

\boldmath\begin{eqnarray*}\bigl(B+k+X_e^2-2 Y_e X_e\bigr)^2/.sprendinysPusiausvyroje==0\end{eqnarray*}


Kadangi $ k\neq 0 $, bedimensį šios lygties sprendinį (pastarasis akivaizdus) galime pavadinti stabilumo-nestabilumo ribos parametru (jis skiria $ \mathrm{Re}(s)<0 $ ir $ \mathrm{Re}(s)>0 $ atvejus):

\boldmath\begin{eqnarray*}stabilumoParametras=\Bigl(\frac{B}{k}-\frac{A^2}{k^3}-1\Bigr)\end{eqnarray*}


Kaip matyti iš pilnutinių $ s $ išraiškų, briuseliatorius bus stabilus (sąlyga $ \mathrm{Re}(s)<0 $ išpildyta) tuo atveju, kai stabilumo-nestabilumo ribos parametras tenkins nelygybę $ \mathbf{stabilumoParametras}<0 $. Priešingai, sistema bus nestabili, jei $ \mathbf{stabilumoParametras}>0 $. Įstatę nagrinėtame pavyzdyje naudotas parametrų vertes matome, kad sistema nėra stabili, nes

\boldmath\begin{eqnarray*}stabilumoParametras/.oscParametrai\end{eqnarray*}


Kiti briuseliatoriaus sprendiniai

Ištirkime, kaip elgsis briuseliatorius, esant įvairioms stabilumo-nestabilumo parametro vertėms.

  1. \boldmath$stabilumoParametras$ teigiamas ir artimas nuliui dydis.
    \boldmath$mazasTeigiamas=\bigl\{$
    \boldmath$A$\boldmath$\to$pradinės medžiagos $ A $ tankis
    \boldmath$B$\boldmath$\to$pradinės medžiagos $ B $ tankis
    \boldmath$k$\boldmath$\to$spartos koeficientų santykis
    \boldmath$\bigr\}$

      

    \boldmath\begin{eqnarray*}stabilumoParametras/.mazasTeigiamas\end{eqnarray*}


    \boldmath\begin{eqnarray*}&&sprendinysMT=\mathrm{NDSolve}[\{X'[t]==Y[t] X[t]^2-B X[t]-k X[t]+A,Y'[t]==B X[t]-X[t]^2 Y[t],\\&&\quad X[0]==\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{X_0}$},Y[0]==\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{Y_0}$}\}/.mazasTeigiamas,\{X[t],Y[t]\},\{t,0,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{tend}$}\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys}$}]\end{eqnarray*}


    parinktys =
    $ X_0 $=
    $ Y_0 $=
    tend=

      

    \boldmath$\mathrm{GraphicsRow}[~]$ padeda grafikus surikiuoti vieną šalia kito.
    \boldmath\begin{eqnarray*}&&\mathrm{Show}[\mathrm{GraphicsRow}[\{\\&&\qquad \mathrm{Plot}[\mathrm{Evaluate}[X[t]/.sprendinysMT], \{t, 0, tend\},      \mathrm{AxesLabel}\to \{t, X\}, \mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys1}$}],\\&&\qquad \mathrm{Plot}[\mathrm{Evaluate}[Y[t]/.sprendinysMT], \{t, 0, tend\},      \mathrm{AxesLabel}\to \{t, Y\}, \mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys2}$}]\\&&\}],\mathrm{DisplayFunction}\to \$\mathrm{DisplayFunction},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys3}$}];\end{eqnarray*}

    Pirmojo grafiko piešimo parinktys:
    Antrojo grafiko piešimo parinktys:
    Grafikų grupės piešimo parinktys:

      

  2. Palyginę su brėžiniais skyriaus pradžioje, matome, kad stabilumo-nestabilumo parametrui artėjant prie nulio osciliacijų amplitudė mažėja, o periodas ilgėja.
  3. \boldmath$stabilumoParametras$ neigiamas.
    \boldmath$neigiamas=\bigl\{$
    \boldmath$A$\boldmath$\to$pradinės medžiagos $ A $ tankis
    \boldmath$B$\boldmath$\to$pradinės medžiagos $ B $ tankis
    \boldmath$k$\boldmath$\to$spartos koeficientų santykis
    \boldmath$\bigr\}$

      

    \boldmath\begin{eqnarray*}stabilumoParametras/.neigiamas\end{eqnarray*}


    \boldmath\begin{eqnarray*}&&sprendinysNeig=\mathrm{NDSolve}[\{X'[t]==Y[t] X[t]^2-B X[t]-k X[t]+A,Y'[t]==B X[t]-X[t]^2 Y[t],\\&&\quad X[0]==\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{X_0}$},Y[0]==\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{Y_0}$}\}/.neigiamas,\{X[t],Y[t]\},\{t,0,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{tend}$}\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys}$}]\end{eqnarray*}


    parinktys =
    $ X_0 $=
    $ Y_0 $=
    tend=

      

    \boldmath$\mathrm{GraphicsRow}[~]$ padeda grafikus surikiuoti vieną šalia kito.
    \boldmath\begin{eqnarray*}&&\mathrm{Show}[\mathrm{GraphicsRow}[\{\\&&\qquad \mathrm{Plot}[\mathrm{Evaluate}[X[t]/.sprendinysNeig], \{t, 0, tend\},      \mathrm{AxesLabel}\to \{t, X\}, \mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys1}$}],\\&&\qquad \mathrm{Plot}[\mathrm{Evaluate}[Y[t]/.sprendinysNeig], \{t, 0, tend\},      \mathrm{AxesLabel}\to \{t, Y\}, \mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys2}$}]\\&&\}],\mathrm{DisplayFunction}\to \$\mathrm{DisplayFunction},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys3}$}];\end{eqnarray*}

    Pirmojo grafiko piešimo parinktys:
    Antrojo grafiko piešimo parinktys:
    Grafikų grupės piešimo parinktys:

      

Iš pusiausvirojo sprendinio \boldmath$sprendinysPusiausvyroje$ tankių dar rasime stacionariuosius sprendinius.

\boldmath\begin{eqnarray*}stabilumoParametras/.neigiamas\end{eqnarray*}


Kaip matyti, esant tokioms $ A $, $ B $ ir $ k $ parametrų vertėms, stacionariuosius sprendinius briuseliatorius pasiekia, padaręs vos keletą pereinamųjų osciliacijų. Jei juos pavaizduotume medžiagų tankių $ X-Y $ plokštumoje, uždaro ciklo jau negautume. Sistema relaksuoja iš vienos būsenos į kitą atvira trajektorija:

\boldmath\begin{eqnarray*}&&\mathrm{ParametricPlot}[\\&&\ \ \mathrm{Evaluate}[\{\{X[t],Y[t]\}/.sprendinysNeig],\{t,0,\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{t_{max}}$}\},\mathrm{AxesLabel}\to \{X, Y\},\\&&\qquad \mathrm{Epilog} \to \{\mathrm{PointSize}[0.03], \mathrm{Point}[\{X_0,Y_0\}]\},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys}$}]\end{eqnarray*}

Piešti iki \boldmath$t_{max}$=
Piešimo parinktys:

  

Be briuseliatoriaus egzistuoja ir kitos panašios sistemos, pavyzdžiui \index{osciliatorius!Volterra-Lotka'os}Volterra-Lotka'os osciliatorius, pasiūlytas dar 1920 metais [Lotka20]. Jį nusako reakcijos

\[ \begin{split}A+X\rightarrow 2  X,\notag \\X+Y\rightarrow 2  Y,\\Y\rightarrow E \notag .\end{split}\tag{5} \]
Jei medžiagos $ A $ tankis šiose reakcijose yra išlaikomas pastovus, tada tarpinių medžiagų $ X $ ir $ Y $ tankiai vykstant reakcijai osciliuoja. Volterra detaliai išnagrinėjo šias lygtis ir pritaikė jas grobuonio ir aukos populiacijų kitimui modeliuoti [Volterra31]. Taigi, spartuminių lygčių metodas, kurį šioje knygoje taikėme cheminėms (briuseliatorius) ir fizikinėms (generaciniai-rekombinaciniai procesai) sistemoms modeliuoti taip pat puikiai tinka (mikro)biologinių ir ekologinių sistemų elgesiui aprašyti.

Skyriuje labai trumpai aptartą diferencialinių lygčių sprendinių stabilumo tyrimo metodą nesunkiai galima apibendrinti ir aukštesnės eilės lygtims, t.y. lygtims, kurių jakobianas yra didesnis nei $ 2\times 2 $ matrica.

Atnaujinta 2010.12.03

Literatūra

I. Prigogine and G. Nicolis, "On symmetry breaking instabilities in dissipative systems", J. Chem. Phys. V.46, p.3542-3550 (1967)

I. Prigogine and R. Lefever, "Symmetry breaking instabilities in dissipative systems", J. Chem. Phys. Vol.48, No4, p1695-1700 (1968)

J. J. Tyson, "Some further studies of nonlinear oscillations in chemical systems", J. Chem. Phys. Vol.58, No.9, 3919-3929 (1973)

F. C. Moon, "ChaoticVibrations', Wiley-Interscience Publications, New York, (1987)

A. J. Lotka, J. Am. Chem Soc., Vol.42, p.1595 (1920)

V. Volterra, "Leçon sur la théorie mathématique de la lutte pour la vie", Gauthier-Villars, Paris (1931)

spausdinti