//
//  equation de Laplace avec conditions de Dirichlet sur [0,1]^2, second membre constat
//
N=evstr(x_dialog('nombre de points intérieurs','10'));
fval=evstr(x_dialog('valeur de f(x) constante','10'));
//
M=N;
h=1/(N+1);
//////////////////////////////////////////
//  résolution éléments finis Q1
//  
////////////////////////////////////////////
//
//---------------------------------
//  matrice A (cas M=N)
L1=toeplitz([-1/3 -1/3 zeros(1,N-2)]);
L2=toeplitz([8/3 -1/3 zeros(1,N-2)]);	//tridiagonale par blocs
A=zeros(N*N,N*N);
for i=1:(N-1)
  i1=(i-1)*N;
  i2=i*N;
  A((i1+1):(i1+N),(i1+1):(i1+N)) = L2;
  A((i2+1):(i2+N),(i1+1):(i1+N)) = L1;
  A((i1+1):(i1+N),(i2+1):(i2+N)) = L1;
end;
i1=(N-1)*N;
A((i1+1):(i1+N),(i1+1):(i1+N)) =L2;
//
//-------------------------------
// second membre
// 
b=zeros(N*N,1);
for j=1:N
 for i=1:N
      b((j-1)*N+i)=h^2*fval; 
    end
end
//----------------------------------
// resolution de AU=b
//
U=inv(A)*b;
//----------------------------------
//
// reecriture sous forme matrice de U (rajout des bords nuls)
//
umat=zeros(N+2,N+2)
for j=2:N+1
    for i=2:N+1
        umat(i,j)=U((j-2)*N+i-1)  // colonne par colonne (adaptée à plot3d)
    end
end
////////////////////////////////////////////////////
/// Différences finies
////////////////////////////////////////////////////
//---------------------------------
//  matrice A (cas M=N)
L1=-eye(N,N);
L2=toeplitz([4 -1 zeros(1,N-2)]);	//tridiagonale par blocs
A=zeros(N*N,N*N);
for i=1:(N-1)
  i1=(i-1)*N;
  i2=i*N;
  A((i1+1):(i1+N),(i1+1):(i1+N)) = L2;
  A((i2+1):(i2+N),(i1+1):(i1+N)) = L1;
  A((i1+1):(i1+N),(i2+1):(i2+N)) = L1;
end;
i1=(N-1)*N;
A((i1+1):(i1+N),(i1+1):(i1+N)) =L2;
A=1/(h^2)*A;
//
//-------------------------------
// second membre
//
b=zeros(N*N,1);
for j=1:N
 for i=1:N
      b((j-1)*N+i)=fval; 
    end
end
//----------------------------------
// resolution de AU=b
//
U=inv(A)*b;
//----------------------------------
//
// reecriture sous forme matrice de U (rajout des bords nuls)
//
umatDF=zeros(N+2,N+2)
for j=2:N+1
    for i=2:N+1
        umatDF(i,j)=U((j-2)*N+i-1)  // colonne par colonne (adaptée à plot3d)
    end
end
//////////////////////////////////
//----------------------------------
figure(1)
clf()
//  tracé
subplot(131)
plot3d(0:h:1,0:h:1,umat) 
xtitle('Poisson par éléments finis Q1')      
subplot(132)
plot3d(0:h:1,0:h:1,umatDF) 
xtitle('Poisson par différences finies') 
subplot(133)
plot3d(0:h:1,0:h:1,abs(umatDF-umat))
xtitle('écart entre les deux méthodes') 
////////////////////////////////////
disp('ecart maximal entre les deux méthodes:');disp(max(abs(umatDF-umat)));
