Matlab code Shape Functions 1D

% 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 + ai2x
% N1(x1) = 1, N1(x2) = 0
% N2(x1) = 0, N2(x2) = 1
% a11 + a12 x1 = 1
% a11 + a12 x2 = 0
% a21 + a22 x1 = 0
% a21 + a22 x2 = 1
% |1   x1    0    0||a11|  |1|
% |1   x2    0    0||a12|  |0|
% |0    0    1   x1||a21|= |0|
% |0    0    1   x2||a22|  |1|
%
% |1   x1 ||a11|= |1|
% |1   x2 ||a12|  |0|
%
% |1   x1||a21|= |0|
% |1   x2||a22|  |1|
%
% a11=x2    a12=-1      a21=-x1     a22=1
%
% det= x2-x1 = l
% N1(x)=(x2-x)/det=(x2-x)/(x2-x1)
% N2(x)=(-x1+x)/det=(x-x1)/(x2-x1)
%

x1=0;
x2=1;

det=x2-x1;
N1=(x2-x)/det;
N2=(x-x1)/det;
Ni=[N1;N2];

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


% Quadratic case
% Three nodes: x1, x2, x3
% Ni(x) = ai1 + ai2•x + ai3•x^2
% N1(x1) = 1, N1(x2) = 0, N1(x3) = 0
% N2(x1) = 0, N2(x2) = 1, N2(x3) = 0
% N3(x1) = 0, N3(x2) = 0, N3(x3) = 1
%
% a11 + a12 x1 + a13 x1^2 = 1
% a11 + a12 x2 + a13 x2^2 = 0
% a11 + a12 x3 + a13 x3^2 = 0
% |1   x1    x1^2||a11|  |1|
% |1   x2    x2^2||a12|= |0|
% |1   x3    x3^2||a13|  |0|
% det = (x2•x3^2 - x3•x2^2) - (x1•x3^2 - x3•x1^2) + (x1•x2^2 - x2•x1^2)
% a11 = x2•x3^2 - x3•x2^2
% a12 = -(x3^2-x2^2)
% a13 = x3 - x2
%
% |1   x1    x1^2||a21|  |0|
% |1   x2    x2^2||a22|= |1|
% |1   x3    x3^2||a23|  |0|
% det = (x2•x3^2 - x3•x2^2) - (x1•x3^2 - x3•x1^2) + (x1•x2^2 - x2•x1^2)
% a21 = -(x1•x3^2 - x3•x1^2)
% a22 = (x3^2-x1^2)
% a23 = -(x3 - x1)
%
% a31 = (x1•x2^2 - x2•x1^2)
% a32 = -(x2^2-x1^2)
% a33 = (x2 - x1)

x1=0;
x2=0.5;
x3=1;

det=(x2*x3^2 - x3*x2^2) - (x1*x3^2 - x3*x1^2) + (x1*x2^2 - x2*x1^2);

a11=x2*x3^2 - x3*x2^2;
a12=-(x3^2-x2^2);
a13=x3 - x2;
N1=(a11 + a12*x + a13* x.^2)/det;

a21=-(x1*x3^2 - x3*x1^2);
a22=(x3^2-x1^2);
a23=-(x3 - x1);
N2=(a21 + a22*x + a23* x.^2)/det;

a31=x1*x2^2 - x2*x1^2;
a32=-(x2^2-x1^2);
a33=x2 - x1;
N3=(a31 + a32*x + a33* x.^2)/det;

Ni=[N1;N2;N3];

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
% N1(x1) = 1, N1(x2) = 0, N1(x3) = 0, N1(x4) = 0
% N2(x1) = 0, N2(x2) = 1, N2(x3) = 0, N2(x4) = 0
% N3(x1) = 0, N3(x2) = 0, N3(x3) = 1, N3(x4) = 0
% N4(x1) = 0, N4(x2) = 0, N4(x3) = 0, N4(x4) = 1
%
% a11 + a12 x1 + a13 x1^2 + a14 x1^3 = 1
% a11 + a12 x2 + a13 x2^2 + a14 x2^3 = 0
% a11 + a12 x3 + a13 x3^2 + a14 x3^3 = 0
% a11 + a12 x4 + a13 x4^2 + a14 x4^3 = 0
% |1   x1    x1^2    x1^3||a11|  |1|
% |1   x2    x2^2    x2^3||a12|= |0|
% |1   x3    x3^2    x3^3||a13|  |0|
% |1   x4    x4^2    x4^3||a13|  |0|
% det = x2*(x3^2*x4^3-x3^3*x4^2)-x3*(x2^2*x4^3-x2^3*x4^2)+x4*(x2^2*x3^3-x2^3*x3^2)
% det = det - (x1*(x3^2*x4^3-x3^3*x4^2)-x3*(x1^2*x4^3-x1^3*x4^2)+x4*(x1^2*x3^3-x1^3*x3^2))
% det = det + (x1*(x2^2*x4^3-x2^3*x4^2)-x2*(x1^2*x4^3-x1^3*x4^2)+x4*(x1^2*x2^3-x1^3*x2^2))
% det = det - (x1*(x2^2*x3^3-x2^3*x3^2)-x2*(x1^2*x3^3-x1^3*x3^2)+x3*(x1^2*x2^3-x1^3*x2^2))
% a11 = x2*(x3^2*x4^3-x3^3*x4^2)-x3*(x2^2*x4^3-x2^3*x4^2)+x4*(x2^2*x3^3-x2^3*x3^2)
% a12 = 
%   |1   x2^2    x2^3|
%  -|1   x3^2    x3^3|
%   |1   x4^2    x4^3|
% a12 = -((x3^2*x4^3-x3^3*x4^2)-(x2^2*x4^3-x2^3*x4^2)+(x2^2*x3^3-x2^3*x3^2))
% a13 =
% |1   x2    x2^3|
% |1   x3    x3^3|
% |1   x4    x4^3|
% a13 = ((x3*x4^3-x3^3*x4)-(x2*x4^3-x2^3*x4)+(x2*x3^3-x2^3*x3))
% a14 =
%  |1   x2    x2^2|
% -|1   x3    x3^2|
%  |1   x4    x4^2|
% a14 =-((x3*x4^2-x3^2*x4)-(x2*x4^2-x2^2*x4)+(x2*x3^2-x2^2*x3))

