Matlab code Shape Functions 1D Lagrangian

From KratosWiki
Jump to: navigation, search
% Shape Functions for 1D problems
% polynomial Ni(x) = a1 + a2x + a3x^2 +...
% Ni(xj)=1 if i=j, Ni(xj)=0 if i<>j 

Np=101;
x=linspace(0,1,Np);

% Linear case
% Two nodes: x1, x2
% Ni(x) = ai1 + ai2•x

xi(1)=0;
xi(2)=1;

Ni(1,:)=shfunc(x,1,xi);
Ni(2,:)=shfunc(x,2,xi);

mdibuja(x,Ni,'x','N1, N2','Linear Shape Funtions'); 


% Quadratic case
% Three nodes: x1, x2, x3
% Ni(x) = ai1 + ai2•x + ai3•x^2

xi(1)=0;
xi(2)=0.5;
xi(3)=1;

for num=1:length(xi)
   Ni(num,:)=shfunc(x,num,xi);
end

mdibuja(x,Ni,'x','N1, N2, N3','Quadratic Shape Funtions'); 


% Cubic case
% Four Nodes: x1, x2, x3, x4
% Ni(x) = ai1 + ai2•x + ai3•x^2 + ai4•x^3

xi(1)=0;
xi(2)=1/3;
xi(3)=2/3;
xi(4)=1;

for num=1:length(xi)
   Ni(num,:)=shfunc(x,num,xi);
end

mdibuja(x,Ni,'x','N1, N2, N3, N4','Cubic Shape Funtions'); 
hold on;
plot(xi,xi*0,'o');


% Five Nodes: x1, x2, x3, x4, x5
% Ni(x) = ai1 + ai2•x + ai3•x^2 + ai4•x^3 + ai5•x^4

xi(1)=0;
xi(2)=1/4;
xi(3)=2/4;
xi(4)=3/4;
xi(5)=1;


for num=1:length(xi)
   Ni(num,:)=shfunc(x,num,xi);
end

mdibuja(x,Ni,'x','N1, N2, N3, N4, N5','Five Node Shape Funtions'); 
hold on;
plot(xi,xi*0,'o');


% Nine Nodes: x1, x2, x3, x4, x5, x6, x7, x8, x9

xi(1)=0;
xi(2)=1/8;
xi(3)=1/4;
xi(4)=3/8;
xi(5)=1/2;
xi(6)=5/8;
xi(7)=3/4;
xi(8)=7/8;
xi(9)=1;



for num=1:length(xi)
   Ni(num,:)=shfunc(x,num,xi);
end

mdibuja(x,Ni,'x','N1, N2, N3, N4, N5, N6, N7, N8, N9','Nine Node Shape Funtions'); 
hold on;
plot(xi,xi*0,'o');
 



% Shape Functions for Lagrange polynomial
% Ni(x)=num/den
% num=(x-x1)*(x-x2)...(x-x{i-1})*(x-x{i+1})...(x-xn)
% den=(xi-x1)*(xi-x2)...(xi-x{i-1})*(xi-x{i+1})...(xi-xn)

function ni=shfunc(x,i,xi)

num=1;
den=1;

for j=1:length(xi)
    if(i~=j)
      num=num.*(x-xi(j));
      den=den.*(xi(i)-xi(j));
    end
end

ni=num./den;
    
 




% plot including labels, grid and ajusting axis...
function mdibuja(x,y,labx,laby,tit)

plot(x,y);
miny=min(1.1*min(min(y)),0.9*min(min(y)));
maxy=max(1.1*max(max(y)),0.9*max(max(y)));
axis([min(x) max(x) floor(miny) ceil(maxy)]);
grid;
xlabel(labx);
ylabel(laby);
title(tit);
 


Download the Source code
Back to One-D Shape Functions
Personal tools
Categories