//
//  programme 6
//
write(%io(2),'programme 6: resolution de l''equation de Laplace sur [0,1]^2')
write(%io(2),'//')
//
xbasc()
//
deff('mat=laplace(d)',[   					 //matrice du laplacien
     'L1=eye(d,d);'						//tridiagonale par blocs
     'L2=toeplitz([-4 1 zeros(1,d-2)]);'
     'mat=zeros(d*d,d*d);'
     'for i=1:(d-1)'
	'i1=(i-1)*d;'
	'i2=i*d;'
	'mat((i1+1):(i1+d),(i1+1):(i1+d)) = L2;'
	'mat((i2+1):(i2+d),(i1+1):(i1+d)) = L1;'
	'mat((i1+1):(i1+d),(i2+1):(i2+d)) = L1;'
    'end;'
    'i1=(d-1)*d;'
   ' mat((i1+1):(i1+d),(i1+1):(i1+d)) =L2'])
//
d=evstr(x_dialog('pas de discretisation (>=5)?','12'));
//
npt= d+2; h=   1.0/(npt-1);Lapl=(1/h^2)*laplace(d);
//
// second membre (fonction creneau)
//
zone=ceil(d/2);
Smm=zeros(d,d);
Smm(zone-2:zone+2,zone-1:zone+1)=-4*(1/(h^2))*ones(5,3);
Smcol=Smm(:);
//
Solcol=Lapl\Smcol;   //resolution du systeme lineaire
//
Solm=zeros(d,d);
Solm(:)=Solcol;
Solm=[zeros(1,d);Solm;zeros(1,d)];
Solm=[zeros(d+2,1) Solm zeros(d+2,1)];
//
// trace du resultat
//
t=1:d+2;
xsetech([0,0,0.5,1])
plot3d(t,t,Solm);
xtitle('solution de l''equation de Laplace')
xsetech([0.5,0.,0.5,1])
contour(t,t,Solm,5);
xtitle('lignes de niveau')



