How to Fix: Index Exceeds Matrix Dimensions

I am getting the error, 'Index Exceeds Matrix Dimensions' and I'm not sure why.

Code:

clear
clc
close all

%Constants:
m_b = 0.3; %[kg]
m_f = 0.025; %[kg]
m_r = 0.1; %[kg]
l_f = 0.05; %[m]
l_r = 0.13; %[m]
k_tf = 4; %[N*m/rad]
k_tr = 10; %[N*m/rad]
k_f = 5; %[N/m]
k_r = 20; %[N/m]
zeta_1 = 0.06;
zeta_2 = 0.02;

%Initial Conditions:
y_0 = 0.01; %[m]
y_dot_0 = 0.5; %[m/s]
theta_f_0 = pi/15; %[rad]
theta_f_dot_0 = -0.01; %[rad/s]
theta_r_0 = pi/5; %[rad]
theta_r_dot_0 = 0.4; %[rad/s]

t = 0:0.01:1;
T_f = 0.1*sin(2*pi*t-0.1); %[N*m]
T_r = 200*sin(10*pi*t); %[N*m]

%Matrices:
M = [m_b+m_f+m_r m_f*l_f/2 m_r*l_r/2;m_f*l_f/2 m_f*l_f^2/3 0;m_r*l_r/2 0 m_r*l_r^2/3];
K = 2*[k_f+k_r k_f*l_f k_r*l_r;k_f*l_f k_f*l_f^2+k_tf 0;k_r*l_r 0 k_r*l_r^2+k_tr];
F = [T_f+T_r;T_f;T_r];
zeta = [zeta_1;zeta_2];

%Mass Normalized Stiffness:
K_tilde = inv(sqrtm(M))*K*inv(sqrtm(M));

[P,Lambda]=eigs(K_tilde);

S = inv(sqrtm(M))*P;

for i=1:length(t)
    R(:,i) = (P'*inv(sqrtm(M)))*F(:,i);
end

omega = sqrtm(Lambda);

AB_constants = [1/(2*omega(1,1)) omega(1,1)/2;1/(2*omega(2,2)) omega(2,2)/2];
AB = AB_constants\zeta;

alpha = AB(1);
beta = AB(2);

C = alpha*M + beta*K;

for i=1:length(zeta(:,1))
    omega_d(i) = omega(i,i)*sqrt(1-(zeta(i))^2);
end      

X_0 = [y_0;theta_f_0;theta_r_0];
X_dot_0 = [y_dot_0;theta_f_dot_0;theta_r_dot_0];

r_0 = (S^-1)*X_0;
r_dot_0 = (S^-1)*X_dot_0;

A_ss = [zeros(3,3) eye(3,3);-M^-1*K -M^-1*C];
B_ss = [zeros(3,3);eye(3,3)];
C_ss = [eye(3,3) zeros(3,3)];
D_ss = zeros(3,3);

system = ss(A_ss,B_ss,C_ss,D_ss);

X_matrix = [X_0;X_dot_0];

[Y_SS, T_SS, X_SS] = initial(system,X_matrix,t);

subplot(3,1,1)
plot(T_SS,Y_SS(:,1))
ylabel('y [m]')
hold on
subplot(3,1,2)
plot(T_SS,Y_SS(:,2))
ylabel('\theta_f [rad]')
hold on
subplot(3,1,3)
plot(T_SS,Y_SS(:,3))
xlabel('t [s]')
ylabel('\theta_r [rad]')
hold on

for j=1:length(r_0)
       R_du(j,:) = DampedSystem_ArbitraryForce_DuhamelIntegral_Function(t,t,R(j,:),0,omega(j,j));
end

for j=1:length(r_0)
     r(j,:) = (exp(-zeta(j)*omega(j,j)*t).*(cos(omega_d(j)*t)+(zeta(j)/sqrt(1-(zeta(j))^2))*sin(omega_d(j)*t))*r_0(j)+((1/(omega_d(j)))*exp(-zeta(j)*omega(j)*t).*sin(omega_d(j)*t))*r_dot_0(j)+R_du(j,:));
end

for i=1:length(t)
    X(:,i) = S*r(:,i);
end

figure(1)
subplot(2,1,1)
plot(t,X(1,:))
hold on
subplot(2,1,2)
plot(t,X(2,:))
hold on
legend('State Space Response','Modal Response')

Error:

Index exceeds matrix dimensions.

Error in Homework6Prob2 (line 103)
     r(j,:) =
     (exp(-zeta(j)*omega(j,j)*t).*(cos(omega_d(j)*t)+(zeta(j)/sqrt(1-(zeta(j))^2))*sin(omega_d(j)*t))*r_0(j)+((1/(omega_d(j)))*exp(-zeta(j)*omega(j)*t).*sin(omega_d(j)*t))*r_dot_0(j)+R_du(j,:)); 
>>

Answer 1

[Walter Roberson]

You set zeta = [zeta1;zeta2] when zeta1 and zeta2 are both scalar constants. That leaves zeta as a 2 by 1.

In the line that is having problems, you access zeta(j,j) where j is up to length(t) which is 101, so you try to access up to zeta(101,101)

[Neil Solan]

The section giving the error is as follows:

for j=1:length(r_0)
     r(j,:) = (exp(-zeta(j)*omega(j,j)*t).*(cos(omega_d(j)*t)+(zeta(j)/sqrt(1-(zeta(j))^2))*sin(omega_d(j)*t))*r_0(j)+((1/(omega_d(j)))*exp(-zeta(j)*omega(j)*t).*sin(omega_d(j)*t))*r_dot_0(j)+R_du(j,:));
end

Am I not calling j to be up to length of r_0? How else can I change the length of j? I thought it was supposed to be done in the 1st line of the for loop, as I wrote it?

[Walter Roberson]

Sorry, I miswrote above. But length(r_0) is 3, but there are only two entries zeta.