I keep receiving an error when I try to run this code for a 3DOF mass-spring system "Index in position 1 exceeds array bounds. Index must not exceed 1".I am confused as to why
1 view (last 30 days)
Show older comments
Fiyinfoluwa Adeseha
on 30 Apr 2022
Commented: Fiyinfoluwa Adeseha
on 30 Apr 2022
I am trying to run this code and everytime I do, I always get an error message that says thus: "Index in position 1 exceeds array bounds. Index must not exceed 1". I am clearly doing something wrong and cannot seem to figure out what it is.
Here is my code below:
clear all;
close all;
n=3; %degree of freedom
m1= 45400; %mass for first floor in (kg)
m2= 45400; %mass for second floor in (kg)
m3= 45400; %mass for third floor in (kg)
k1= 120000; %first steel stiffness in (N/m)
k2= 120000; %second steel stiffness in (N/m)
k3= 120000; %third steel stiffness in (N/m)
c1= 7381.06; %first damper in (N.s/m^2)
c2= 7381.06; %second damper in (N.s/m^2)
c3= 7381.06; %third damper in (N.s/m^2)
M= [m1 0 0;0 m2 0;0 0 m3]; %mass matrix
K= [k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3]; %stiffness matrix
C= [c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3]; % damping matrix
[u,s]=eig(K,M);
u1=u;
s1=s;
x0=[0 0 0]';
xd0=[0 0 0]';
Kh=M^(-0.5)*K*M^(-0.5);
[P,s2]=eig(Kh);
w=sqrt(diag(s2));
S=M^(-0.5)*P;
r0=inv(S)*xd0;
rd0=inv(S)*x0;
%alp=0;
%bet=0.1;
zeta=0.05;
wd=w.*((1-zeta^2)).^0.5;
F=[0 0 840]';
f=P'*M^(-.5)*F;
wr=[0 0 0]';
XX=(f./(((w.^2-wr.^2)+((2*zeta.*w.*wr).^2)).^0.5));
th=atan2((2*zeta.*w.*wr),((w.^2)-(wr.^2)));
phir = atan2(((wd.*(r0-XX.*cos(th)))),(rd0+(zeta.*w.*(r0-XX.*cos(th))))-(wr.*XX.*sin(th)));
Ar=(r0-(XX.*cos(th)))./sin(phir);
t=(0:0.001:10)';
for i=1:n
rt(i,:)=Ar(i,1).*exp(-zeta(i,1)*w(i,1)*t).*sin(wd(i,1)*t+phir(i,1))+XX(i,1)*cos(wr(i,1)*t-th(i,1));
end
%xt=S*rt;
%xt=xt';
figure(1)
subplot(3,1,1)
plot(t,rt(1,:))
title ('Response of Mass #1');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
subplot(3,1,2)
plot(t,rt(1,:))
title ('Response of Mass #2');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
subplot(3,1,3)
plot(t,rt(1,:))
title ('Response of Mass #3');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
0 Comments
Accepted Answer
Cris LaPierre
on 30 Apr 2022
The problem is that zeta is a scalar, not a vector. The fix is to not index this variable.
n=3; %degree of freedom
m1= 45400; %mass for first floor in (kg)
m2= 45400; %mass for second floor in (kg)
m3= 45400; %mass for third floor in (kg)
k1= 120000; %first steel stiffness in (N/m)
k2= 120000; %second steel stiffness in (N/m)
k3= 120000; %third steel stiffness in (N/m)
c1= 7381.06; %first damper in (N.s/m^2)
c2= 7381.06; %second damper in (N.s/m^2)
c3= 7381.06; %third damper in (N.s/m^2)
M= [m1 0 0;0 m2 0;0 0 m3]; %mass matrix
K= [k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3]; %stiffness matrix
C= [c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3]; % damping matrix
[u,s]=eig(K,M);
u1=u;
s1=s;
x0=[0 0 0]';
xd0=[0 0 0]';
Kh=M^(-0.5)*K*M^(-0.5);
[P,s2]=eig(Kh);
w=sqrt(diag(s2));
S=M^(-0.5)*P;
r0=inv(S)*xd0;
rd0=inv(S)*x0;
%alp=0;
%bet=0.1;
zeta=0.05;
wd=w.*((1-zeta^2)).^0.5;
F=[0 0 840]';
f=P'*M^(-.5)*F;
wr=[0 0 0]';
XX=(f./(((w.^2-wr.^2)+((2*zeta.*w.*wr).^2)).^0.5));
th=atan2((2*zeta.*w.*wr),((w.^2)-(wr.^2)));
phir = atan2(((wd.*(r0-XX.*cos(th)))),(rd0+(zeta.*w.*(r0-XX.*cos(th))))-(wr.*XX.*sin(th)));
Ar=(r0-(XX.*cos(th)))./sin(phir);
t=(0:0.001:10)';
for i=1:n
rt(i,:)=Ar(i,1).*exp(-zeta*w(i,1)*t).*sin(wd(i,1)*t+phir(i,1))+XX(i,1)*cos(wr(i,1)*t-th(i,1));
% ^^ removed (i,1)
end
%xt=S*rt;
%xt=xt';
figure(1)
subplot(3,1,1)
plot(t,rt(1,:))
title ('Response of Mass #1');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
subplot(3,1,2)
plot(t,rt(1,:))
title ('Response of Mass #2');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
subplot(3,1,3)
plot(t,rt(1,:))
title ('Response of Mass #3');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
More Answers (0)
See Also
Categories
Find more on Automated Driving Toolbox 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!