You'll need to download the following files.
They should be placed in the same directory as the piece of code included below. Good luck!
Disclaimer: No liability is accepted in any event for any damages,
including incidental or consequential damage, lost profits or otherwise
in connection with the use of this program.
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
%
% FEM approximation for the ordinary differential equation
%
% x^2 u'' - 2x u' - 4u = x^2, 10 < x < 20
% u(10) = 0 and u(20) = 100
%
% using 10 linear elements
%
% VARIABLES DESCRIPTION
% k = element matrix
% f = element vector
% kk = system matrix
% ff = system vector
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofs in 'bcdof'
%
%//////////////////////////////////////////////////////////////////////////////////
%
%
% INPUT DATA FOR CONTROL PARAMETERS
%
%
nel=10; % NUMBER OF ELEMENTS
nnel=2; % NUMBER OF NODES PER ELEMENT
ndof=1; % NUMBER OF DOFS PER ELEMENT
nnode=11; % TOTAL NUMBER OF NODES IN THE SYSTEM
sdof=nnode*ndof; % TOTAL SYSTEM DOFS
%
% INPUT DATA FOR NODAL COORDINATE VALUES
%
gcoord(1)=10; gcoord(2)=11; gcoord(3)=12; gcoord(4)=13;
gcoord(5)=14; gcoord(6)=15; gcoord(7)=16; gcoord(8)=17;
gcoord(9)=18; gcoord(10)=19; gcoord(11)=20;
%
% INPUT DATA FOR NODAL CONNECTIVITY FOR EACH ELEMENT
%
nodes(1,1)=1; nodes(1,2)=2; nodes(2,1)=2; nodes(2,2)=3;
nodes(3,1)=3; nodes(3,2)=4; nodes(4,1)=4; nodes(4,2)=5;
nodes(5,1)=5; nodes(5,2)=6; nodes(6,1)=6; nodes(6,2)=7;
nodes(7,1)=7; nodes(7,2)=8; nodes(8,1)=8; nodes(8,2)=9;
nodes(9,1)=9; nodes(9,2)=10; nodes(10,1)=10; nodes(10,2)=11;
%
% INPUT DATA FOR BOUNDARY CONDITIONS
%
bcdof(1)=1; % first node is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=11; % 11th node is constrained
bcval(2)=100; % whose described value is 100
%
% INITIALISATION OF MATRICES AND VECTORS
%
ff=zeros(sdof,1); % initialization of system force vector
kk=zeros(sdof,sdof); % initialization of system matrix
index=zeros(nnel*ndof,1); % initialization of index vector
%
% COMPUTATION OF ELEMENT MATRICES AND VECTORS AND THEIR ASSEMBLY
%
for iel=1:nel % loop for the total number of elements
nl=nodes(iel,1); nr=nodes(iel,2); % extract nodes for (iel)-th element
xl=gcoord(nl); xr=gcoord(nr); % extract nodal coord values for the element
eleng=xr-xl; % element length
index=feeldof1(iel,nnel,ndof); % extract system dofs associated with element
k=feodex2l(xl,xr); % compute element matrix
f=fefx2l(xl,xr); % compute element vector
[kk,ff]=feasmbl2(kk,ff,k,f,index); % assemble element matrices and vectors
end
%
% APPLY BOUNDARY CONDITIONS
%
[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);
%
% SOLVE THE MATRIX EQUATION
%
fsol=kk\ff;
%
% ANALYTICAL SOLUTION
%
esol(1)=0.0;
for i=2:nnode
x=gcoord(i);
esol(i)=0.00102*x^4-0.16667*x^2+64.5187/x;
end
%
% PRINT OUT BOTH THE EXACT AND THE FEM SOLUTIONS
%
num=1:1:sdof;
store=[num' fsol esol']