MATLAB: Does: “Subscript indices must either be real positive integers or logicals.” error occur? This appears in the equation where the effective load vector is calculated.

numerical integrationnumerical methodvibrationswilson

function [u,v,a]= Wilson(m,c,k,u0,v0,dt,f,wn,zj)
% Input data:
% ----------
% m  = Mass matrix of oscillator [consistent units].
% dt = Time interval [sec].
% u0 = Initial displacement  [Length].
% v0 = Initial velocity [Length/sec].
% f  = Vector with the sampled excitation matrix [consistent units].
% t  = Vector with discrete time matrix [sec].
% c  = Damping coefficient (N.m.s/rad)
% k  = Stiffness coefficient (N.m/rad)
% Output data:
% -----------
% u  = Displacement time history.
% v  = Velocity time history.
% a  = Acceleration time history.
tf = 10 ;                          %Final time
t = 0:dt:tf ;                      %Time intervals 
nt = length(t); 
%Initialize: Matrices for displacement, velocity, acceleration, effective force, displacement at x i+o 
 Z  = zeros(3, nt) ;
fi  = zeros(1,nt)  ;
%Initial conditions
Z(1,1) = u0                                ;
Z(2,1) = v0                                ;
Z(3,1) = ( (f(1)-c*Z(2,1)-k*Z(1,1)) / m )  ;      %30.5081 m/s^2              
%Constants 
o  = 1.4              ;            %Theta, can be changed. 
c0 = 6/(o*dt)^2       ;            %76530612.2449
c1 = 3/(o*dt)         ;            %10714.2857
c2 = 6/o*dt           ;            %18750
c3 = (o*dt)/2         ;            %0.00014
c4 = 6/((o^3)*(dt^2)) ;            %54664723.0321
c5 = 6/(o^2)*dt       ;            %15306.1224
c6 = (1-(3/o))        ;            %-1.1429
c7 = dt/2             ;            %0.0001
c8 = (dt^2)/6         ;            %6.67e-9
%Calculations
for i = 2:nt
%Effective load vector
fi(1,i) = f(1,i-1)+o(f(1,i)-f(1,i-1))+m*(c0*Z(1,i-1)+c2*Z(2,i-1)+2*Z(3,i-1))+c*(c1*Z(1,i-1)+2*Z(2,i-1)+c3*Z(3,i-1)) ;  
%Displacement at time t i+dt
Z(1,i+1) = k/fi(i)                                                                   ;
%Acceleration at time t i+dt
Z(3,i+1) = c4*(Z(1,i+1)-Z(1,i))+c5*Z(2,i)+c6*Z(3,i)                                  ;
%Velocity at time t i+dt
Z(2,i+1) = Z(2,i)+c7*(Z(3,i+1)-Z(3,i))                                               ;
%Displacement t i+dt
Z(1,i+1)= Z(1,i)+dt*Z(2,i)+c8*(Z(3,i+1)+2*Z(3,i))                                    ;
end
u = Z(1,:);
v = Z(2,:);
a = Z(3,:);

Best Answer

You have
fi(1,i) = f(1,i-1)+o(f(1,i)-f(1,i-1))+m*(c0*Z(1,i-1)+c2*Z(2,i-1)+2*Z(3,i-1))+c*(c1*Z(1,i-1)+2*Z(2,i-1)+c3*Z(3,i-1)) ;
Notice the o(f(1,i)-f(1,i-1)) part. f(1,i)-f(1,i-1) is going to be be a floating point value that is not likely to be a positive integer, but you try to index the scalar o there. Perhaps you forgot a multiplication between o and the term after that.
Related Question