/////////////////////
// TP1, exercice 2 //
/////////////////////


//Choix du problème à étudier 
num_probleme=3;//À choisir entre 1 et 4, pour les quatres seconds membres différents

//On ouvre une nouvelle fenêtre
scf(0);clf();

//Le second membre
function y=second_membre_fct(x)
  if num_probleme==1
	y=x.^2;
  elseif num_probleme==2
	y=-6*abs(x-1/2);
  elseif num_probleme==3
	y=-2*(x>1/2)+3*(x>3/4);
  else
	y=50*(x>0.49).*(x<0.51);//ce n'est qu'un approximation, dont on peut
                            //en fait se passer
  end;
endfunction;

//La solution exacte
function y=sol_ex_fct(x)
  if num_probleme==1
	y=(x.^4-x)/12;
  elseif num_probleme==2 
	y=1/8-abs(x-1/2)^3;
  elseif num_probleme==3
	y=(5*x).*(x<1/2)+(-32*x.^2+37*x-8).*(x>=1/2).*(x<3/4)+(16*x.^2-35*x+19).*(x>=3/4);
	y=y/32;
  else
	y=(abs(2*x-1)-1)/4;
  end;
endfunction;


for exposant=1:10
  //Le nombre de pas de discrétisation
  N=2^exposant;
  
  //Le pas d'espace
  dx=1/(N+1);
  
  //La grille d'espace
  grille=[0:N+1]'*dx;
  
  //Définition de la matrice du problème
  A=-2*diag(ones(1,N))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
  A=A/(dx^2);
  
  //Le second membre
  f=second_membre_fct([1:N]'*dx);//Ce n'est pas vraiment le second
                                 //membre donné par la méthode des
                                 //éléments finis
  
  //La solution exacte
  x=linspace(0,1,1000);
  sol_ex=sol_ex_fct(x);
  
  
  ///////////////////////////////////////////////
  
  
  //On resoud le schéma numérique
  sol_app=A\f;
  
  //On rajoute les conditions au bord
  sol_app=[0 ; sol_app ; 0];
  
  
  //////////////////////////////////////////////
  
  
  //Tracé de la solution approchée
  plot2d(grille,sol_app,5);//5 = en rouge
  
end;
//Tracé de la solution exacte
plot2d(x,sol_ex,13);//13 = en vert
