
///// Resolution de l'eq de fisher Ut= f(U)) + Uxx sur [a,b]
////  par differences finies. Schema Euler implicite pour diffusion
///  explicite pour terme NL
///  Guillaume Chiavassa, mars 2017

//  avec f(U)= 1-u si r<u<1  et f(u)=-u si 0<u<r  avec r entre 0 et 1

clear;


// matrice tri diagonale du laplacien en dim 1, differences finies centrees
//    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;


//terme nonlineaire, r,K scalaires, U vecteur colonne


function [sol]= bistable(r,U)
 for i=1:length(U),
    if U(i) < r then 
        sol(i)=-U(i);
    else
        sol(i)= 1-U(i);
    end,
    
end;
        
endfunction;


// condition initiale une gaussienne centree en x0

 function U=cond_init(x0,x)

  U=exp(-(x-x0).^2)';
  
endfunction;

function U=solexacte(c,r,x,t)

  lambdaR=0.5*(c+sqrt(c^2+4));
  lambdaL=0.5*(-c+sqrt(c^2+4));
  z=x-c*t;
  
  for i=1:length(z)
      if z(i)<0 then
          U(i)=1-(1-r)*exp(lambdaL*z(i));
      else
          U(i)=r*exp(-lambdaR*z(i));
      end,
      
  end
  
endfunction;



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


 r=0.3;

 a=-35; b=35; 

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

//temps final
T=17;

dt=0.01;



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

// les variables

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

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

 coef=dt/dx^2;

 B=laplacien(coef,N-1);

//calcul de l'inverse pour l'implicite
  
  M=eye(N-1,N-1)-B;

// initialisation de VN

  VN=U0(2:N);

// BOUCLE EN TEMPS
t=0;
iter=1;

while t<T,
  
    // terme NonLin
      
      NL= bistable(r,VN);
    
      V=VN + dt*NL;

   // iteration
      VN1=linsolve(M,-V);

   time(iter)=t;
   //enregistrement pour calcul de la vitesse
  [rien(1,iter),xvit(1,iter)]=min(diff(VN1)  );   

  t= t+dt;
  iter=iter+1;

  VN=VN1;
  drawlater();
  clf();
  plot2d(x(2:N),VN);
  drawnow();
end;


// trace de la solution
xset("window",1)
clf();
plot2d(x(2:N),VN,-1);

ctheo=(1-r/(1-r))*sqrt((1-r)/r);


Utheo=solexacte(ctheo,r,x,t);

plot2d(x,Utheo,2);
title('SOLUTION (noir); onde progressive theorique (bleu)');
xlabel('x');

// trace la position du front en fct de t

xset("window",2)
clf();
plot2d(time,xvit*dx+a);

// estimation de la vitesse sur les 20 dernieres iterations
[vit,toto]=reglin(time($-40:$)',xvit($-40:$)*dx+a)

plot2d(time, (a+b)/2 + time*ctheo,2)

title('position du Front en fct du temps. calcul (noir), theorique (bleu) ');
xlabel('temps');


disp(ctheo,'vitesse theorique: ')
disp(vit,'vitesse mesurée: ')





