function [ ErrL2, ErrH1 ] = solveFEM_assign2( N )

dx = 1/N; % size of an element
x = 0:dx:1; % mesh

% form the Laplacian left hand side
[ L ] = solveFEM_assign2_laplacian( N );
% estimate the condition number of the matrix
c = condest( L );
disp([' The condition number of the matrix is: ',num2str(c)]);

% form the right hand side
[ R ] = solveFEM_assign2_R( N );

% solve the problem
u = L \ R;

% compute the exact solution on a very dense grid
xdense = 0:0.1*dx:1;
uopt = cos( 0.5*pi*xdense' );

%plot the solution, note that we have to add the boundary conditon 1 as the
%first entry of the solution
uu = [ 1; u; ];
de = N* uu(1:N);
de = de - N* uu(2:N+1);
xx = dx/2:dx:1-dx/2;
duopt = -0.5*pi*sin( 0.5*pi*xx' );

plot( x, [ 1; u], xdense, uopt, xx, -de, xx, duopt );

ErrL2 = solveFEM_assign2_L2( N, [ 1; u; ] );
ErrH1 = solveFEM_assign2_H1( N, [ 1; u; ] );


function [ L ] = solveFEM_assign2_laplacian( N )
% this is dealing with MATLAB and sparse matrixes
% lookup compressed sparse format
pntr = zeros( 3*N - 2, 1 );
indx = zeros( 3*N - 2, 1 );
vals = zeros( 3*N - 2, 1 );
% 
pntr( 1:3:end ) = 1:N;
pntr( 2:3:end ) = 1:N-1;
pntr( 3:3:end ) = 2:N;
% 
indx( 1:3:end ) = 1:N;
indx( 2:3:end ) = 2:N;
indx( 3:3:end ) = 1:N-1;
% 
vals( 1:3:end ) = 2 * ones( N, 1 );
vals( 2:3:end ) = -ones( N-1, 1 );
vals( 3:3:end ) = -ones( N-1, 1 );
% 
L = -N * sparse( pntr, indx, vals, N, N, 3*N-2 );

function [ R ] = solveFEM_assign2_R( N )
dx = 1/N;
mk = 0.5*dx : dx : 1; % get the mid-points of the elements
%R = dx*0.5* (-pi^2/4) * cos( 0.5*pi*x );
R = -dx*0.125*pi^2*cos( 0.5*pi*mk' );
R(1:N-1) = R(1:N-1) + R(2:N);
% Add the boundary conditions
R(1) = R(1) - N;
R(N) = R(N) + 0.5*pi;

function [ ErrL2 ] = solveFEM_assign2_L2( N, u )
% compute the L2 error for the computer solution u
% use mid point quadrature rule on each element
dx = 1 /N;
mk = 0.5*dx : dx : 1;
% get the values of the exact solution at each point
uopt = cos( 0.5*pi*mk' );
% get the difference
e = uopt - 0.5*u(1:N);
e = e - 0.5*u(2:N+1);
% square and get the error
ErrL2 = sqrt( sum( sum( dx*e.^2 ) ) );

function [ ErrH1 ] = solveFEM_assign2_H1( N, u )
% compute the H1 error for the computer solution u
% use mid point quadrature rule on each element
%dx = 1 /N;
%mk = 0.5*dx : dx : 1;
% get the values of the exact solution at each point
%uopt = cos( 0.5*pi*mk' );
% get the difference
%e = uopt - 0.5*u(1:N);
%e = e - 0.5*u(2:N+1);
% get the error in the derivative
%duopt = -0.5*pi*sin( 0.5*pi*mk' );
%de = duopt + N* u(1:N);
%de = de - N* u(2:N+1);
% square and get the error
%ErrH1 = sqrt( sum( sum( dx*(e.^2) + dx*(de.^2) ) ) );

%%% Second method for computing the H1 error %%%
% the error in the derivative is computed analytically, without the use of a
% quadrature
ErrL22 = solveFEM_assign2_L2( N, u )^2;

x = 0:1/N:1;
I = 0;
for i = 1:N
    c = N*( u(i+1) - u(i) );
    I = I + 0.25*pi^2*( x(i+1)/2 - x(i)/2 + (sin(pi*x(i))/2 - sin(pi*x(i+1))/2)/pi ) - 2*c*( cos( 0.5*pi*x(i+1) ) - cos( 0.5*pi*x(i) ) ) + c^2 * ( x(i+1) - x(i) );
end;
ErrH1 = sqrt( ErrL22 + I );
