
///// Resolution du systeme proie predateur avec diffusion
// systeme de 2 eq: U proies, V predateurs
// Ut= f1(U)  + D*Uxx et Vt=f2(V) + Vxx
///   sur [a,b]
////  par differences finies. Schema Euler implicite pour diffusion
///  explicite pour terme NL
///  Guillaume Chiavassa, mars 2017

clear;


// matrice tri diagonale du laplacien en dim 1, differences finies centrees
// CL nulles au bord  coeff*[-2 1; 1 -2 1;....., 1 -2] de taille M*M

function [B]=laplacien(coef,M)

 B=zeros(M,M);
 // diagonale
 for i=1:M,
  B(i,i)=-2*coef;
 end;

 //sur diagonale
 for i=1:M-1,
  B(i,i+1)= coef;
 end;
 // sous diagonale
 for i=2:M,
  B(i,i-1)=coef;
 end;

endfunction;


//termes nonlineaires, r,K scalaires, U,V vecteurs colonnes

// f1= r1*U*(1-U-K1*V)
function [sol]= f1(r,K,U,V)
  sol=r*U.*(1-U-K*V);
endfunction;

//f2=r2*V*(1-V-K2*U)
function [sol]= f2(r,K,U,V)
 sol=r*V.*(1-V-K*U);
endfunction;


// condition initiale
// gaussienne pour chaque espece, centree en x0

 function U=cond_init(x0,x)
  U=0.5*exp(-(x-x0).^2/2)' 
 endfunction;


//----------------------------------------------------------------------
// programme principal

// parametres

 r1=0.6; K1=0.6;  
 r2=0.6; K2=0.4;
 D=.5;

 a=-75; b=75; //domaine

N=150;
dx=(b-a)/N;
x=[a:dx:b];

//temps final

T=50.;

dt=0.02;


// trace la condition initiale
xset("window",0)
clf();
U0=cond_init(a+(b-a)/2,x);
plot2d(x,U0,1);
V0=cond_init(a+(b-a)/2,x);
plot2d(x,V0,2);
title('Condition initiale, proies (noir) predateurs (bleu)');

// les variables


B=zeros(N-1,N-1);
NL=zeros(N-1,1);

//-------  remplissage de B

 coef=dt/dx^2;

 B1=laplacien(D*coef,N-1);
 B2=laplacien(coef,N-1);
 
//calcul de l'inverse pour l'implicite
  
  M1=eye(N-1,N-1)-B1;
  M2=eye(N-1,N-1)-B2;

// initialisation de VN

   UN=U0(2:N);
   VN=V0(2:N);

// BOUCLE EN TEMPS
clf()



t=0;
iter=1;

while t<T,
  
    // terme NonLin

       NL1= f1(r1,K1,UN,VN);
       NL2= f2(r2,K2,UN,VN);

       V1=UN + dt*NL1;
       V2=VN + dt*NL2;
    // iteration

   UN1=linsolve(M1,-V1);

   VN1=linsolve(M2,-V2);


   time(iter)=t;
   proie(iter)=UN1(N/2);
   pred(iter)=VN1(N/2);   

  t= t+dt;
  iter=iter+1;
  
  UN=UN1;
  VN=VN1;
   

 drawlater();
 clf();
  plot2d(x(2:N),UN,1,rect=[a,0,b,4]);
  plot2d(x(2:N),VN,2,rect=[a,0,b,4]);
  
 
 drawnow();   
end;  // fin boucle en temps




// trace de la solution
xset("window",1)
clf();

plot2d(x(2:N),UN,1);
plot2d(x(2:N),VN,2);
title('U (noir) et V (bleu) au temps final ');
xlabel('x');


xset("window",2)
clf();
plot2d(time,proie,1);
plot2d(time,pred,2);
xlabel('temps');
title('U et V en 0 ');

disp((1-K2)/(1-K1*K2),'U equilibre Theo');
disp(pred($),'U equilibre numerique');
disp((1-K1)/(1-K1*K2),'V equilibre Theo');
disp(proie($),'V equilibre numerique');