x1=0;
x2=1/3;
x3=2/3;
x4=1;

det = x2*(x3^2*x4^3-x3^3*x4^2)-x3*(x2^2*x4^3-x2^3*x4^2)+x4*(x2^2*x3^3-x2^3*x3^2);
det = det - (x1*(x3^2*x4^3-x3^3*x4^2)-x3*(x1^2*x4^3-x1^3*x4^2)+x4*(x1^2*x3^3-x1^3*x3^2));
det = det + (x1*(x2^2*x4^3-x2^3*x4^2)-x2*(x1^2*x4^3-x1^3*x4^2)+x4*(x1^2*x2^3-x1^3*x2^2));
det = det - (x1*(x2^2*x3^3-x2^3*x3^2)-x2*(x1^2*x3^3-x1^3*x3^2)+x3*(x1^2*x2^3-x1^3*x2^2));

a11 = x2*(x3^2*x4^3-x3^3*x4^2)-x3*(x2^2*x4^3-x2^3*x4^2)+x4*(x2^2*x3^3-x2^3*x3^2);
a12 = -((x3^2*x4^3-x3^3*x4^2)-(x2^2*x4^3-x2^3*x4^2)+(x2^2*x3^3-x2^3*x3^2));
a13 = ((x3*x4^3-x3^3*x4)-(x2*x4^3-x2^3*x4)+(x2*x3^3-x2^3*x3));
a14 = -((x3*x4^2-x3^2*x4)-(x2*x4^2-x2^2*x4)+(x2*x3^2-x2^2*x3));

N1=(a11 + a12*x + a13* x.^2+a14*x.^3)/det;

a21 = -(x1*(x3^2*x4^3-x3^3*x4^2)-x3*(x1^2*x4^3-x1^3*x4^2)+x4*(x1^2*x3^3-x1^3*x3^2));
a22 = ((x3^2*x4^3-x3^3*x4^2)-(x1^2*x4^3-x1^3*x4^2)+(x1^2*x3^3-x1^3*x3^2));
a23 = -((x3*x4^3-x3^3*x4)-(x1*x4^3-x1^3*x4)+(x1*x3^3-x1^3*x3));
a24 = ((x3*x4^2-x3^2*x4)-(x1*x4^2-x1^2*x4)+(x1*x3^2-x1^2*x3));

N2=(a21 + a22*x + a23* x.^2+a24*x.^3)/det;

a31 = (x1*(x2^2*x4^3-x2^3*x4^2)-x2*(x1^2*x4^3-x1^3*x4^2)+x4*(x1^2*x2^3-x1^3*x2^2));
a32 = -((x2^2*x4^3-x2^3*x4^2)-(x1^2*x4^3-x1^3*x4^2)+(x1^2*x2^3-x1^3*x2^2));
a33 = ((x2*x4^3-x2^3*x4)-(x1*x4^3-x1^3*x4)+(x1*x2^3-x1^3*x2));
a34 = -((x2*x4^2-x2^2*x4)-(x1*x4^2-x1^2*x4)+(x1*x2^2-x1^2*x2));

N3=(a31 + a32*x + a33* x.^2+a34*x.^3)/det;

a41 = -(x1*(x2^2*x3^3-x2^3*x3^2)-x2*(x1^2*x3^3-x1^3*x3^2)+x3*(x1^2*x2^3-x1^3*x2^2));
a42 = ((x2^2*x3^3-x2^3*x3^2)-(x1^2*x3^3-x1^3*x3^2)+(x1^2*x2^3-x1^3*x2^2));
a43 = -((x2*x3^3-x2^3*x3)-(x1*x3^3-x1^3*x3)+(x1*x2^3-x1^3*x2));
a44 = ((x2*x3^2-x2^2*x3)-(x1*x3^2-x1^2*x3)+(x1*x2^2-x1^2*x2));

N4=(a41 + a42*x + a43* x.^2+a44*x.^3)/det;

Ni=[N1;N2;N3;N4];

mdibuja(x,Ni,'x','N1, N2, N3, N4','Cubic Shape Funtions'); 
nodes=[0,1/3,2/3,0];
hold on;
plot(nodes,nodes*0,'o');
 

% 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