Méthode de Monte-Carlo cinétique

À chaque étape, le système dispose de plusieurs états accessibles, les taux de transferts vers chacun des états étant supposés connus.

L'algorithme de Monte-Carlo cinétique est utilisé pour modéliser les sauts aléatoires d'un système d'un état initial vers un ensemble d'états finaux possibles, le nombre total N_t de sauts effectués à un instant t étant modélisé par un processus de Poisson. À l'issue de chaque saut, on évalue les conséquences de la transition entre l'état initial et le nouvel état : arrêt de la simulation, nouvelle étape de simulation avec ou sans changement des taux de transfert, etc.

Choix de l'état final : une variable aléatoire est tirée entre 0 et Γ_tot ; la probabilité que le système évolue vers l'état i est proportionnelle à Γ_i.

Pour se familiariser avec l'algorithme, on suppose dans cette section que les taux de transfert restent constants pendant l'intervalle de temps séparant deux sauts successifs, ce qui correspond à un processus de Poisson homogène. L'algorithme est donné par les étapes suivantes :

1. On note t₀ le temps actuel.
2. On forme une liste de l'ensemble des taux de transfert Γ_i du système, i = 1,..., N, et l'on note Γ_tot = ∑ Γ_i le taux de transfert total.
3. Le temps T au bout duquel le système quitte son état initial est donné par la formule T = – ln(U)/Γ_tot, où U est une variable aléatoire uniformément répartie entre 0 et 1.
4. On choisit un deuxième nombre aléatoire U' uniformément réparti entre 0 et Γ_tot qui va servir à tirer au hasard l'état final de la transition (voir la figure ci-contre).
5. Si U' est compris entre 0 et Γ₁, le système transite vers l'état 1 ; entre Γ₁ et Γ₁+Γ₂ le système transite vers l'état 2 ; et si U' est entre S_i-1 et S_i, où S₀ = 0 et ${\displaystyle S_{i}=\textstyle \sum _{j=1}^{i}\Gamma _{j}}$, le système transite vers l'état i.
6. Prendre en compte cette transition en effectuant les calculs correspondants.
7. Changer le temps en t₀+T et recommencer l'étape 1.

Justification de l'étape 3

La variable aléatoire T représente le temps d'attente avant que le système quitte l'état initial. Sa densité de probabilité est ${\displaystyle \rho _{T}(t)=\Gamma _{\text{tot}}\exp(-\Gamma _{\text{tot}}t)}$ pour t ∈ [0, +∞[. La formule utilisée dans l'étape 3 correspond au changement de variable U = exp( – Γ_totT), ce qui donne bien la même densité de probabilité : ${\displaystyle \rho _{T}(t)=\rho _{U}(u)\left|{\frac {du}{dt}}\right|}$ si l'on considère ρ_U(u) = 1 sur ]0, 1].

Justification de l'étape 5

La probabilité que l'état final soit l'état i est p_i = ^Γ_i⁄_Γtot, ce qui est égal à

${\displaystyle p_{i}={\frac {\Gamma _{i}}{\Gamma _{\text{tot}}}}={\frac {S_{i}-S_{i-1}}{\Gamma _{\text{tot}}}}=\int _{S_{i-1}}^{S_{i}}{\frac {1}{\Gamma _{\text{tot}}}}du'=P(S_{i-1}<U'<S_{i})}$

Cet algorithme simule l'évolution du système au cours du temps et n'a pas a priori de condition d'arrêt. En pratique, on arrête la boucle lorsque le temps dépasse une certaine valeur ou lorsque le système a effectué un certain nombre de sauts ; la suite des états i suivis par le système au cours du temps est appelée une trajectoire Monte-Carlo ; on répète ensuite cette boucle de manière à obtenir un nombre important de trajectoires, auxquelles on peut appliquer une analyse statistique. Puisque les trajectoires sont indépendantes les unes des autres, il est facile de paralléliser l'algorithme de Monte-Carlo cinétique.

L'algorithme est essentiellement le même si l'on suppose que les Γ_i dépendent du temps, mais les formules reliant U et U' à T et à l'état final i sont différentes :

• T est donné par la formule ${\displaystyle \int _{t_{0}}^{t_{0}+T}\sum _{i=1}^{N}\Gamma _{i}(t)dt=-\ln u}$, qu'il est nécessaire d'inverser pour trouver T ;
• Les sommes S_i sont à évaluer en t = t₀+T. Ceci provient du fait que jusqu'au temps t₀+T le système est encore dans l'état initial (Chotia 2008). Cette absence de "mémoire" est caractéristique d'un processus de Markov.