Systèmes dynamiques : modèles biologiques et chimiques
Les réponses

Stéphane Binois
(last modification date: October 10, 2017)
Version pdf de ce document
Version sans bandeaux

1 Réactions en chaîne

Contents

1 Réactions en chaîne
 1.1 Script du fichier chaine.sci
 1.2 Script du fichier chaine.sce

1 Réactions en chaîne

1.1 Script du fichier chaine.sci

Dans ce fichier, nous codons le système dynamique à proprement parler, c’est-à-dire les équations d’évolution des différents éléments présents. On y inclut également les fonctions d’évolution exactes, issues d’une analyse “mathématique” du problème.

approche numérique du système dynamique

  function xdot=chaine(t,x)
    xdot(1)=-k1*x(1);
    xdot(2)=k1*x(1)-k2*x(2);
    xdot(3)=k2*x(2);
  endfunction

évolution exacte du système

  function a=aexact(t)
    a=a0*exp(-k1*t);
  endfunction
  function b=bexact(t)
    b=1/(k2-k1)*(k1*a0*exp(-k1*t)+(k2*b0-k1*(a0+b0))*exp(-k2*t));
  endfunction
  function c=cexact(t)
    c=a0-aexact(t)-bexact(t);
  endfunction

1.2 Script du fichier chaine.sce

Ce fichier est le fichier de commandes, qui appelle les fonctions que l’on vient de définir. Il s’agit principalement de fonctions graphiques, pour tracer les différentes courbes demandées.

  a0=evstr(x_dialog('valeur initiale de A','1'));
  b0=evstr(x_dialog('valeur initiale de B','0'));
  k1=evstr(x_dialog('taux de reaction 1','1'));
  k2=evstr(x_dialog('taux de reaction 2','10'));
  pas=evstr(x_dialog('pas de temps de calcul','0.1'));
  tf=evstr(x_dialog('durée de la simulation','6'));
  X0=[a0;b0;0];
  exec('chaine.sci');
  t=0:pas:tf;

solution exacte

  A=aexact(t);
  B=bexact(t);
  C=cexact(t);
  xset("window",0)
  xbasc()
  leg='A@B@C';
  plot2d([t',t',t'],[A',B',C'],[7,9,3],'121',leg);
  xtitle('Concentration théorique des espèces A B et C','t')

solution ode

  X=ode(X0,0,t,chaine);
  xset("window",1)
  xbasc()
  leg='A@B@C';
  plot2d([t',t',t'],[X(1,:)',X(2,:)',X(3,:)'],[7,9,3],'121',leg);
  xtitle('Concentration calculée par ode des espèces A B et C','t')

Différences

  xset("window",2)
  xbasc()
  leg='A@B@C';
  plot2d([t',t',t'],[X(1,:)'-A',X(2,:)'-B',X(3,:)'-C'],[7,9,3],'121',leg);
  xtitle('Différence entre les concentrations calculées par ode et les concentrations théoriques', ...
         't')