//
// Equation de la chaleur 
// methode des éléments finis
/////////////////////////////////
/////////////////////////////////
function y=fcta(x)
    y=ones(x)
//    ind=find(x-floor(x)>1/2)   // fonction créneau
//    y(ind)=2*ones(ind)     
   y=2*ones(x)+cos(2*%pi*x)      // fonction cosinus
endfunction
//////////////////////////
function y=fctna(x)  
    y=ones(x)
//    ind=find(n*x-floor(n*x)>1/2) 
//    y(ind)=2*ones(ind) 
    y=2*ones(n*x)+cos(2*n*%pi*x) 
endfunction
/////////////////////////////////
function y=fct1sura(x)
    y=ones(x)./fcta(x)     
endfunction
/////////////////////////////////
function y=fctalim(x)
   y=alim
endfunction
/////////////////////////////////
function y=fctamean(x)
   y=amean
endfunction
/////////////////////////////////
function y=fctf(x)
   y=1    
endfunction  
/////////////////////////////////
/////////////////////////////////
// paramètres
n=12222;             // variations coefficient a
Nt=400003;        // discrétisation méthode rectangles   
N=129;h=1/N;        // discrétisation EF
//
alim=1/intg(0,1,fct1sura);   
amean=intg(0,1,fcta);
///////////////////////////////// 
//////////////////////////////////
// matrices A
//
A=zeros(N-1,N-1);
//
for i=1:N-2
////////////////////////////////////
    s1=0;xs1=(i-1)/N:1/Nt:(i+1)/N;l1=length(xs1)
    for k=1:l1;
       s1=s1+fctna(xs1(k))
     end   
     i1=s1*2/(N*l1)
//
    s2=0;xs2=(i)/N:1/Nt:(i+1)/N;l2=length(xs2)
    for k=1:l2;
       s2=s2+fctna(xs2(k))
     end   
     i2=s2/(N*l2)
//
    s3=0;xs3=(i-1)/N:1/Nt:(i)/N;l3=length(xs3)
    for k=1:l3;
       s3=s3+fctna(xs3(k))
     end   
     i3=s3/(N*l3)
///
    A(i,i)=N^2*i1;  // calcul des intégrale (avec rectangles)
    A(i,i+1)=-N^2*i2;  
    A(i+1,i)=-N^2*i3;    
////////////////////   
//    A(i,i)=N^2*intg((i-1)/N,(i+1)/N,fctna);// calcul des intégrale (avec intg)
//    A(i,i+1)=-N^2*intg((i)/N,(i+1)/N,fctna);  
//    A(i+1,i)=-N^2*intg((i-1)/N,(i)/N,fctna);  
end
s4=0;xs4=(N-2)/N:1/Nt:1;l4=length(xs4)
for k=1:l4
    s4=s4+fctna(xs4(k))
end   
i4=s4*2/(N*l4)  
A(N-1,N-1)=N^2*i4;
//A(N-1,N-1)=N^2*intg((N-2)/N,1,fctna);
//////////////////////////////////////////////////////
Alim=zeros(N-1,N-1);
//
for i=1:N-2
    Alim(i,i)=N^2*intg((i-1)/N,(i+1)/N,fctalim);
    Alim(i,i+1)=-N^2*intg((i)/N,(i+1)/N,fctalim);  
    Alim(i+1,i)=-N^2*intg((i-1)/N,(i)/N,fctalim);  
end
Alim(N-1,N-1)=N^2*intg((N-2)/N,1,fctalim);
/////////////////////////////////////////////////////
Amean=zeros(N-1,N-1);
//
for i=1:N-2
    Amean(i,i)=N^2*intg((i-1)/N,(i+1)/N,fctamean);
    Amean(i,i+1)=-N^2*intg((i)/N,(i+1)/N,fctamean);  
    Amean(i+1,i)=-N^2*intg((i-1)/N,(i)/N,fctamean);  
end
Amean(N-1,N-1)=N^2*intg((N-2)/N,1,fctamean);
///////////////////////////////////////////////////////
// second membre
//
F=zeros(N-1,1);
for i=1:N-1
  F(i)= 1/N; 
end
//
////////////////////////////////////:
u=inv(A)*F;
ulim=inv(Alim)*F;
umean=inv(Amean)*F;
//////
x=0:h:1;
x2=0:1/5656:1;
figure(0)
clf()
plot(x,[0;u;0],'b*-')              
plot(x,[0;ulim;0],'g*-') 
plot(x,[0;umean;0],'r*-')
plot(x2,1/(2*amean)*x2.*(ones(x2)-x2),'black');
/////////////////////////////////////
figure(1)
clf()

plot2d(x2,fctna(x2),rect=[0,0,1,3]);  // variation de a
//
