% Solution using the Galerkin methods of the boundary
% value problem:
%   - u'' + u = f   in (0,1)
%     u(0) = 0; u'(1) = b;

function fem_01
  
% Size of the problem;
N1 = 10;

% Viewing the basis functions
Nx = 100;
dx = 1.0 / (Nx - 1);
X = 0.0 : dx : 1.0;

figure(1); 
plot( X, phi( N1, 1, X ), X, phi( N1, 2, X ), X, phi( N1, 3, X ), ...
      X, phi( N1, 10, X ) );
legend( 'order 1', 'order 2', 'order 3', 'order 10' );


% Solving the system;
dx1 = 1.0 / N1;
X1 = 0.0 : dx1 : 1.0;

A = Assemble( N1 );
rhs = RHS( N1, 10 * N1 );

y = A \ rhs';

% plotting the solution
figure(2);
plot(X, f( X ), '-o', X, sol(X, y, N1), '-x', X, sol0( X ), '-');
legend('RHS', 'Approx. Sol', 'Exact Sol');

% visualizing the system matrix;
%figure(3); surf( log(abs(A)) );
figure(3); surf( (abs(A)) );

% Analyzing the condition number;
Nmin = 1e1;
Nmax = 1e3;

n = Nmin;
i = 1;
while ( n < Nmax )  
  C(i,1) = n;
  C(i,2) = cond( Assemble( n ) );
  C(i,3) = cond( Assemble0( n ) );
  n = n * 2;
  i = i + 1;
end;

figure(4); loglog(C(:,1), C(:,2),'-x', C(:,1), C(:,3),'-o' );
legend('FEM','x^k*(1-x)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The RHS function
function [fun] = f( x )  
  
  fun = x;
  
% The RHS vector in the basis chosen;  
function [G] = RHS( N, Ni )

dxi = 1.0 / Ni;
Xi = 0.0 : dxi : 1.0;
   
  for k=1:N
    s = 0.0;
    for m=2:Ni;      
     s = s + f( Xi(m) ) * phi( N, k, Xi(m) );
    end;
    G(k) = ( s + 0.5 * (f( Xi(1) ) * phi( N, k, Xi(1) ) + f( Xi(Ni+1) ...
						  ) * phi( N, k, Xi(Ni+1) ) ) ) * dxi;    
  end;

% Adding the boundary term;
  G(N) = G(N) + 1.0;
  
  

% The basis functions
function [fun] = phi( N, k, x )
  
  Nx = length( x );
  
  h = 1.0 / N;
  for i=1:Nx
    if ( (x(i) > (k-1) * h) & (x(i) < k * h) )
      fun(i) = ( x(i) - (k-1) * h ) / h;
    elseif ( (x(i) > k * h) & (x(i) < (k+1) * h) )
      fun(i) = ( (k+1) * h - x(i) ) / h;
    else
      fun(i) = 0.0;
    end;
  end;
  

% Assembling the matrix;
function [A] = Assemble( N )

  A(1:N,1:N) = 0.0;
  h = 1.0 / N;
  for k=1:N-1
% The main diagonal
    A(k,k)   = 2.0 * h / 3.0 + 2.0 / h;
% The upper diagonal
    A(k,k+1) = h / 6.0 - 1.0 / h;
% The lower diagonal
    A(k+1,k) = h / 6.0 - 1.0 / h;    
  end;
% THe last element;
  A(N,N) = h / 3.0 + 1.0 / h;
  

  
% Calculating the solution based on the expansion coeffients in the
% given basis
function [y] = sol( X, G, N )
  
  y = 0.0;
  for k=1:N
    y = y + G(k) * phi( N, k, X );
  end;
  
  
% Analytical solution for the problem f=-x;
function [y] = sol0( X )
  
  y = X;
  
  
% Assembling the matrix (using the pathological basis functions);
function [A] = Assemble0( N )
  
  for k=1:N
    for m=1:N
% $$$       A(j,k) = (j+1)*(k+1)/(j+k+1) - ((j+1)*(k+2)+(j+2)*(k+1))/(j+k+2) ...
% $$$ 	       + (j+2)*(k+2)/(j+k+3);
    
%      A(k,m) = 2*k*m / (k^3+3*k^2*m+3*k*m^2-k-m+m^3);
      A(k,m) = 2*(k*m^3+m^2+k^3*m-k+8*k*m+5*k*m^2-m+k^2+2*k^2*m^2+5*k^2*m)/(-6*k-6*m+5*k^3+5*m^3+m^5-10*k*m+k^5-5*k^2-5*m^2+20*k^3*m+30*k^2*m^2+15*k^2*m+15*k*m^2+20*k*m^3+10*k^3*m^2+5*k^4*m+10*k^2*m^3+5*k*m^4+5*k^4+5*m^4);
    end;
  end;
