%%% Solution probabiliste du problème de Dirichlet sur un carré
clear
cote=10;
nbreal=100;
%
Bord=zeros(cote,cote);
Bord(1,1:cote)=1;
Bord(cote,1:cote)=1;
Bord(1:cote,1)=1;
Bord(1:cote,cote)=1;
%
F=zeros(cote,cote);
V=cumsum(ones(1,cote));
F(1,1:cote)=V;
F(cote,1:cote)=fliplr(V);
F(1:cote,1)=V';
F(1:cote,cote)=flipud(V');
%
MaxL(1:cote,1:cote)=zeros(cote,cote);
%
Sol=F;
%
for x=2:cote-1,
  for y=2:cote-1,
	N=zeros(cote,cote);
	for m=1:nbreal,
	  sortie=0;
	  n=1;
	  Pos(n,1:2)=[x,y];
	  while (sortie==0)
		axe=1+(rand<1/2);
		increment=2*(rand<1/2)-1;
		Pos(n+1,1:2)=Pos(n,1:2);
		Pos(n+1,axe)=Pos(n,axe)+increment;
		MaxL(x,y)=max(MaxL(x,y),n);
		if Bord(Pos(n+1,1),Pos(n+1,2))==1,
		  sortie=1;
		  N(Pos(n+1,1),Pos(n+1,2))=N(Pos(n+1,1),Pos(n+1,2))+1;
		else
		  n=n+1;
		end
	  end
	end
	Sol(x,y)=sum(sum(F.*N))/nbreal;
  end
end
%
figure(1); clf; 
mesh(Sol);
figure(2); clf; 
imagesc(Sol,ceil(300/cote));
figure(3); clf; 
imagesc(MaxL,ceil(300/cote));
%
Lap=zeros(cote,cote);
for x=2:cote-1,
  for y=2:cote-1,
	Lap(x,y)=Sol(x-1,y)+Sol(x+1,y)+Sol(x,y-1)+Sol(x,y+1);
	Lap(x,y)=Lap(x,y)/4-Sol(x,y);
  end
end
%
figure(4); clf;
mesh(Lap(2:cote-1,2:cote-1));
figure(5); clf; 
imagesc(Lap,ceil(300/cote));
%