 X = [];
 Y = [];




//  on charge les programmes dont on a besoin
//
exec('P1.sci'); exec('gauss.sci');
//  nbre de points pour l'integration numerique
Ng = 3;
[XG,PG] = gauss(Ng); 
//
//
//Nelt = 20; // nbre d'elements
//---------------------------------------------
for Nelt=5:5:100

Nno = 2; // nbre de noeud par element

Npt = Nelt * (Nno-1) +1; // nbre de points

Dx = 1./(Npt-1);

COOR = matrix((0:Dx:1),Npt,1);  //  Coordonnees des points

CONNEC = zeros(Nelt,Nno);
for i=1:Nelt
   for j=1:Nno
      CONNEC(i,j) = j + (i-1);  // matrice de connectivite
   end;                        
end;
//
//    Fin initialisation des tableaux

A_glob = zeros(Npt,Npt); F_glob = zeros(Npt,1);

for iel = 1:Nelt    
   // calcul du syteme elementaire sur chaque element
   A_el = zeros(Nno,Nno);
   F_el = zeros(Nno,1);

   x1 = COOR(CONNEC(iel,1)); x2 = COOR(CONNEC(iel,2));
   jac = (x2 - x1)/2.;   // le jacobien est constant sur l'element

   for ig=1:Ng    // boucle pour l'integration numerique
      [FP,DxFP] = P1(XG(ig));
      
      DxFP = DxFP/jac;

       for i=1:Nno
          F_el(i) = F_el(i) + PG(ig)*FP(i)*jac;
          for j=1:Nno
             A_el(i,j) = A_el(i,j) + PG(ig)*(DxFP(j)*DxFP(i) + FP(i)*FP(j))*jac;
          end;
       end;
    end;
   // fin calcul systeme elementaire pour l'element iel

// on doit assembler

   for i=1:Nno
      i_glob = CONNEC(iel,i);
      F_glob(i_glob) = F_glob(i_glob) + F_el(i);
      for j=1:Nno
         j_glob = CONNEC(iel,j);
         A_glob(i_glob,j_glob) = A_glob(i_glob, j_glob) + A_el(i,j);
      end;
    end;
end;  // fin boucle sur les elements

//
// on impose les deux conditions aux limites de Diriclet
//   U(0) = U(1) = 0

A = [zeros(1,Npt); A_glob(2:(Npt-1),:); zeros(1,Npt)];
A(1,1) = 1; A(Npt,Npt) = 1;
F = [0; F_glob(2:(Npt-1)); 0];

u_cal = A\F;

// comparaison avec la solution exacte
//
x = COOR;
u_exact = 1- (exp(x) + exp(1-x))/(1+exp(1));
erreur = norm(u_cal-u_exact)/norm(u_exact);

X = [X Npt];
Y = [Y erreur];
end;

xbasc();
//plot2d('ll',X,Y,rect=[1 10^(-2) 100 10^(-5)]);
//plot2d('ll',X,Y);
plot2d(log10(X),log10(Y),rect=[1 -2 2 -6]);
xtitle('Elements P1','Log(Npt)','Log(Erreur)');

