MOKSLASplius.lt

Monte Carlo metodas

Monte Carlo vardas primena lošimus ir azartą. Fizikams ir matematikams jis asocijuojasi dar ir su kompiuteriniu skaičiavimo metodu, kuriame uždaviniui spręsti naudojamas atsitiktinių skaičių generatorius. Monte Carlo metodas atsirado 1949 metais, kartu su pirmaisiais elektroniniais kompiuteriais, ir siejamas su dviejų matematikų — Johno von Neumanno ir Stanislawo Ulamo — vardais. Šie mokslininkai pirmieji pritaikė tikimybių teoriją sudėtingiems procesams atominiuose reaktoriuose modeliuoti. Fizikoje Monte Carlo metodas dažnai yra naudojamas atsitiktiniams procesams, tokiems kaip Browno dalelių judėjimas skystyje ar elementariųjų dalelių virsmai, modeliuoti. Šiame skyriuje skaitytoją pirmiausia supažindinsime su metodo esme. Po to Monte Carlo metodą pritaikysime eksitonų — kvazidalelių, sudarytų iš elektronų ir skylių, — stochastiniam atsiradimui ir išnykimui aprašyti. Sumodeliuosime lazerio generuojamų eksitonų tankio kitimo laikinę priklausomybę puslaidininkyje ar izoliatoriuje ir ją palyginsime su deterministiniu modeliu.

Nustatymai

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

  

Metodo esmė

Įsivaizduokime, kad iš patrankos apšaudome teritoriją, kurioje tarp žolių yra pasislėpęs netaisyklingos formos ežeras. Sviediniui pataikius į ežerą, virš jo aukštai pakyla iš tolo matomi purslai. Ar galima šaudant patranka išmatuoti ežero plotą, stebint tik virš ežero trykštantį vandenį? Pasirodo, taip. Tik reikia turėti pakankamai sviedinių ir labai blogą šaulį. Suskaičiavę patrankos šūvius ir kiek kartų pasirodė virš ežero purslai, galėsime apytiksliai įvertinti ežero plotą. Jį nustatysime tuo tiksliau, kuo daugiau kartų iššausime. Panagrinėkime konkretų pavyzdį. Monte Carlo metodu apskaičiuosime plotą tarp sinuso funkcijos vieno pusperiodžio kreivės ir abscisių ašies. Kaip žinoma, jis yra tiksliai lygus dviems:

\boldmath\begin{eqnarray*}&&\int_0^\pi \mathrm{Sin}[x]\dd x\end{eqnarray*}

Pradžioje nubraižykime reikiamą sinuso funkcijos kreivės dalį. Po to ją apveskime stačiakampiu, kurio aukštis lygus vienetui, o pagrindas — $ \pi $. Skaičiuojamą plotą nuspalvinkime pilkai.

