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 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Fiyinfoluwa Adeseha
le 30 Avr 2022
Commenté : Fiyinfoluwa Adeseha
le 30 Avr 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 commentaires
Réponse acceptée
Cris LaPierre
le 30 Avr 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
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Automated Driving Toolbox dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!