Code covered by the BSD License

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

Dr. Ramin Shamshiri (view profile)

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=================================

```