// reacteur biologique
/////////////////////////
// 1er schema: Euler + trapezes
//
kappa=0.5
//
uinit=0.1;
//
N=evstr(x_dialog('nombre de points','100'))
T=evstr(x_dialog('temps final','20'))
//
h=T/N;t=0:h:T;
//
u0=uinit;u=[uinit];
S=0;
for i=1:N
   uold=u0; 
   u0=u0+h*u0/kappa*(1-u0-S)
   S=S+h/2*(uold+u0)
   u=[u,u0];
end
//
figure(1)
clf()
plot2d(t,u,1)
////////////////////////////
// 2 eme schema: ode et système differentiel
////////////////////////////
//
deff('xdot=f(t,x)','xdot=[x(2),x(2)*(1-x(2)-x(1))/kappa]')
//
xsol=ode([0;uinit],0,t,f)
plot2d(t,xsol(2,:),2)
legends(['solution schema EDO','solution systeme ode'],[1,2])
//
////////////////////////////////////////`
// portrait de phase
///////////////////////////////////////
//
figure(2)
xf=0:0.1:2;yf=0:0.1:2;
clf()
fchamp(f,0,xf,yf)
for i=1:10
   uinitb=2*i/10;
   xsol=ode([0;uinitb],0,t,f)
   plot2d(xsol(1,:),xsol(2,:),i)   
end
plot([0,1],[1,0],'b-')
xtitle('portait de phase')
//////////////////////
// calcul de la production de biogaz
//////////////////////
//
disp('production totale (méthode 1: avec le schema EDO)');disp(S)
disp('production totale (méthode 2: avec la formule (13))');
deff('y=phi(x)','y=(1+kappa-x)-(1+kappa-uinit)*exp(-x/kappa)');
prod=fsolve(1,phi);
disp(prod)
