Code covered by the BSD License  

Highlights from
Galerkins method over "ne" elements for solving 2nd-order homogeneous, c.c BVP

image thumbnail

Galerkins method over "ne" elements for solving 2nd-order homogeneous, c.c BVP

by

 

Implement Galerkin method over "ne" individual elements for solving 2nd order BVPs

BVP_Galerkin2(a,b,c,t1,t2,x1,x2,ne)
%The purpose of this program is to implement Galerkin method over "ne"  
% individual element for solving the following general 2nd order, 
% homogeneous, Boundary Value problem (BVP) with constant coefficients, and
% then comparing the answer with the exact solution.
%==========================================================================
%               ax"(t)+bx'(t)+cx(t)=0 for t1<=t<=t2
%                   BC: x(t1)=x1 and x(t2)=x2
% >> BVP_Galerkin(a,b,c,t1,t2,x1,x2,ne)
% where "ne" is the number of elements
% The output of this program is 
% 1- The approximated x(t) vs. exact  x(t)
% 2- The approximated x'(t) vs. exact x'(t)
% 3- The approximated x"(t) vs. exact x"(t)
% ======Example============================================================
% Equation                  x"(t)+ 0.5x'(t)+ 10x(t)=0
% Boundary values           x(1)=2,     x(10)=0;
% Solution: We have:
%                           a=1;    b=2;   c=3;
%                           t1=1;       t2=10;          
%                           x1=2;       x2=0;
% Using ne=64 elements,
% >>BVP_Galerkin2(1,2,3,1,10,2,0,64)
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%                   Sep.26 to Dec26, 2010. Ramin Shamshiri  
%                       ramin.sh@ufl.edu
%            Doctoral student at the University of Florida
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
%=============================Program Begin================================
function BVP_Galerkin2(a,b,c,t1,t2,x1,x2,ne) % Declare function
 
%% Begin Approximate solution via Galerkin method over each individual element 
 
% Define element length
le = (t2-t1)/ne; 
% Define t matrix: t=t1:t2
t_GM=t1:le:t2;
t_GM=t_GM';
 
% Building matrices arrays
K1 = (a/le) * [1,-1;-1,1];      % K1 matrix corresponding to the x"(t)
K2 =  b*[-1/2 1/2;-1/2 1/2];    % K2 matrix corresponding to the x'(t)
K3 =  -c*le*[1/3 1/6;1/6 1/3];  % K3 matrix corresponding to the x(t)
Ke = K1+K2+K3;                  % Element stiffness Matrix
 
%*****************Begin Assembly Global stiffness matrix*******************
k = zeros(ne+1);
for i=1:ne+1
  for j=1:ne +1
      if (i==j)
          if(i==1)
              k(i,j)=Ke(1,1);
          elseif(i==ne+1)
              k(i,j)=Ke(2,2);
          else
              k(i,j)=Ke(1,1)+Ke(2,2);
          end
      elseif(i==j+1)
          k(i,j)=Ke(1,2);
      elseif(j==i+1)
          k(i,j)=Ke(2,1);
      else
          k(i,j)=0;
      end
  end
end
%********************End Assembly Global stiffness matrix******************
 
%The Global f Matrix
f = zeros(ne+1,1);
 
%BC apply x(t1) = x1
f(1,1) = x1;
%BC apply x(t2) = x2
f(ne+1,1) = x2;
 
% Display the Global stiffness matrix before striking row
K_Global=k;
 
%Striking first and last rows
for i=2:ne+1
  k(1,i) = 0;
  k(ne+1,i) = 0;
end
k(1,1) = 1;
k(ne+1,ne+1) = 1;
 
 
% Display the solvable stiffness matrix 
disp('The solvable stiffness matrix is')
K_strike=k;
 
%solving the result and finding the displacement matrix, {x}
x_GM=inv(k)*f;
 
% calculating x'(t)
dx_GM=zeros(ne+1,1);
for i=1:ne
    dx_GM(i,1)=(x_GM(i+1)-x_GM(i))/le;
end
 
% calculating x"(t)
d2x_GM=zeros(ne+1,1);
for i=2:ne
    d2x_GM(i,1)=(dx_GM(i+1)-dx_GM(i))/le;
end
 
% Plotting the Approximate solution x(t) in red (Galerkin method over each individual element)
figure1 = figure('Color',[1 1 1],'units','normalized','outerposition',[0 0 1 1]);
subplot(2,2,1)
plot(t_GM,x_GM, '--', 'Color','r', 'LineWidth',3);  hold on
xlabel('Time (sec)','FontSize',12);
ylabel('x(t) (m)','FontSize',12);
 
% Plotting the approximate dx
subplot(2,2,2)
plot(t_GM,dx_GM, '--', 'Color','r', 'LineWidth',3); hold on
xlabel('Time (sec)','FontSize',12);
ylabel('dx/dt (m/s)','FontSize',12);
 
% Plotting the approximate d2x
subplot(2,2,3)
plot(t_GM,d2x_GM, '--', 'Color','r', 'LineWidth',3); hold on
xlabel('Time (sec)','FontSize',12);
ylabel('d^2x/dt^2 m/s^2','FontSize',12);
 
