%-----------------------------------------------------------------
% Calculation of a simply supported (Euler-Bernoulli) beam with 
% a constant distributed load
% (analytical solution and weighted residuals)
% (c) 10/2004 by A. Schmidt, IAM, Uni Stuttgart
%-----------------------------------------------------------------
clear all

% Variables
l = 1.0;                            % length of the beam
EI = 1.0;                           % (bending) stiffness
q = 1.0;                            % distributed load

n = 501;                            % number of sample points
                                    % (for plot)
				    
% Parameters
a = q*l^4/EI;                       % abbreviation
c1 = 1.0/(2.0*pi^3.0);              % constant for subdomain method
c2 = 1.0/(pi^4.0);                  % constant for collocation method
c3 = 4.0/(pi^5.0);                  % constant for least squares method
c4 = 4.0/(pi^5.0);                  % constant for Galerkin's method

%-----------------------------------------------------------------
% calculation of the curves
for i=1:1:n
  x(i) = (i-1)/(n-1);                             % sample points (x/l)
  w_a(i) = (x(i)^4.0-2.0*x(i)^3.0+x(i))*a/24.0;   % analytical solution
  %
  w = sin(pi*x(i));                               % shape function value
  w_1(i) = c1*a*w;                                % solution subdomain method
  w_2(i) = c2*a*w;                                % solution collocation method
  w_3(i) = c3*a*w;                                % solution least squares method
  w_4(i) = c4*a*w;                                % solution Galerkin's method
end

%-----------------------------------------------------------------
% plot
figure(1)
plot(x,-w_a,'k',x,-w_1,'r',x,-w_2,'b',x,-w_3,'g',x,-w_4,'m-.')
xlabel('x')
ylabel('deflection of the beam')
title(['black: analytical   red: subdomain   blue: collocation   green: least squares   magenta: Galerkin'])