%-----------------------------------------------------------------
% Exercise #3 'Discretization Methods'
% Calculation of a simply supported (Euler-Bernoulli) beam with 
% a linear (0 at x=0 and q at x=l) distributed load
% (analytical solution and weighted residuals)
% (c) 11/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/(4.0*pi^3.0);              % constant for subdomain method
c2 = 1.0/(2.0*pi^4.0);              % constant for collocation method
c3 = 1.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) = (3.0*x(i)^5.0-10.0*x(i)^3.0+7.0*x(i))*a/360.0;   % analytical solution
  %
  w1 = sin(pi*x(i));                              % shape function 1 value
  w2 = sin(2.0*pi*x(i));                          % shape function 2 value
  w_1(i) = c1*a*(w1-w2/32.0);                     % solution subdomain method
  w_2(i) = c2*a*(w1-(3.0-sqrt(2.0))/32.0*w2);     % solution collocation method
  w_3(i) = c3*a*(2.0*w1-w2/16.0);                 % 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')
xlabel('x')
ylabel('deflection of the beam')
title(['black: analytical   red: subdomain   blue: collocation   green: Galerkin'])