Méthode de Runge–Kutta pour les EDS

Considérons la diffusion Itō ${\displaystyle X}$ satisfaisant l'équation différentielle stochastique Itō suivante

${\displaystyle {{d}X_{t}}=a(X_{t})\,{d}t+b(X_{t})\,{d}W_{t},}$

avec comme condition initiale ${\displaystyle X_{0}=x_{0}}$, où ${\displaystyle W_{t}}$ représente un processus de Wiener, et supposons que nous souhaitons résoudre cette EDS sur un certain intervalle de temps ${\displaystyle [0,T]}$ . Alors, l' approximation Runge – Kutta de la vraie solution ${\displaystyle X}$ est la chaîne de Markov ${\displaystyle Y}$ défini comme suit[1]:

• découper l'intervalle ${\displaystyle [0,T]}$ dans ${\displaystyle N}$ sous-intervalles de largeur ${\displaystyle \delta =T/N>0}$ :

${\displaystyle 0=\tau _{0}<\tau _{1}<\dots <\tau _{N}=T;}$

• soit ${\displaystyle Y_{0}:=x_{0}}$ ;
• calculer récursivement ${\displaystyle Y_{n}}$ pour ${\displaystyle 1\leq n\leq N}$ par

${\displaystyle Y_{n+1}:=Y_{n}+a(Y_{n})\delta +b(Y_{n})\Delta W_{n}+{\frac {1}{2}}\left(b({\hat {\Upsilon }}_{n})-b(Y_{n})\right)\left((\Delta W_{n})^{2}-\delta \right)\delta ^{-1/2},}$

où ${\displaystyle \Delta W_{n}=W_{\tau _{n+1}}-W_{\tau _{n}}}$ et ${\displaystyle {\hat {\Upsilon }}_{n}=Y_{n}+a(Y_{n})\delta +b(Y_{n})\delta ^{1/2}.}$ Les variables aléatoires ${\displaystyle \Delta W_{n}}$ sont des variables aléatoires normales indépendantes et distribuées de manière identique avec une espérance nulle et une variance ${\displaystyle \delta }$ .

Ce schéma a un ordre fort 1, ce qui signifie que l'erreur d'approximation de la solution réelle avec une échelle de temps fixe est proportionnelle au pas de temps ${\displaystyle \delta }$ . Il y a également un ordre faible 1, ce qui signifie que l'erreur sur les statistiques de la solution évolue avec le pas de temps ${\displaystyle \delta }$ . Voir les références pour des déclarations complètes et exactes.

Les fonctions ${\displaystyle a}$ et ${\displaystyle b}$ peut varier dans le temps sans aucun problème. La méthode peut être généralisée au cas de plusieurs équations couplées ; le principe est le même mais les équations s'allongent.

Un nouveau schéma Runge-Kutta également d'ordre fort 1 se réduit directement au schéma Euler amélioré pour les ODE déterministes[2]. Considérons le processus stochastique vectoriel ${\displaystyle {\vec {X}}(t)\in \mathbb {R} ^{n}}$ qui satisfait une EDS îto :

${\displaystyle d{\vec {X}}={\vec {a}}(t,{\vec {X}})\,dt+{\vec {b}}(t,{\vec {X}})\,dW,}$

où ${\displaystyle {\vec {a}}}$ est la dérive et ${\displaystyle {\vec {b}}}$ la diffusion, sont des fonctions suffisamment lisses. Compte tenu du pas de temps ${\displaystyle h}$, et étant donné que ${\displaystyle {\vec {X}}(t_{k})={\vec {X}}_{k}}$, on estime ${\displaystyle {\vec {X}}(t_{k+1})}$ par ${\displaystyle {\vec {X}}_{k+1}}$ pour le temps ${\displaystyle t_{k+1}=t_{k}+h}$ via

${\displaystyle {\begin{array}{rl}&{\vec {K}}_{1}=h{\vec {a}}(t_{k},{\vec {X}}_{k})+(\Delta W_{k}-S_{k}{\sqrt {h}}){\vec {b}}(t_{k},{\vec {X}}_{k}),\\&{\vec {K}}_{2}=h{\vec {a}}(t_{k+1},{\vec {X}}_{k}+{\vec {K}}_{1})+(\Delta W_{k}+S_{k}{\sqrt {h}}){\vec {b}}(t_{k+1},{\vec {X}}_{k}+{\vec {K}}_{1}),\\&{\vec {X}}_{k+1}={\vec {X}}_{k}+{\frac {1}{2}}({\vec {K}}_{1}+{\vec {K}}_{2}),\end{array}}}$

• où ${\displaystyle \Delta W_{k}={\sqrt {h}}Z_{k}}$ avec ${\displaystyle Z_{k}}$ la loi normale centrée réduite ${\displaystyle Z_{k}\sim N(0,1)}$ ;
• et où ${\displaystyle S_{k}=\pm 1}$, chaque valeur étant choisie avec probabilité ${\displaystyle 1/2}$ .

Ce qui précède décrit un seul pas de temps. Il faut répéter ce processus afin d'intégrer la EDS sur tout l'intervalle de temps.

Des schémas d'ordre supérieur existent également, mais deviennent de plus en plus complexes. Rößler a développé de nombreuses solutions pour les EDS de type Ito[3] - [4], tandis que Komori a développé des solutions pour les EDS de Stratonovich[5] - [6] - [7]. Rackauckas a étendu ces schémas pour permettre un pas de temps adaptatif via l'échantillonnage de rejet avec mémoire (RSwM), ce qui entraîne des augmentations de l'efficacité de plusieurs ordres de grandeur dans des modèles biologiques[8] , ainsi qu'une optimisation des coefficients pour une stabilité améliorée.