plotting error, different curve for same code
2 views (last 30 days)
Show older comments
rakesh kumar
on 29 Jan 2022
Answered: Priyanka Kondapalli
on 2 Feb 2022
i am getting two different curves for code attached below. the two figures that I have attached are for the same code, please help why am i getting this error, in actual practice I should get only one type of curve as shown in figure(i).
clc;
clear all;
% EX891.m: MATLAB program to solve a static beam deflection using
% Hermitian beam elements
% Variable descriptions
% ? = element stiffness matrix
% kk = system stiffness matrix
% ff = system force 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 dofe in bcdof
nel=5; % number of elements
nnel=2; % number of nodes per element
ndof=2; % number of dofs per node
nnode=(nnel-1)*nel+1; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
el=1; % elastic modulus
xi=1; % moment of inertia of cross-section
tleng=1; % length of a half of the beam
leng=tleng/nel; % element length of equal size
area=1; % cross-sectional area of the beam
rho=1; % mass density (arbitrary value for this problem because
% it is not used for the static problem)
ipt=1; % option for mass matrix (arbitrary value and not used here)
bcdof(1)=1; % first dof (deflection at left end) is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=2; % 12th dof (slope at the symmetric end) is constrained
bcval(2)=0; % whose described value is 0
displacement=zeros(sdof/2,1);
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
ff(2*nel+1)=-1; % because a half of the load is applied due to symmetry
v=1:nel;
N=nel;
p=permn(v,N);
[l r]=size(p);
mm=zeros(sdof,sdof);
force=zeros(sdof, 1);
index=zeros(nnel*ndof,1);
%%
%%
i=1:l;
kk(:,:,i)=zeros(sdof, sdof,l);
mm(:,:,i)=zeros(sdof, sdof,l);
for l=1:l
for iel=1:nel % loop for the total number of elements
edof = nnel*ndof;
start = (iel-1)*(nnel-1)*ndof;
%
for i=1:edof
index(i)=start+i;
end
sigma=0.55;
tleng=1;
d=tleng/nel;
i=1:nel;
j=1:nel;
ee=zeros(nel,nel);
R=zeros(nel,nel);
e(i)=tleng/2*nel+(tleng/nel)*(i-1);
for i=1:nel
for j=1:nel
ee(i,j)=e(i)-e(j);
R(i,j)=sigma^2*exp(-abs(ee(i,j)/d));
end
end
C1=chol(R,'lower');
P1=normrnd(0,1,[1,50]);
Q1=[0.18066 0.83684 -1.4424 1.3025 -0.48065 1.3233]
R1=rand(1,nel);
Z1=R1(1:nel);
alpha=C1*transpose(Z1);
E0=1;
el4=E0*(1+alpha);
el2= 1 + zeros( 1, 100 );
el=el4(1:nel);
x=p(l,iel);
c(iel)=el(iel)*xi/(leng^3);
k(:,:,iel)=c(iel)*[12 6*leng -12 6*leng; 6*leng 4*leng^2 -6*leng 2*leng^2; -12 -6*leng 12 -6*leng;6*leng 2*leng^2 -6*leng 4*leng^2];
edof = length (index);
for i=1:edof
ii=index(i) ;
for j=1:edof
jj=index(j) ;
kk(ii,jj,l)=kk(ii,jj,l)+k(i,j,x);
end
end
if ipt==1
%
c2(iel)=el2(iel)*rho*area*leng/420 ;
m(:,:,l)=c2(iel)*[156 22*leng 54 -13*leng;...
22*leng 4*leng^2 13*leng -3*leng^2;...
54 13*leng 156 -22*leng;...
-13*leng -3*leng^2 -22*leng 4*leng^2];
%
% lumped mass matrix
%
elseif ipt==2
m=zeros(4,4);
mass=rho*area*leng;
m=diag([mass/2 0 mass/2 0]);
%
% diagonal mass matrix
%
else
%
m=zeros(4,4);
mass=rho*area*leng;
m=mass*diag([l/2 leng^2/78 1/2 leng^2/78]);
%
end
edof = length (index);
for i=1:edof
ii=index(i) ;
for j=1:edof
jj=index(j) ;
mm(ii,jj,l)=mm(ii,jj,l)+m(i,j,x);
end
end
n=length(bcdof);
sdof=size(kk);
for i=1:n
c=bcdof(i);
for j=1:sdof
kk(c,j,l)=0;
kk(j,c,l)=0;
mm(c,j,l)=0;
mm(j,c,l)=0;
end
%
mm(c,c,l)=1;
end
end
kkk(:,:,l)=kk(3:end,3:end,l);
mmm(:,:,l)=mm(3:end,3:end,l);
[V(:,:,l) D(:,:,l)]=eig(kkk(:,:,l),mmm(:,:,l));
omega_1(:,:,l)=diag(D(:,:,l));
omega(:,:,l)=omega_1(:,:,l).^0.5;
f(:,:,l)=omega(:,:,l)/(2*pi);
v1(:,:,l)=V(:,[1 2 3],l);
end
g1=reshape(v1,2*nel,[]);
b=zeros(2,size(g1,2));
g=([b;g1]);
x=linspace(0,1);
x1=0;
x2=1/nel;
x3=1;
L1=(x.^2-x*(x2+x3)+x2*x3)/((x1-x2)*(x1-x3));
L2=(x.^2-x*(x1+x3)+x1*x3)/((x2-x1)*(x2-x3));
L3=(x.^2-x*(x1+x2)+x1*x2)/((x3-x1)*(x3-x2));
L1d=(2*x-(x2+x3))/((x1-x2)*(x1-x3));
L2d=(2*x-(x1+x3))/((x2-x1)*(x2-x3));
L3d=(2*x-(x1+x2))/((x3-x1)*(x3-x2));
L1da=(2*x1-(x2+x3))/((x1-x2)*(x1-x3));
L2db=(2*x2-(x1+x3))/((x2-x1)*(x2-x3));
L3dc=(2*x3-(x1+x2))/((x3-x1)*(x3-x2));
h1=(1-2*(x-x1)*L1da).*L1.*L1;
h2=(x-x1).*L1.*L1;
h3=(1-2*(x-x2)*L2db).*L2.*L2;
h4=(x-x2).*L2.*L2;
h5=(1-2*(x-x3)*L3dc).*L3.*L3;
h6=(x-x3).*L3.*L3;
for m=3:3:100
u=h1*g(1,m)+h2*g(2,m)+h3*g(3,m)+h4*g(4,m)+h5*g((end-1),m)+h6*g(end,m);
plot(x,u)
hold on
v=zeros(length(u),1);
plot(x,v,'r')
title('Mode shapes of the Cantilever beam ;Number of element=5')
end
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/877665/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/877670/image.jpeg)
figure(i) figure(ii)
1 Comment
DGM
on 29 Jan 2022
This is not a plotting error. The values in V(:,1:3,x) get flipped for some particular cases of x. Why the eigenvectors are different, I don't know. I don't see any dramatic changes to the numbers further up the code. Currently, you're the expert in your own calculations.
To anyone else trying to help with troubleshooting, you'll need permn() from the FEX:
Accepted Answer
Priyanka Kondapalli
on 2 Feb 2022
Hi,
Please execute the below code where the last section of the code is out side the loop
clc;
clear all;
% EX891.m: MATLAB program to solve a static beam deflection using
% Hermitian beam elements
% Variable descriptions
% ? = element stiffness matrix
% kk = system stiffness matrix
% ff = system force 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 dofe in bcdof
nel=5; % number of elements
nnel=2; % number of nodes per element
ndof=2; % number of dofs per node
nnode=(nnel-1)*nel+1; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
el=1; % elastic modulus
xi=1; % moment of inertia of cross-section
tleng=1; % length of a half of the beam
leng=tleng/nel; % element length of equal size
area=1; % cross-sectional area of the beam
rho=1; % mass density (arbitrary value for this problem because
% it is not used for the static problem)
ipt=1; % option for mass matrix (arbitrary value and not used here)
bcdof(1)=1; % first dof (deflection at left end) is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=2; % 12th dof (slope at the symmetric end) is constrained
bcval(2)=0; % whose described value is 0
displacement=zeros(sdof/2,1);
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
ff(2*nel+1)=-1; % because a half of the load is applied due to symmetry
v=1:nel;
N=nel;
p=permn(v,N);
[l r]=size(p);
mm=zeros(sdof,sdof);
force=zeros(sdof, 1);
index=zeros(nnel*ndof,1);
%%
%%
i=1:l;
kk(:,:,i)=zeros(sdof, sdof,l);
mm(:,:,i)=zeros(sdof, sdof,l);
for l=1:l
for iel=1:nel % loop for the total number of elements
edof = nnel*ndof;
start = (iel-1)*(nnel-1)*ndof;
%
for i=1:edof
index(i)=start+i;
end
sigma=0.55;
tleng=1;
d=tleng/nel;
i=1:nel;
j=1:nel;
ee=zeros(nel,nel);
R=zeros(nel,nel);
e(i)=tleng/2*nel+(tleng/nel)*(i-1);
for i=1:nel
for j=1:nel
ee(i,j)=e(i)-e(j);
R(i,j)=sigma^2*exp(-abs(ee(i,j)/d));
end
end
C1=chol(R,'lower');
P1=normrnd(0,1,[1,50]);
Q1=[0.18066 0.83684 -1.4424 1.3025 -0.48065 1.3233]
R1=rand(1,nel);
Z1=R1(1:nel);
alpha=C1*transpose(Z1);
E0=1;
el4=E0*(1+alpha);
el2= 1 + zeros( 1, 100 );
el=el4(1:nel);
x=p(l,iel);
c(iel)=el(iel)*xi/(leng^3);
k(:,:,iel)=c(iel)*[12 6*leng -12 6*leng; 6*leng 4*leng^2 -6*leng 2*leng^2; -12 -6*leng 12 -6*leng;6*leng 2*leng^2 -6*leng 4*leng^2];
edof = length (index);
for i=1:edof
ii=index(i) ;
for j=1:edof
jj=index(j) ;
kk(ii,jj,l)=kk(ii,jj,l)+k(i,j,x);
end
end
if ipt==1
%
c2(iel)=el2(iel)*rho*area*leng/420 ;
m(:,:,l)=c2(iel)*[156 22*leng 54 -13*leng;...
22*leng 4*leng^2 13*leng -3*leng^2;...
54 13*leng 156 -22*leng;...
-13*leng -3*leng^2 -22*leng 4*leng^2];
%
% lumped mass matrix
%
elseif ipt==2
m=zeros(4,4);
mass=rho*area*leng;
m=diag([mass/2 0 mass/2 0]);
%
% diagonal mass matrix
%
else
%
m=zeros(4,4);
mass=rho*area*leng;
m=mass*diag([l/2 leng^2/78 1/2 leng^2/78]);
%
end
edof = length (index);
for i=1:edof
ii=index(i) ;
for j=1:edof
jj=index(j) ;
mm(ii,jj,l)=mm(ii,jj,l)+m(i,j,x);
end
end
n=length(bcdof);
sdof=size(kk);
for i=1:n
c=bcdof(i);
for j=1:sdof
kk(c,j,l)=0;
kk(j,c,l)=0;
mm(c,j,l)=0;
mm(j,c,l)=0;
end
%
mm(c,c,l)=1;
end
end
kkk(:,:,l)=kk(3:end,3:end,l);
mmm(:,:,l)=mm(3:end,3:end,l);
[V(:,:,l) D(:,:,l)]=eig(kkk(:,:,l),mmm(:,:,l));
omega_1(:,:,l)=diag(D(:,:,l));
omega(:,:,l)=omega_1(:,:,l).^0.5;
f(:,:,l)=omega(:,:,l)/(2*pi);
v1(:,:,l)=V(:,[1 2 3],l);
end
g1=reshape(v1,2*nel,[]);
b=zeros(2,size(g1,2));
g=([b;g1]);
x=linspace(0,1);
x1=0;
x2=1/nel;
x3=1;
L1=(x.^2-x*(x2+x3)+x2*x3)/((x1-x2)*(x1-x3));
L2=(x.^2-x*(x1+x3)+x1*x3)/((x2-x1)*(x2-x3));
L3=(x.^2-x*(x1+x2)+x1*x2)/((x3-x1)*(x3-x2));
L1d=(2*x-(x2+x3))/((x1-x2)*(x1-x3));
L2d=(2*x-(x1+x3))/((x2-x1)*(x2-x3));
L3d=(2*x-(x1+x2))/((x3-x1)*(x3-x2));
L1da=(2*x1-(x2+x3))/((x1-x2)*(x1-x3));
L2db=(2*x2-(x1+x3))/((x2-x1)*(x2-x3));
L3dc=(2*x3-(x1+x2))/((x3-x1)*(x3-x2));
h1=(1-2*(x-x1)*L1da).*L1.*L1;
h2=(x-x1).*L1.*L1;
h3=(1-2*(x-x2)*L2db).*L2.*L2;
h4=(x-x2).*L2.*L2;
h5=(1-2*(x-x3)*L3dc).*L3.*L3;
h6=(x-x3).*L3.*L3;
for m=3:3:100
u=h1*g(1,m)+h2*g(2,m)+h3*g(3,m)+h4*g(4,m)+h5*g((end-1),m)+h6*g(end,m);
plot(x,u)
hold on
title('Mode shapes of the Cantilever beam ;Number of element=5')
end
v=zeros(length(u),1);
plot(x,v,'r')
0 Comments
More Answers (0)
See Also
Categories
Find more on Dynamic System Models in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!