//Probleme Dirichlet/Neuman homogene 2D et zone penalisee
// Dirichlet sur x=0 et y=0, Neuman ailleurs
// correction TP 3
//
// taille du domaine [0,L]x[0,H]
L=2;
H=3;
// nombre de points de discretisation et pas de maillage
N=130;
M=150;
// on definit le bord x=L entre l'avant dernier et le dernier point de grille.
//Le but est de definir une DF centree d'ordre 2 pour la CL Neuman avec le dernier et l'avant dernier point
// soit l*(N-0.5)=L
l=L/(N+0.5);
h=H/(M+0.5);
// vecteur discretisation selon x
x=[l:l:L-0.5*l];
// vecteur discretisation selon y
y=[h:h:H-0.5*h];
//
xx=x'*ones(1,M);// repetition de x sur un tableau de taille (N,M) ie xx(i,j)=x(i)
yy=ones(N,1)*y;// ici yy(i,j)=y(j)
//
//construction de la matrice A par ses diagonales 
//
D0=(2/l^2+2/h^2)*sparse(ones(N*M,1));// diagonale principale
D0(N:N:N*M)=D0(N:N:N*M)-1/l^2// modif pour Neuman bord droit
D0(N*M-N+1:N*M)=D0(N*M-N+1:N*M)-1/h^2// modif pour Neuman bord haut
D1=-1/l^2*sparse(ones(N*M-1,1));// surdiagonale
D1(N:N:N*M-1,1)=0//correction de la surdiagonale CL
DN=-1/h^2*sparse(ones(N*M-N,1));
A=diag(D0,0)+diag(D1,1)+diag(D1,-1)+diag(DN,N)+diag(DN,-N);// A est pentadiagonale
//
// on rajoute a cette matrice une matrice diagonale pour penaliser a 0 la solution sur un disque au centre du domaine
// definition de la fonction caracteristique du disque
ray=max(L,H)/7;//rayon du disque
DISC=(xx-L/2).^2+(yy-H/2).^2-ray^2// negatif dans le disque positif sinon
DISC=-DISC./(2*abs(DISC)+1.d-12)+0.5;// tableau d'evaluation de la fonction caracteristique du disque
DISCV=matrix(DISC,N*M,1);
// ajout de penalisation
A=A+1.d10*diag(sparse(DISCV));
//construction d'un second membre a l'aide de la fonction f(x,y)=sin(x)*cos(y) (exemple)
ff=sin(xx).*cos(yy);// ff a la taille de xx et yy et f(i,j)=sin(x(i))*cos(y(j))
f=matrix(ff,N*M,1);// equivalent de reshape en matlab, 
                 //permet de construire f comme un vecteur colonne 
                 //avec la numerotation usuelle
//
//resolution du systeme lineaire creux
//
u=A\f;
uu=matrix(u,N,M)// permet de reconstruire uu comme un tableau a 2 indices, 
                //uu(i,j)=u(k), avec k la numerotation de TD
// representation graphique de la solution
uubord=zeros(N+1,M+1);
uubord(2:N+1,2:M+1)=uu;
surf(uubord)
