%Finite element method code for solving bvp nonlinear ODEs%
% u''+uu'-u=exp(2x), u(0)= 1, u(1)=e %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FEM_Code()
clear all; close all; clc
n=5; % NO of element
nn=n+1; % No of nodes
lgth=1; % Domain length
he=lgth/n; % lenth of each elemnet
x=[0:he:lgth]; % Data point for independant variable
AC=0.00005; % Accuracy
F=zeros(nn,1); % Initialization
F(1)=exp(0); F(nn)=exp(1); % Boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Direct Iterative process to handle nonlinear problem
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c=1.0;
count=0; % Initializations for count for iterations
tic % Time start
while (c>0)
[F1]=assembly(F,n,he);
c=0.0;
for i=1:nn
if (abs(F(i)-F1(i))>AC)
c=c+1;
break;
end
end
F=F1;
count=count+1;
end
disp('Hence solution=:');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output for prinmary and secondary variables %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff=abs(F-exp(x)');
fprintf('No of element=%d\134n',n)
disp(' x FEM Exact Error')
disp([x',F,exp(x)',diff])
fprintf('No of iterations=%d\134n',count)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Ploting of primary variable %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(x,F,'--rs','Linewidth',2)
xlabel('x')
ylabel('u(x)')
title('solution plot to given BVP')
toc % given totlal time
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Derivative of element matrix and Assembly%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [F1]=assembly(F,n,he)
nn=n+1;
k = zeros(nn,nn); % Initialization of main Matrix
R = zeros(nn,1); % Initialization of RHS Matrix
syms x % x as symbolic variable
s=[1-x/he,x/he]; % linear shape function
ds=diff(s,x); % Differentiations of shape function
lmm =[];
for i=1:n
lmm=[lmm;[i,i+1]]; % connectvity Matrix
end
for i=1:n
lm=lmm(i,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Generation of Element Matrix k11 and RHS Matrix f1%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k11=-int(ds'*ds,x,0,he)+(int(s'*ds*s(1),x,0,he)*F(lm(1))...
+int(s'*ds*s(2),x,0,he)*F(lm(2)))-int(s'*s,x,0,he);
f1 = int(exp(2*(x+(i-1)*he))*s',x,0,he);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Assembly accroding to connectivity Matrix%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k(lm,lm) = k(lm,lm) + k11;
R(lm) = R(lm) + f1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Imposing Boundary Conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k(1,:) = 0.0; k(nn,:) = 0.0;
k(1,1) = 1.0; k(nn,nn) = 1.0;
R(1,1) = F(1); R(nn,1) = F(nn);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Solution of equations (F1) %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d = k\134R; % better than using inverse k*R
F1 = d;
end