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

function galerkin_01
  

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

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



% Solving the system;
N1 = 50;
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 ), '-');

% visualizing the system matrix;
figure(3); surf( log(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 ) );
  n = n * 2;
  i = i + 1;
end;

figure(4); loglog(C(:,1), C(:,2),'-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( k, Xi(m) );
    end;
    G(k) = ( s + 0.5 * (f( Xi(1) ) * phi( k, Xi(1) ) + f( Xi(Ni+1) ) * phi( k, Xi(Ni+1) ) ) ) * dxi;    
  end;
  
  

% The basis function
function [fun] = phi( k, x )
  
  fun = x.^k .* ( 1.0 - x );
  

% Assembling the matrix;
function [A] = Assemble( 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);
    end;
  end;

  
% 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( k, X );
  end;
  
  
% Analytical solution for the problem f=-x;
function [y] = sol0( X )
  
  y = 1.0 / 6.0 * X .* ( X.^2 - 1.0 );
