write(%io(2),'reaction equilibree A2+B2<--> 2AB ([LBL], page 470)');
a0=evstr(x_dialog('valeur initiale de A','8'));
b0=evstr(x_dialog('valeur initiale de B','1'));
k=evstr(x_dialog('taux de reaction directe','1'));
r=evstr(x_dialog('taux de reaction inverse','2'));
tpas=evstr(x_dialog('pas de temps de calcul','0.04'));
x0=0;
getf('C:\Mes Documents\Agreg\oraux\equadiff\reac1.sci');
//
// solution exacte obtenue avec MAPLE (puis Copier/coller)
//
deff('x=exact(t)', 'x = -9/7+1/7*sqrt(305)*tanh(1/4*sqrt(305)*t+1/2*log((9+sqrt(305))/(-9+sqrt(305))))')
//
tfin=min(14/abs(a0-b0+0.0001),15);
t=0:tpas:tfin;
//
// solution ode
x=ode(x0,0,t,list(reac1,a0,b0,k,r));
//
// solution Euler
//
npas=int(tfin/tpas);xeul=zeros(npas+1,1);xeul(1)=x0;
for i=1:npas
   xeul(i+1)=xeul(i)+tpas*reac1(tpas,xeul(i),a0,b0,k,r);
end
//
//  solution Rung Kutta 4
//
xrk4=zeros(npas+1,1);xrk4(1)=x0;
for i=2:npas+1
  p1=reac1(tpas,xrk4(i-1),a0,b0,k,r);
  t2=t(i-1)+0.5*tpas;
  y2=xrk4(i-1)+0.5*tpas*p1;
  p2=reac1(tpas,y2,a0,b0,k,r);
  y3=xrk4(i-1)+0.5*tpas*p2;
  p3=reac1(tpas,y3,a0,b0,k,r);
  y4=xrk4(i-1)+tpas*p3;
  p4=reac1(tpas,y4,a0,b0,k,r);
  xrk4(i)=xrk4(i-1)+tpas*((1/6)*p1+(2/6)*p2+(2/6)*p3+(1/6)*p4);
end;
//
xinfty=(a0+b0-sqrt((a0-b0)^2+16*a0*b0*r/k))/(1-4*r/k);
xbasc()
xxx=exact(t);
plot2d([t',t',t',t',t'],[x',xinfty*ones(x'),xeul,xrk4,xxx'],[1,2,3,4,5],'121','AB solveur Scilab@AB infinie @AB Euler@AB RK4@AB exacte');
xtitle('concentration de l''espece AB','t')
//