function [ conerr, XX, YY, merrs, merrv ] = fdmerr ( N );
% FDMErr
%
% This m-file is part of the material of the lecture
%   Theory and Numerics of Partial Differential Equations
% at the Math Department of Halle University, Germany.
%
% Author :      Prof. Dr. M. Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Jun 3, 2005
%
% Consistency and convergence of 5-point-stencil FDM approximation (Error)
%
% Parameters:
%   N      (input)  : number of grid points in FD approximation
%   conerr (output) : vector of errors
%                       conerr(1) = \|A_h*U-q_h\|_\infty
%                       conerr(2) = \|A_h*U-q_h\|_h
%                       conerr(3) = \|u_h-U\|_\infty
%                       conerr(4) = \|u_h-U\|_h
%   XX     (output) : x-grid for array output
%   YY     (output) : y-grid for array output
%   merrs  (output) : error A_h*U-q_h in array format
%   merrv  (output) : error u_h-U in array format
%
% Examples:
%   Nvec = zeros ( 6, 1 );    Nvec(1) = 4;
%   conerr = zeros ( length(Nvec), 4 );
%   for iN=1:length(Nvec),  
%       [ conerr(iN,:), XX, YY, merrs, merrv ] = fdmerr ( Nvec(iN) );  
%       if iN<length(Nvec),  Nvec(iN+1) = 2*Nvec(iN);  end;
%   end;
%   loglog ( 1./Nvec, conerr );

% -> finite difference approximation of Poisson problem in 2D
h = 1/N;
x = linspace ( 0, 1, N+1 );
[ XX, YY ] = meshgrid ( x );

% -> exact solution
e  = ones ( (N-1)^2, 1 );
ee = ones ( (N-1)^2, 1 );
ee(N-1:N-1:(N-1)^2) = 0;
L = spdiags ( [ e ee ] / h^2, [ -(N-1) -1 ], (N-1)^2, (N-1)^2 );
A = - ( L + spdiags ( -4*e / h^2, 0, (N-1)^2, (N-1)^2 ) + L' );

qmat = 5*pi^2 * sin(pi*XX) .* sin(2*pi*YY);
qvec = zeros ( (N-1)^2, 1 );
for i=1:N-1,
    qvec((i-1)*(N-1)+(1:N-1)) = qmat(1+(1:N-1),1+i);
end;

xnum = A \ qvec;
Unum = zeros ( N+1, N+1 );
Unum(1+(1:N-1),1+(1:N-1)) = reshape ( xnum, N-1, N-1 );

Usol  = sin(pi*XX) .* sin(2*pi*YY);
merrv = Unum - Usol;

xsol = zeros ( (N-1)^2, 1 );
for i=1:N-1,
    xsol((i-1)*(N-1)+(1:N-1)) = Usol(1+(1:N-1),1+i);
end;

errs  = A*xsol - qvec;
merrs = zeros ( N+1, N+1 );
merrs(1+(1:N-1),1+(1:N-1)) = reshape ( errs, N-1, N-1 );

conerr = [   norm(errs,inf) 
           h*norm(errs) 
             norm(xsol-xnum,inf) 
           h*norm(xsol-xnum) ]';
