% 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 );