next up previous contents index
Next: Génération de nombres aléatoires Up: Équations différentielles Previous: 5.1.2 Van der Pol

5.1.3 Un peu plus d'ode

Dans ce deuxième exemple, nous allons utiliser la primitive ode avec un second membre qui admet un paramètre supplémentaire et nous allons fixer nous même les tolérances pour la gestion du pas de temps du solveur. Voici notre nouvelle équation différentielle (Le Brusselator) :

\begin{displaymath}
\left\{
\begin{array}
{l}
\frac{du_1}{dt} = 2 - (6 + \epsilo...
 ...ac{du_2}{dt} = (5 + \epsilon) u_1 - u_1^2 u_2\end{array}\right.\end{displaymath}

qui admet comme seul point critique $P_{stat} = (2 , (5+\epsilon)/2)$. Lorsque le paramètre $\epsilon$ passe d'une valeur strictement négative à une valeur positive, ce point stationnaire change de nature (de stable il devient instable avec pour $\epsilon = 0$ un phénomène de bifurcation de Hopf). On s'intéresse aux trajectoires avec des conditions initiales voisines de ce point. Voici la fonction calculant ce second membre :
function [f] = Brusselator(t,u,eps)
  //
  f(1) = 2 - (6+eps)*u(1) + u(1)^2*u(2)
  f(2) = (5+eps)*u(1) - u(1)^2*u(2)
Pour faire << passer >> le paramètre supplémentaire, on remplace dans l'appel à ode le nom de la fonction (ici Brusselator) par une liste constituée du nom de la fonction et du ou des paramètres supplémentaires :
  [x] = ode(x0,t0,t,list(MonSecondMembre, par1, par2, ...))
Dans notre cas :
  [x] = ode(x0,t0,t,list(Brusselator, eps))
et l'on procède de même pour tracer le champ avec fchamp.

Pour fixer les tolérances sur l'erreur locale du solveur on rajoute les paramètres rtol et atol, juste avant le nom de la fonction second membre (ou de la liste formée par celui-ci et des paramètres supplémentaires de la fonction). À chaque pas de temps, $t_{k-1} \rightarrow t_k = t_{k-1} + \Delta t_k$, le solveur calcule une estimation de l'erreur locale e (c-a-d l'erreur sur ce pas de temps en partant de la condition initiale v(tk-1) = U(tk-1)) :

\begin{displaymath}
e(t_k) \simeq U(t_k) - \left( \int_{t_{k-1}}^{t_k} f(t,v(t)) dt + U(t_{k-1}) \right)\end{displaymath}

(le deuxième terme étant la solution exacte partant de la solution numérique U(tk-1) obtenue au pas précédent) et compare cette erreur à la tolérance formée par les deux paramètres rtol et atol :

\begin{displaymath}
tol_i = rtol_i*\vert U_i(t_k)\vert + atol_i, \; 1 \le i \le n\end{displaymath}

dans le cas où l'on donne deux vecteurs de longueur n pour ces paramètres et :

\begin{displaymath}
tol_i = rtol*\vert U_i(t_k)\vert + atol, \; 1 \le i \le n\end{displaymath}

si on donne deux scalaires. Si $\vert e_i(t_k)\vert \le tol_i$ pour chaque composante, le pas est accepté et le solveur calcule le nouveau pas de temps de sorte que le critère sur la future erreur ait une certaine chance de se réaliser. Dans le cas contraire, on réintègre à partir de tk-1 avec un nouveau pas de temps plus petit (calculé de sorte que le prochain test sur l'erreur locale soit aussi satisfait avec une forte probabilité). Comme les méthodes mise en jeu sont des méthodes << multipas[*] >> le solveur, en plus du pas de temps variable, joue aussi avec l'ordre de la formule pour obtenir une bonne efficacité informatique... Par défaut les valeurs utilisées sont rtol=10-5 et atol=10-7 (sauf lorsque type selectionne une méthode de Runge Kutta). Remarque importante : le solveur peut trés bien échouer dans l'intégration...

Voici un script possible, la seule fioriture supplémentaire est un marquage du point critique avec un petit carré noir que j'obtiens avec la primitive graphique xfrect :

// Brusselator
eps = -4
P_stat = [2 ; (5+eps)/2];
// limites pour le trace du champ de vecteur
delta_x = 6; delta_y = 4;
x_min = P_stat(1) - delta_x; x_max = P_stat(1) + delta_x;
y_min = P_stat(2) - delta_y; y_max = P_stat(2) + delta_y;
n = 20;
x = linspace(x_min, x_max, n);
y = linspace(y_min, y_max, n);
// 1/ trace du champ de vecteurs
xbasc()
fchamp(list(Brusselator,eps),0,x,y,1,[x_min,y_min,x_max,y_max],"031")
xfrect(P_stat(1)-0.08,P_stat(2)+0.08,0.16,0.16) // pour marquer le point critique
xselect()

// 2/ resolution de l'equation differentielle
m = 500 ; T = 5 ;
rtol = 1.d-09; atol = 1.d-10; // tolerances pour le solveur
t = linspace(0,T,m);
couleurs = [21 2 3 4 5 6 19 28 32 9 13 22 18 21 12 30 27]
num = -1
while %t
  [c_i,c_x,c_y]=xclick();
  if c_i == 0 then
    plot2d(c_x, c_y, -9, "000")  // un petit o pour marquer la C.I.
    u0 = [c_x;c_y];
    [u] = ode(u0, 0, t, rtol, atol, list(Brusselator,eps));
    num = modulo(num+1,length(couleurs));
    plot2d(u(1,:)',u(2,:)',couleurs(num+1),"000")
  elseif c_i == 2 then
    break
  end
end


  
Figure 5.2: Quelques trajectoires dans le plan de phase pour le Brusselator ($\epsilon = -4$)
\begin{figure}
\begin{center}

\includegraphics [width=18cm]{brusselator1}\end{center}\end{figure}


next up previous contents index
Next: Génération de nombres aléatoires Up: Équations différentielles Previous: 5.1.2 Van der Pol
Pincon Bruno
6/23/2000