
#tests of FDM solver for -Laplacian u = f in [a,b]^2
1==1;
function [u,x]=FDsolver2d(f,g=@(x,y) 0.*f(x,y),N=40,a=0,b=1,c=0)
#FD solver for -Laplacian u+cu=f on [a,b]^2 u(s)=g on boundary
#Input: f - function handle to f(x,y) should work for matrix 
#a,b- ends of interval
#N - no of mesh points including boundary -1 x_k k=0,..,N x_0 x_N on bnd  
#c - constant
#g = function handle for bnd condition
#output:
#x- mesh points in x (or y) direction 
#u - FD solution on these mesh (interior mesh points - 5 point stencil)
#Examples:
#f=@(x,y) 2*pi*pi*sin(pi*x).*sin(pi*y); u=FDsolver2d(f);mesh(u);
#f=@(x,y) 2*pi*pi*sin(pi*x).*sin(pi*y);g=@(x,y) x+y;u=FDsolver2d(f,g);mesh(u)
#g=@(x,y) exp(x).*sin(pi*y); f=@(x,y) (-1.0+pi*pi)*g(x,y);[u,x]=FDsolver2d(f,g);mesh(u);
#[XX,YY]=meshgrid(x,x);mesh(g(XX,YY)-u)
if(c<0)
 printf("Warning: c=%g<0. The BVP -u''+cu=f may have no solution!\n",c);
endif
h=(b-a)/N;

x=a+h*(0:N)';
dd=-ones(N-1,3);dd(:,2)=-2*dd(:,2);
#1d laplacian  (N-1)x(N-1) (scaled by h^2)
A=spdiags(dd,[-1,0,1],N-1,N-1);
Id=spdiags(ones(N-1,1),0,N-1,N-1);
L=kron(A,Id)+kron(Id,A);

ff=zeros(N-1,N-1);
xi=x(2:N);
[XX,YY]=meshgrid(xi,xi);
ff=f(XX,YY);
[XX,YY]=meshgrid(x,x);
gg=g(XX,YY);

ff*=h*h;
#upper bnd y=b
ff(:,N-1)+=gg(2:N,N+1);
#lower bnd y=a
ff(:,1)+=gg(2:N,1);
#left bnd x=a
ff(1,:)+=gg(1,2:N);
#right bnd x=b
ff(N-1,:)+=gg(N+1,2:N);

F=reshape(ff,(N-1)*(N-1),1);
#G=reshape(gg,(N+1)*(N+1),1);


U=L\F;
u=[gg(2:N,1),reshape(U,N-1,N-1),gg(2:N,N+1)];
u=[gg(1,:);u;gg(N+1,:)];
#u=[al; u; be];
endfunction





#tests of convergence order
f=@(x,y) 2*pi*pi*sin(pi*x).*sin(pi*y); g=@(x,y) 0*f(x,y);
sol=@(x,y) sin(pi*x).*sin(pi*y);
Nt=5;
N=10;erm=erL2=0;
errm=errL2=zeros(Nt,1);
for k=1:Nt,
  [u,x]=FDsolver2d(f,g,N);
  [XX,YY]=meshgrid(x,x);
  S=sol(XX,YY);
  ermp=erm;erL2p=erL2;
  errm(k)=erm=norm(reshape(S-u,(N+1)*(N+1),1),'inf');  
  errL2(k)=erL2=norm(reshape(S-u,(N+1)*(N+1),1),2)/N; 
  N*=2;
endfor
errm
 errm(1:Nt-1)./errm(2:end)
 errL2
 errL2(1:Nt-1)./errL2(2:end)