\boldmath\begin{eqnarray*}&&sinusoPusperiodis=\mathrm{Plot}[\{\sin [x],1\},\{x,0,\pi \},\mathrm{Filling}\to \mathrm{Bottom},\\&&\quad\mathrm{AxesLabel}\to \{x,y\},\mathrm{Ticks}\to \{\{\{0,"0"\},\{ \frac{\pi }{2},"\frac{\pi }{2}\},\{ \pi, "\pi"\}\},\mathrm{Automatic}\},\\&&\qquad\mathrm{DisplayFunction}\to \mathrm{Identity}];\\&&staciakampis=\mathrm{Line}[\{\{ 0, 0\},\{ 0 , 1 \},\{ \pi  , 1 \},\{ \pi  , 0 \},\{ 0 , 0\}\}];\\&&\mathrm{Show}[\{sinusoPusperiodis,\mathrm{Graphics}[staciakampis]\}, \mathrm{DisplayFunction}\to \$\mathrm{DisplayFunction}]\end{eqnarray*}

Piešimo parinktys:

  

Parinktis \boldmath$\mathrm{DisplayFunction}\to\mathrm{Identity}$ sustabdo braižymą, o \boldmath$\mathrm{DisplayFunction}\to\$\mathrm{DisplayFunction}$ vėl ją įjungia. Stačiakampė $ 1\times \pi  $ ploto atskaitos sritis patogi tuo, kad joje galima atsitiktinai ,,barstyti'' taškus, tiesiog generuojant vieną nuo kitos nepriklausančias taško $ x $ ir $ y $ koordinates. Atsitiktinę $ x $ arba $ y $ koordinatę gausime atsitiktinių skaičių generatoriumi — komanda \boldmath$\mathrm{Random}[~]$. Pavyzdžiui, komanda \boldmath$\mathrm{Random}[\mathrm{Real}, \{0, \pi\}]$ generuoja atsitiktinius skaičius intervale tarp nulio ir $ \pi $ ($ x $ koordinatė), o atsitiktinius skaičius tarp nulio ir vieneto duoda komanda \boldmath$\mathrm{RandomReal}[\{0, 1\}]$. Paskutiniuoju atveju galima rašyti paprasčiau, tiesiog \boldmath$\mathrm{Random}[~]$. \boldmath$\mathrm{RandomReal}[\ ]$ komanda yra greitesnė už bendresnę komandą \boldmath$\mathrm{Random}[\ ]$.

\boldmath\begin{eqnarray*}&&\{\mathrm{RandomReal}[\{0,\pi \}],\mathrm{Random}[~]\}\end{eqnarray*}

Ploto skaičiavimo algoritmas Monte Carlo metodu yra labai paprastas. Stačiakampyje atsitiktinių skaičių generatoriumi barstome taškus $ (x,y) $, kurių $ x $ ir $ y $ koordinatės yra atsitiktinės. Jei taškas pakliuvo į ieškomą plotą, t.y. tarp sinuso kreivės ir abscisių ašies, komanda\boldmath$n++$ padidiname skaitliuko \boldmath$n$ turinį vienetu. Priešingu atveju skaitliuko parodymo nekeičiame. Pačią Monte Carlo procedūrą apipavidalinsime kaip modulį.\boldmath$\mathrm{Graphics}^\backprime\mathrm{FilledPlot}^\backprime$.

\boldmath\begin{eqnarray*}&&plotas[taskuSkaicius\_]:=\mathrm{Module}[\{n=0,x,y\},\mathrm{Do}[x=\mathrm{RandomReal}[\{0,\pi \}];\\&&\quad   y=\mathrm{RandomReal}[~];\mathrm{If}[y<\sin [x],n++],\{i,1, taskuSkaicius\}]; sinPlotas=\mathrm{N}[\frac{\pi  n}{taskuSkaicius}]]\end{eqnarray*}
\boldmath$plotas[$ \boldmath$]$

  

Pamėginkite plotą apskaičiuoti nsudodami pavyzdžiui $ 10, 100 $ ir $ 1000 $ taškų. Iš jų matysite, kad apskaičiuoto ploto tikslumas labai priklauso nuo išnaudotų ,,sviedinių'' skaičiaus \boldmath$taskuSkaicius$. Todėl išsiaiškinkime, kaip tikslėja atsakymas, kai generuojame vis daugiau ir daugiau taškų $ (x,y) $. Tam sudarykime du sąrašus — vieną generuotų taškų skaičiui, \boldmath$nList$, o kitą — Monte Carlo metodu surastų plotų \boldmath$plotasList$ vertėms saugoti. Skaičiavimo pradžioje sąrašai yra tušti. Kiekvieno skaičiavimo ciklo metu komanda \boldmath$\mathrm{AppendTo}[~]$ juos papildome naujais duomenimis. Gautą atsakymą — ploto priklausomybę nuo taškų skaičiaus — pavaizduojame grafiškai.

\boldmath$plotasList=\{\};i=0;kiekTasku=$
\boldmath\begin{eqnarray*}&&\mathrm{Do}[\mathrm{If}[\mathrm{Random}[]<\sin [\mathrm{Random}[\mathrm{Real},\{0,\pi \}]],i+=1];\\&&\mathrm{AppendTo}[plotasList,\mathrm{N}[\frac{i \pi }{n}]],\{n,1,kiekTasku\}]\end{eqnarray*}
Piešimo parinktys:

  

Konkreti tokiu būdu surasto ploto vertė priklauso nuo komanda \boldmath$\mathrm{Random}[~]$ sugeneruotų atsitiktinių skaičių, o kiekvienas naujas skaičiavimas duos vis kitokią plotas-taškai kreivę. Pavaizduota vieno skaičiavimo kreivė vadinama realizacija. Statistiniais metodais parodoma, kad jei kiekvienoje iš daugelio realizacijų buvo sugeneruota $ n $ atsitiktinių taškų, Monte Carlo integravimo metodo santykinė paklaida yra $ \Delta n/n = 1/\sqrt{n} $. Pasinaudojus formule nesunku įvertinti, kad, pavyzdžiui, sugeneravus $ 2000 $ taškų, santykinė taip apskaičiuoto ploto paklaida turėtų būti apie $ 0{,}022 $. Tokiu būdu Monte Carlo integravimo tikslumas auga daug lėčiau negu atsitiktinai barstomų stačiakampyje taškų skaičius $ n $. Antra vertus, skaičiavimo laikas auga proporcingai $ n $. Todėl praktikoje parenkant optimalų ,,šūvių'' skaičių tenka ieškoti kompromiso tarp pageidaujamo tikslumo ir priimtino skaičiavimo laiko. Be to, kaip matyti iš pirmojo brėžinio, stačiakampio, o bendresniu atveju — daugiamačio atskaitos stačiakampio ,,tūris'' turi būti kiek galima mažesnis. Optimaliu atveju jis turėtų būti vos didesnis už ieškomąjį: tada bus mažiau ir tuščių, nepakliūnančių į skaičiuojamąjį tūrį taškų.

Kam reikalingas toks keistas skaičiavimo algoritmas, kai yra žinoma visa eilė puikių deterministinių metodų, kurių tikslumas bei skaičiavimo laikas yra gerai žinomi? Be abejo, ,,gerai'' besielgiančių, vieno kintamojo tolydžių funkcijų atveju, kokia ir yra čia aptarta \boldmath$\mathrm{Sin}[x]$ funkcija, tikimybių teorija paremto integravimo metodo niekas ir nenaudoja. Kvantų mechanikoje, signalų analizėje ir daugelyje kitų sričių susiduriama su labai ,,blogomis'' funkcijomis, priklausančiomis nuo daugelio kintamųjų ir turinčiomis daugybę smailių. Štai tada ir atsiskleidžia Monte Carlo metodo pranašumas, kuriam iš tiesų visiškai nesvarbi nei integruojamos funkcijos forma, nei integruojamų dimensijų skaičius. Deterministiniai gi metodai vienu ar kitu būdu visada atsižvelgia į funkcijos kitimo pobūdį, todėl dažnai siaurų smailių (ypač pavienių) jie gali ir ,,nepastebėti''. Integruojant daugelio kintamųjų funkcijas žinomi deterministiniai algoritmai veikia dar ir žymiai lėčiau.

Eksitonų susidarymo ir išnykimo kinetika

Atsitiktinių skaičių generatoriai plačiai taikomi stochastiniuose modeliuose. Tokiais vadinami modeliai, turintys atsitiktinę struktūrą arba nagrinėjantys atsitiktinius laike procesus. Čia trumpai aptarsime nesudėtingą modelį, aprašantį atsitiktinį eksitonų susidarymą ir jų išnykimą. Eksitonu vadinama kvazidalelė arba sužadinimas stebimas puslaidininkiuose bei dielektrikuose. Jį sudaro laidumo juostos elektronas ir valentinės juostos skylė. Kadangi elektronas $ e $ turi neigiamą krūvį, o skylė $ h $ — teigiamą, tai dalelės vieną kitą traukia ir gali susidaryti surišta būsena, eksitonas. Eksitonas labai primena tikrų dalelių sistemą — pozitronijų. Pastarąjį, kaip žinoma, galima įsivaizduoti kaip ,,molekulę'', kurioje elektronas ir pozitronas sukasi apie bendrą abiejų dalelių masės centrą. Eksitono susidarymą galima pavaizduoti chemine reakcija $ e+h\rightarrow X $, kurioje raidė $ X $ žymi susidariusią kvazidalelę arba eksitoną. Eksitono gyvavimo trukmę lemia keletas priešingo veikimo procesų.

  1. Dėl elektrostatinės traukos gali susidaryti eksitonas: $ e+h\rightarrow X $.
  2. Veikiant gardelės virpesiams, eksitonas gali suirti: $ X\rightarrow e+h $.
  3. Kaip ir pozitronijus, eksitonas gali anihiliuoti, dėl ko vietoje kvazidalelės atsiranda energijos $ \hbar \omega  $ ($ \hbar $ — Plancko pastovioji, $ \omega $ — kvanto dažnis) šviesos kvantas: $ X\rightarrow \hbar \omega  $.

Kiekvienos iš trijų paminėtų reakcijų efektyvumą nusako dalelės atsiradimo ir išnykimo sparta (rate). Šiomis spartomis galima išreikšti elektronų, skylučių ir eksitonų tankių kitimą laike. Reakcijų spartos yra proporcingos dalyvaujančių reakcijoje dalelių skaičiui ar jų tankiui. Laikysime, kad bet kuriuo laiko momentu yra vienodas elektronų $ n_e $ ir skylučių $ n_h $ tankis, kurį žymėsime $ n_{eh}=n_e=n_h $. Jei $ n_X $ yra eksitonų tankis, tada trijų aprašytų reakcijų spartos užrašomos taip:

\[ \begin{array}{rl}r(1)&=\gamma  n_{eh}^2,\notag  \\r(2)&=n_X/T_{eh}, \tag{1}\\r(3)&=n_X/T_\nu .\notag\end{array} \]
Proporcingumo koeficientai lygtyse, būtent, $ \gamma  $ — eksitono susidarymo koeficientas, $ T_{eh} $ — eksitono suirimo į $ e-h $ porą vidutinė gyvavimo trukmė ir $ T_\nu $ — eksitono anihiliacijos vidutinė trukmė, yra apskaičiuojami arba teoriškai (pavyzdžiui, kvantinės mechanikos metodais), arba išmatuojami eksperimente. Kiekvieną iš aptartųjų spartų Mathematica kalboje nusakysime tokiu būdu:
\boldmath\begin{eqnarray*}&&r[1]:=\gamma  n_{eh}^2;\\&&r[2]:=\frac{n_X}{T_{eh}};\\&&r[3]:=\frac{n_X}{T_{\nu }}\end{eqnarray*}

$ i- $osios reakcijos tikimybė $ p(i) $ yra proporcinga jos spartos ir visų reakcijų spartų sumos santykiui:

\boldmath\begin{eqnarray*}&&p[i\_]:=\frac{r[i]}{\sum _{k=1}^3 r[k]}\end{eqnarray*}

Kadangi tikimybės yra proporcingos spartoms, o pastarosios priklauso nuo dalyvaujančių reakcijose dalelių skaičiaus, tai vykstant reakcijai tikimybės gali keistis. Pavyzdžiui, eksitono anihiliacijos tikimybė duotu laiko momentu priklauso nuo elektronų, skylių ir pačių eksitonų skaičiaus tokiu būdu:

\boldmath$p[$\boldmath$]$

  

Į eksitono susidarymą ir jo suirimą žiūrėsime kaip į stochastinį (atsitiktinį) procesą. Stochastiniame procese elementarios reakcijos vyksta atsitiktinai su anksčiau aprašytomis tikimybėmis $ p(1) $, $ p(2) $ ir $ p(3) $. Tam tikru laiko momentu jos, pavyzdžiui, gali įgyti tokias vertes: $ p(1)=0{,}2 $; $ p(2)=0{,}3 $; $ p(3)=0{,}5 $. Algoritme konkrečią reakciją išrinksime, pasinaudodami didėjančiomis iki vieneto tikimybių sumomis. Pateiktoms tikimybėms šias sumas galima pasirinkti, pavyzdžiui, tokias:

\[ \begin{array}{rl}p(1)=&0{,}2 \notag \\p(1)+p(2)=&0{,}5 \tag{2}\\p(1)+p(2)+p(3)=&1{,}0\notag\end{array} \]
Kaip gi išrenkama konkreti reakcija? Pirmasis stochastinio modeliavimo principas remiasi tokiu reakcijos išrinkimo algoritmu. Atsitiktinių skaičių generatoriumi generuojame atsitiktinį skaičių $ p_1 $, telpantį tarp $ 0 $ ir $ 1 $. Jei sugeneruotas skaičius tenkina nelygybę $ 0<p_1\le p(1) $, laikome, kad įvyko pirmoji reakcija. Jei jis yra iš intervalo $ p(1)< p_1<p(1)+p(2) $, tada laikome, kad įvyko antroji reakcija, o jei $ p(1)+p(2)\le p_1<p(1)+p(2)+p(3) $ — trečioji reakcija.
\boldmath\begin{eqnarray*}&&kuriReakcija:=\mathrm{Switch}[\mathrm{Random}[],\\&&\qquad \_?(#1\leq p[1]\&), 1,\quad \_?(p[1]<#1<p[1]+p[2]\&),2,\quad \_?(#1\geq p[1]+p[2]\&),3]\end{eqnarray*}

Šiame algoritme priklausomai nuo konkrečios tikimybės $ p(i) $ reikšmės, kurią sugeneravo \boldmath$\mathrm{Random}[~]$, komanda \boldmath$\mathrm{Switch}[~]$ atrenka turinčios įvykti reakcijos numerį $ i=1 $, $ 2 $ arba $ 3 $. Algoritme vietoje \boldmath$\#$ įstatomas sugeneruotas atsitiktinis skaičius ir iš eilės patikrinama, kuri iš grynosiomis funkcijomis nusakytų sąlygų \boldmath$\_{}?(\#{}\leq p[1]\&{})$ ar \boldmath$\_{}?(p[1]<\#{}<(p[1]+p[2])\&{})$, ar \boldmath$\_{}?(\#{}\geq  p[1]+p[2]\&{})$ yra patenkinama. Nuo to priklauso, koks numeris ($ 1 $, $ 2 $, ar $ 3 $) bus priskirtas identifikatoriui \boldmath$kuriReakcija$. Klaustukas prieš grynąsias funkcijas yra komandos \boldmath$\mathrm{PatternTest}[~]$ santrumpa. Brūkšnelis prieš klaustuką rodo, kad be apribojimų tikrinamas kiekvienas \boldmath$\mathrm{Random}[~]$ sugeneruotas skaičius.

Po to, kai stochastiškai nustatyta, kuri iš reakcijų turi įvykti, pasinaudoję reakcijų taisyklių apibrėžimais, keičiame dalelių skaičių. Mūsų atveju po reakcijos dalelių skaičius gali padidėti, sumažėti vienetu arba iš viso nepasikeisti. Šios taisyklės akivaizdžiai seka iš reakcijų apibrėžimų, būtent:

  • jei išrinkta 1-a reakcija, tada $ \Delta n_{eh}=-1, \quad \Delta n_X=+1 $,
  • jei išrinkta 2-a reakcija, tada $ \Delta n_{eh}=+1, \quad \Delta n_X=-1 $,
  • jei išrinkta 3-a reakcija, tada $ \Delta  n_{eh}=0,\hphantom{-}\quad \Delta  n_X=-1 $.

Pažymėję reakcijos numerį \boldmath$rn$, visas sąlygas užrašysime tokiu pavidalu:

\boldmath\begin{eqnarray*}&&daleles[rn\_]:=\mathrm{Switch}[rn,\\&&\qquad 1,\mathrm{If}[n_{eh}>0,n_{eh}--;n_X++],\quad       2,\mathrm{If}[n_X>0,n_{eh}++;n_X--],\quad 3,\mathrm{If}[n_X>0,n_X--]]\end{eqnarray*}

Čia sąlygos operatoriumi \boldmath$\mathrm{If}[~]$ apsidraudžiame nuo nefizikinio atvejo, kad dalelių skaičius nepasidarytų neigiamas.

Antrasis stochastinio modeliavimo principas teigia, kad laiko tarpas tarp elementarių reakcijų taip pat yra atsitiktinis dydis. Todėl laiko tarpą $ \tau $ tarp dviejų viena po kitos sekančių reakcijų vėl nulems atsitiktinių skaičių generatoriumi sugeneruotas skaičius. Bendru atveju $ \tau $ galėtų būti bet kokio dydžio, $ 0<\tau <\infty $. Tačiau patyrimas rodo, kad mažos $ \tau  $ vertės turi didesnę tikimybę už dideles. Jei procesai nėra tarpusavyje priklausomi, $ \tau  $ pasiskirstymas turi eksponentinį pobūdį. Šis teiginys gerai dera su eksperimentu. Todėl sugeneravę kitą atsitiktinį skaičių $ p_2 $, kuris tenkina nelygybę $ 0<p_2<1 $, laiko tarpą $ \tau  $ tarp dviejų reakcijų apskaičiuosime iš pasiskirstymo formulės:

\[ p_2=\exp\big(-\smash{\sum_i} \,r(i)\,\tau \big).\tag{3} \]
Čia $ \sum_i r(i) $ yra pilnutinė visų reakcijų sparta. Taigi, trumpi tarpai tarp elementarių reakcijų pasitaikys eksponentiškai dažniau nei ilgi. Paėmę užrašytos formulės abiejų pusių logaritmą, laiko tarpą nuo vienos elementarios reakcijos iki kitos elementarios reakcijos skaičiuosime tokiu algoritmu:
\boldmath\begin{eqnarray*}&&tarpas[p2\_]:=-\frac{\log [p2]}{\sum _{k=1}^3 r[k]}\end{eqnarray*}

Sukurkime tris sąrašus, kuriuose atitinkamai kaupsime $ e-h $ porų skaičių, eksitonų skaičių ir laiko tarpą $ \tau $, praėjusį po kiekvienos elementarios reakcijos. Pradžioje sąrašai yra tušti.

\boldmath\begin{eqnarray*}&&nehList=\{\};\quad nXList=\{\};\quad laikasList=\{\};\end{eqnarray*}

Užduosime tokias koeficientų ir pradinių dalelių skaičių reikšmes:

\boldmath$\gamma$\boldmath$=$eksitono susidarymo koeficientas
\boldmath$T_{eh}$\boldmath$=$eksitono vidutinė gyvavimo trukmė prieš suyrant į $ e-h $ porą
\boldmath$T_{\gamma}$\boldmath$=$eksitono vidutinė gyvavimo trukmė prieš virstant $ \gamma $ kvantu
\boldmath$n0_{eh}$\boldmath$=$ Pradiniu laiko momentu susidariusių $ e-h $ porų skaičius
\boldmath$n0_{X}$\boldmath$=$Eksitonų skaičius pradiniu laiko momentu
\boldmath$t0$\boldmath$=$Pradinis laiko momentas
\boldmath$tmax$\boldmath$=$Galinis laiko momentas

  

Čia \boldmath$t0$ ir \boldmath$tmax$ žymi pradinį ir galinį modeliavimo laiką. Laikome, kad pradiniu laiko momentu, pavyzdžiui, tiriamąjį bandinį apšvietus mažos trukmės lazerio impulsu, jame susikuria $ 500 $ $ e-h $ porų. Visus aprašytus algoritmus surenkame į vieną ir inicializuojame pradines sąlygas.

\boldmath\begin{eqnarray*}&&n_{eh}=n0_{eh};\quad n_X=n0_X;\quad t=t0;\\&&\mathrm{While}[t<tmax,daleles[kuriReakcija];\\&&\quad    t+=tarpas[\mathrm{Random}[]];\mathrm{AppendTo}[nehList,n_{eh}];\mathrm{AppendTo}[nXList,n_X];\mathrm{AppendTo}[laikasList,t]]\end{eqnarray*}

Sulig kiekvienu žingsniu komanda \boldmath$\mathrm{AppendTo}[~]$ sąrašus papildo naujais elementais. Kiekvieno ciklo metu sugeneruojami du atsitiktiniai skaičiai: sugeneravus pirmąjį išrenkama viena iš trijų reakcijų (\boldmath$\mathrm{Random}[~]$ yra paslėpta komandoje \boldmath$kuriReakcija$), o sugeneravus antrąjį — nustatomas atsitiktinis laiko tarpas $ \tau  $ iki kitos elementarios reakcijos.Taigi, laikas modelyje ,,teka'' nenuspėjamais diskretiškais šuoliukais.

\boldmath\begin{eqnarray*}&&nXt=\mathrm{Transpose}[\{laikasList,nXList\}];\\&&neht=\mathrm{Transpose}[\{laikasList,nehList\}];\end{eqnarray*}

Vizualizavę apskaičiuotus sąrašus, matysime, kaip eksitonų bei $ e-h $ porų skaičius priklauso nuo laiko.

\boldmath\begin{eqnarray*}&&\\&&nXtStoch=\mathrm{ListPlot}[nXt,      \mathrm{AxesLabel}\to \{t, \mu s,n_X\},\mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys1}$}];\\&&nehtStoch=\mathrm{ListPlot}[neht,\mathrm{AxesLabel}\to \{t, \mu s,n_{eh}\},          \mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 6pt{}}\smash{parinktys2}$}];\\&&\mathrm{Show}[\mathrm{GraphicsArray}[\{nXtStoch,nehtStoch\}] ,    \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 $ e-h $ porų, kurias laiko momentu $ t=0 $ sukūrė trumpas šviesos impulsas, skaičius mažėja laikui bėgant. Tuo tarpu eksitonų skaičius, priešingai, iš pradžių auga, nes didelis $ e-h $ porų skaičius skatina jų susidarymą. Tačiau ilgainiui mažėjant $ e-h $ porų skaičiui, eksitonų skaičius taip pat ima mažėti.

Palyginimas su deterministiniu modeliu

Kadangi pasirinktas dalelių skaičius nėra didelis, brėžinyje parodytų kreivių forma šiek tiek keisis su kiekvienu kompiuteriniu eksperimentu. Didinant dalelių skaičių, fliuktuacijos mažės. Riboje, kai turime be galo daug dalelių, stochastinius rezultatus galima gauti ir remiantis deterministiniu modeliu. Deterministiniai reakcijų modeliai yra aprašomi pirmos eilės diferencialinėmis lygtimis. Fizikoje tokios lygtis yra vadinamos spartuminėmis (rate equations). Kiekviena spartuminė lygtis aprašo vienos rūšies dalelių tankio kitimo kinetiką. Kairėje spartuminių lygčių pusėje rašome konkrečios dalelės tankio išvestinę laiko atžvilgiu, o dešinėje išvardiname visas spartas, kuriose minėta dalelė dalyvauja. Jei reakcija mažina nagrinėjamų dalelių skaičių, spartą rašome su minuso, o jei didina — su pliuso ženklu. Pavyzdžiui, mūsų uždavinyje eksitonų kitimo kinetika aprašoma tokia spartumine lygtimi:

\[ \frac{\mathrm{d} n_X}{\mathrm{d} t}=r(1)-r(2)-r(3).\tag{4} \]
Elektronams ir skylėms aprašyti pakanka vienos diferencialinės lygties, nes pagal padarytą prielaidą jų skaičius bet kuriuo laiko momentu yra vienodas ir nusakomas ta pačia diferencialine lygtimi
\[ \frac{\mathrm{d} {n_{eh}}}{\mathrm{d} t}=r(2)-r(1).\tag{5} \]
Mathematica kalboje lygtys (4), (5) užrašomos tokiu būdu:
\boldmath\begin{eqnarray*}&&\mathrm{Clear}[t];\\&&lygtis1=\bigl( nX'[t]==\gamma * (neh[t])^2-\frac{nX[t]}{T_{eh}}-\frac{nX[t]}{T_{\nu }}\bigr);\\&&lygtis2=\bigl( neh'[t]==\frac{nX[t]}{T_{eh}}-\gamma * (neh[t])^2\bigr);\end{eqnarray*}

Įvedę pradines sąlygas, skaitinius diferencialinių lygčių sprendinius rasime \boldmath$\mathrm{NDSolve}[~]$ komanda.

\boldmath\begin{eqnarray*}&&sol=\mathrm{NDSolve}[\{lygtis1,lygtis2, neh[0]==n0_{eh},nX[0]==n0_X\},\{neh,nX\},\{t,t0,tmax\}]\end{eqnarray*}

Palyginimui gautus deterministinius sprendinius pavaizduosime kartu su stochastiniais.

\boldmath\begin{eqnarray*}&&\\&&nXtDeterm=\mathrm{Plot}[\mathrm{Evaluate}[nX[t]t/.sol], \{t, t0, tmax\},      \mathrm{AxesLabel}\to \{t, \mu s,n_X\},\\[1pt]&&\quad \mathrm{PlotStyle}\to \{\mathrm{RGBColor}[1,0,0],\mathrm{RGBColor}[0,0,0]\},\mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 5pt{}}\smash{parinktys1}$}];\\&&nehtStoch=\mathrm{Plot}[\mathrm{Evaluate}[neh[t]/.sol], \{t, t0, tmax\},      \mathrm{AxesLabel}\to \{t, \mu s,n_X\}, \\[2pt]&&\quad \mathrm{PlotStyle}\to \{\mathrm{RGBColor}[1,0,0],\mathrm{RGBColor}[0,0,0]\},\mathrm{DisplayFunction}\to \mathrm{Identity},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 5pt{}}\smash{parinktys2}$}];\\&&\mathrm{Show}[\mathrm{GraphicsArray}[\{\mathrm{Show}[\{nXtStoch,nXtDeterm\}],\mathrm{Show}[\{nehtStoch,nehtDeterm\}]\}],\\[2pt]&&\quad   \mathrm{DisplayFunction}\to \$\mathrm{DisplayFunction},\fcolorbox[rgb]{1,0,0}{1,1,1}{$\vphantom{\vbox to 4pt{}}\smash{parinktys3}$}];\end{eqnarray*}

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

  