assignin('base', 'x_GM', x_GM);       % Sending x_GM to workspace
assignin('base', 'dx_GM', dx_GM);     % Sending dx_GM to workspace
assignin('base', 'd2x_GM', d2x_GM);   % Sending d2x_GM to workspace
assignin('base', 't_GM', t_GM);       % Sending t_GM to workspace
clc
 
%% End of Approximate solution via Galerkin method over each individual element 
 
 
%% ==============Begin Calculating Exact Soution========================== 
dt=0.001; % increment
r=roots([a b c]);   % Roots of auxiliary equation
r1=r(1);  r2=r(2);
 
% Bulding t-axis values
t=t1:dt:t2; t=t'; 
s=size(t);
 
% Bulding x(t) axis values
x_exact=zeros(s(1),1);
if r1==r2   % For critically damped case
    
    C=inv([exp(r1*t1)  t1*exp(r2*t1);exp(r1*t2)  t2*exp(r2*t2)])*([x1;x2]);
    C1=C(1);
    C2=C(2);
    
    for i=1:s(1)
        x_exact(i)=C1*exp(r1*t(i))+C2*t(i)*exp(r2*t(i));
    end
    
else        % For Under-damped or Over-damped case
    
    C=inv([exp(r1*t1) exp(r2*t1);exp(r1*t2) exp(r2*t2)])*([x1;x2]);
    C1=C(1);
    C2=C(2);
        
    for i=1:s(1)
        x_exact(i)=C1*exp(r1*t(i))+C2*exp(r2*t(i));
    end
end
 
% calculating x'(t)
dx_exact=zeros(s(1),1);
for i=2:s(1)
    dx_exact(i,1)=(x_exact(i)-x_exact(i-1))/dt;
end
 
% calculating x"(t)
d2x_exact=zeros(s(1),1);
for i=3:s(1)
    d2x_exact(i,1)=(dx_exact(i)-dx_exact(i-1))/dt;
end
 
% Plotting the exact solution x(t) in blue
subplot(2,2,1)
plot(t,x_exact, 'Color','b', 'LineWidth',2); hold on
legend1 = legend(subplot(2,2,1),'show');
legend('Approximate x(t)','Exact x(t)')
 
% Plotting the exact solution x'(t) in blue
subplot(2,2,2)
plot(t,dx_exact, 'Color','b', 'LineWidth',2); hold on
legend1 = legend(subplot(2,2,2),'show');
legend('Approximate dx(t)','Exact dx(t)')
 
% Plotting the exact solution x"(t) in blue
subplot(2,2,3)
plot(t,d2x_exact, 'Color','b', 'LineWidth',2); hold on
legend1 = legend(subplot(2,2,3),'show');
legend('Approximate dx(t)','Exact dx(t)')
 
assignin('base', 'x_exact', x_exact);           % Sending x_exact to workspace
assignin('base', 'dx_exact', dx_exact);         % Sending dx_exact to workspace
assignin('base', 'd2x_exact', d2x_exact);       % Sending d2x_exact to workspace
assignin('base', 't', t);                       % Sending t to workspace
clc
 
 
%% =====================End of exact solution============================== 
 
 
%% =========Producing errors between exact x(t) and approximated x(t)======
error_GM=zeros (ne+1,1);
if r1==r2   % For critically damped case
    C=inv([exp(r1*t1)  t1*exp(r2*t1);exp(r1*t2)  t2*exp(r2*t2)])*([x1;x2]);
    C1=C(1);
    C2=C(2);
   
    for j=1:ne+1
        
        error_GM(j)=((C1*exp(r1*t_GM(j))+C2*t_GM(j)*exp(r2*t_GM(j)))-(x_GM(j)))/(C1*exp(r1*t_GM(j))+C2*t_GM(j)*exp(r2*t_GM(j)));
    end
    
else        % For Under-damped or Over-damped case
    
    C=inv([exp(r1*t1) exp(r2*t1);exp(r1*t2) exp(r2*t2)])*([x1;x2]);
    C1=C(1);
    C2=C(2);
    
        
    for j=1:ne+1
        error_GM(j)=((C1*exp(r1*t_GM(j))+C2*exp(r2*t_GM(j)))-(x_GM(j)))/(C1*exp(r1*t_GM(j))+C2*exp(r2*t_GM(j)));
    end
end
 
% Plotting the error of Exact x(t) and approximated x(t)
subplot(2,2,4)
plot(t_GM,error_GM, '-.', 'Color','r', 'LineWidth',2); hold on
xlabel('Time (s)','FontSize',12);
ylabel('Error (%)','FontSize',12);
legend1 = legend(subplot(2,2,4),'show');
legend('Error of approximation')
 
 
assignin('base', 'error_GM', error_GM);     % Sending error to workspace
assignin('base', 't_GM', t_GM);             % Sending t_GM to workspace
assignin('base', 'ne', ne);                 % Sending t_GM to workspace
 
%% ====End of Producing errors between exact x(t) and approximated x(t)====
 
%=============================Program ends=================================


Contact us