Šiame eksperimente aprašyti metodai taikomi ne tik fizikoje, bet ir chemijoje, biologijoje, biofizikoje, farmakologijoje, — cheminėms reakcijoms, katalizei, biocheminiams ciklams ir t.t. modeliuoti. Tam tikslui yra sudarytos įvairių medžiagų spartos koeficientų lentelės, todėl dominančią cheminę ar kitokią reakciją pradžioje verta sumodeliuoti kompiuteriu, o tik po to ją bandyti realizuoti laboratorijoje. Taip yra ir pigiau, ir saugiau. Ypač, jei reakcijose išsiskiria žmonėms ir aplinkai pavojingos medžiagos. Žinant, ko galima tikėtis iš reakcijos, ją galima realizuoti greičiau, na, o eksperimentas galutinai parodys, kiek jūsų modelis atspindi realybę. Modeliavimas Monte Carlo metodu taip pat leidžia nagrinėti fliuktuacijas, nes reakcijoje dalyvaujančių dalelių skaičius visada yra baigtinis. Tokios fliuktuacijos gerai matomos ir čia nupieštuose brėžiniuose. Kitame eksperimente, pasinaudoję spartuminėmis lygtimis, panagrinėsime briuseliatorių — klasikinį nestabilios cheminės reakcijos pavyzdį.

Matematinį įvadą apie Monte Carlo metodą skaitytojas ras knygoje [Sobol68]. Jo taikymas krūvininkų pernašoje pateiktas apžvalginiame staipsnyje [Jacoboni83], o taikymams statistinėje fizikoje yra skirta visa speciali K.mBinderio redaguojama knygų serija [Binder87,Binder88].

Atnaujinta 2010.12.05

Literatūra

И.М. Собол "Mетод Монте-Карло" (серия популярных лекции по математике), Москва, Наука (1968).

C. Jacoboni, L. Reggiani, "The Monte Carlo method for the solution of charge transport in semiconductors with application to covalent materials", Rev. Mod. Phys. Vol. 55, No3., p.645-705 (1983)

K. Binder, " Monte Carlo methods in statistical physics", Springer, Berlin (1987)

K. Binder, D. W. Heermann, " Monte Carlo simulation in statistical physics", Springer, Berlin (1988)

spausdinti