How to use 'trapz' command in a foor loop

main file:

% Analysis of response spectra for given ground motion acceleration data
clear all
close all
clc
% Supershear & Subshear Earthquake Data
load Histories/IND.X0kY-1k.CXY.sema;% Ground-motion acceleration data 
load Histories/IND.X0kY-2k.CXY.sema;% Ground-motion acceleration data
ground_acc(:,:,1) = IND_X0kY_1k_CXY;
ground_acc(:,:,2) = IND_X0kY_2k_CXY;

T_D = 8 ; % duration of earthquake in secs
dt = 0.002; % dt = sampling interval in seconds
xi = 0.02;  % xi = ratio of critical damping
T_series_EQ = [0.002:dt:T_D]'; % Time column for EQ data
g =9.81 ; % accleration due to gravity in m/sec^2

% time period in vector form for response spectrum plot
Time_Period_Spectrum = [0.01:0.01:10];
for j=1:2
[PSA(:,:,j), PSV(:,:,j), SD(:,:,j), SA(:,:,j), SV(:,:,j), OUT(:,:,j),vel(:,:,j)] = responseSpectra(xi, Time_Period_Spectrum, ground_acc(:,:,j), dt);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Different Intensity Measure Calculation for EQ DATA %%%%%%%%%%%%%%%%
%%% Reference : Ground motion intensity measures for seismic probabilistic risk analysis
%% Peak Based Intensity Measure Calculation
for p=1:2
Peak_Ground_Velocity(p) = max(abs(vel(:,:,p)));
Peak_Ground_Accleration(p) = max(abs(ground_acc(:,:,p)));
%%% Pseudo-spectral accleration at fundamental period
T_f_structure = 1.42; % Fundamental Period of structure  in sec 
PSA_TF(p) = interp1 (Time_Period_Spectrum,PSA(:,:,p),T_f_structure); 
%%% Spectral Intensity (related to kinetic energy stored in the structure During Earthquake)
S_v(:,p) = PSV(:,:,p)';
T = Time_Period_Spectrum';
Trange_SI = (T >= 0.1) & (T <= 2.5); %-period interval of Structures 
I_H(p) = trapz(T(Trange_SI), S_v(:,p,(Trange_SI))); %Spectral_Intensity 
end 

Function file: given below

How to use Trapz command in case of loop for multiple eq data?

[Star Strider]

That appears to me to be correct, and it is apparently not throwing errors.

How do you want to change it?

[Sumit Saha]

@Star Strider i'm facing the problem with I_H part

[Star Strider]

What problem?

[Sumit Saha]

@Star Strider Index exceeds matrix dimensions.

Error in Different_IM_Calculation (line 35)

I_H(p) = trapz(T(Trange_SI), S_v(:,p,(Trange_SI))); %Spectral_Intensity

[Sumit Saha]

@Star Strider S_v(:,p,(Trange_SI)) ...this is the problematic area

[Star Strider]

‘Index exceeds matrix dimensions.’

That could be due either to the magnitude of ‘p’ or to ‘Trang_SI’. It is attempting to reference elements of ‘S_v’ that are greater than the dimensions of that array.

For example —

A = [1 2 3; 2 1 0];
B = A(3,4)

Index in position 1 exceeds array bounds (must not exceed 2).

[Sumit Saha]

@Star Strider p basically 1,2 and T is time arraged in column vector. Ignore my last two lines. I just want to extract Trang_SI as T >= 0.1 & T<=2.5 and that Trang_SI should be arranged in column matrix. S_v is the spectral velocity column vector related to T but during integration I need those values lying in between T>=0.1 & T<=2.5. (S_v can be represented in two columns and one columns represents time and other column represents actual S-v values). This intergration should be completed for this for loop. I hope you've understood my point.

[Sumit Saha]

@Star Strider I've attached one pic for your reference

[Star Strider]

I understand what you want to do (at least as best as I can), however the issue of whether the loigical vector reference ‘Trang_SI’ or ‘p’ is throwing the dimension error is still not resolved (at least as far as I can understand this).

Prehaps knowing the dimensions of ‘S_v’ would help.

[Sumit Saha]

@Star Strider Can you please run the matlab script...I'm attaching text files. Kindly run once. What can be the otherway to solve the integration?

[Sumit Saha]

@Star Strider earlier I had attached the matlab scripts

[Star Strider]

In the context of the code, what are ‘1.txt’ and ‘2.txt’?

They are vectors, however everything else seem to be matrices.

[Sumit Saha]

@Star delete these two lines...write load 1.txt && load 2.txt

% Supershear & Subshear Earthquake Data

load Histories/IND.X0kY-1k.CXY.sema;% Ground-motion acceleration data

load Histories/IND.X0kY-2k.CXY.sema;% Ground-motion acceleration data

Answer 1

[Star Strider]

The problem is that ‘S_v’ has only one non-singleton dimension (i.e. it is a vector).

That is throwing the error.

Changing the assignment to —

I_H(p) = trapz(T(Trange_SI), S_v((Trange_SI))); %Spectral_Intensity

no longer throws the error, however the values fo ‘I_H’ aare both identical.

I have no idea what you are doing, so I will let you sort that out.

I ran this in the Online Run Feature—

ground_acc(:,:,1) = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/613135/1.txt');
ground_acc(:,:,2) = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/613140/2.txt');
T_D = 8 ; % duration of earthquake in secs
dt = 0.002; % dt = sampling interval in seconds
xi = 0.02;  % xi = ratio of critical damping
T_series_EQ = [0.002:dt:T_D]'; % Time column for EQ data
g =9.81 ; % accleration due to gravity in m/sec^2
% time period in vector form for response spectrum plot
Time_Period_Spectrum = [0.01:0.01:10];
for j=1:2
[PSA(:,:,j), PSV(:,:,j), SD(:,:,j), SA(:,:,j), SV(:,:,j), OUT(:,:,j),vel(:,:,j)] = responseSpectra(xi, Time_Period_Spectrum, ground_acc(:,:,j), dt);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Different Intensity Measure Calculation for EQ DATA %%%%%%%%%%%%%%%%
%%% Reference : Ground motion intensity measures for seismic probabilistic risk analysis
%% Peak Based Intensity Measure Calculation
for p=1:2
Peak_Ground_Velocity(p) = max(abs(vel(:,:,p)));
Peak_Ground_Accleration(p) = max(abs(ground_acc(:,:,p)));
%%% Pseudo-spectral accleration at fundamental period
T_f_structure = 1.42; % Fundamental Period of structure  in sec 
PSA_TF(p) = interp1 (Time_Period_Spectrum,PSA(:,:,p),T_f_structure); 
%%% Spectral Intensity (related to kinetic energy stored in the structure During Earthquake)
S_v(:,p) = PSV(:,:,p)';
T = Time_Period_Spectrum';
Trange_SI = (T >= 0.1) & (T <= 2.5); %-period interval of Structures 
% Qp = p
% QTrange_SI = Trange_SI
% QsizeT = size(T)
% QsizeS_v = size(S_v)
% QT = T(Trange_SI)
% % QS_v = S_v(:,p,(Trange_SI))
% QS_v = S_v((Trange_SI))
% % I_H(p) = trapz(T(Trange_SI), S_v(:,p,(Trange_SI))); %Spectral_Intensity 
I_H(p) = trapz(T(Trange_SI), S_v((Trange_SI))); %Spectral_Intensity 
QI_H = sprintf('I_H(%d) = %23.15E',p,I_H(p))
end 

QI_H = 'I_H(1) =   4.013137327351848E-01'

QI_H = 'I_H(2) =   4.013137327351848E-01'

function [PSA, PSV, SD, SA, SV, OUT,vel] = responseSpectra(xi, Time_Period_Spectrum, ground_acc, dt)
vel = cumtrapz(ground_acc)*dt;
disp = cumtrapz(vel)*dt;

% Spectral solution
for i = 1:length(Time_Period_Spectrum)
    omega_w = 2*pi/Time_Period_Spectrum(i);
    C = 2*xi*omega_w;
    K = omega_w^2;
    y(:,1) = [0;0];
    A = [0 1; -K -C]; 
    Ae = expm(A*dt); AeB = A\(Ae-eye(2))*[0;1];
    
    for k = 2:numel(ground_acc)
        y(:,k) = Ae*y(:,k-1) + AeB*ground_acc(k);
    end
    
    displ = (y(1,:))';                          % Relative displacement vector (cm)
    veloc = (y(2,:))';                          % Relative velocity (cm/s)
    foverm = omega_w^2*displ;                    % Lateral resisting force over mass (cm/s^2)
    absacc = -2*xi*omega_w*veloc-foverm;         % Absolute acceleration from equilibrium (cm/s^2)
    
    % Extract peak values
    displ_max(i) = max(abs(displ));             % Spectral relative displacement (cm)
    veloc_max(i) = max(abs(veloc));             % Spectral relative velocity (cm/s)
    absacc_max(i) = max(abs(absacc));           % Spectral absolute acceleration (cm/s^2)
    
    foverm_max(i) = max(abs(foverm));           % Spectral value of lateral resisting force over mass (cm/s^2)
    pseudo_acc_max(i) = displ_max(i)*omega_w^2;  % Pseudo spectral acceleration (cm/s)
    pseudo_veloc_max(i) = displ_max(i)*omega_w;  % Pseudo spectral velocity (cm/s)
    
    PSA(i) = pseudo_acc_max(i);                 % PSA (cm/s^2)
    SA(i)  = absacc_max(i);                     % SA (cm/s^2)
    PSV(i) = pseudo_veloc_max(i);               % PSV (cm/s)
    SV(i)  = veloc_max(i);                      % SV (cm/s)
    SD(i)  = displ_max(i);                      % SD  (cm)
   
    % Time series of acceleration, velocity and displacement response of SDF system
    OUT.acc(:,i) = absacc;
    OUT.vel(:,i) = veloc;
    OUT.disp(:,i) = displ;
end
end

I kept the debugging steps in the code, although commented-out. Delete them if you no longer need them, or un-comment or change them if there are still problems and you want to see the interim results.

[Sumit Saha]

@Star Strider S_v (two column vectors) basically corresponding to two eqs but in this line S_v has not been specified. I think that is why it executes with the first column of S_v.

I_H(p) = trapz(T(Trange_SI), S_v((Trange_SI))); %Spectral_Intensity

[Star Strider]

I solved the original problem you requested help with, and explained it.

Since ‘S_v’ is a (1000x1) vector, it has only one column.

Now that the code runs, you need to correct any remaining problems.

[Sumit Saha]

@Star Strider S_v has two columns right. From the begining I'm telling you this...I want to run it for two or multiple files. You've solved the problem but it executes with the one column of S_v. Just tell me when p=2 how will I use data from S_v(:,2) in that I_H(p) line ?

[Star Strider]

‘S_v has two columns right.’

It didn’t previously It does now.

Try this —

ground_acc(:,:,1) = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/613135/1.txt');
ground_acc(:,:,2) = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/613140/2.txt');
T_D = 8 ; % duration of earthquake in secs
dt = 0.002; % dt = sampling interval in seconds
xi = 0.02;  % xi = ratio of critical damping
T_series_EQ = [0.002:dt:T_D]'; % Time column for EQ data
g =9.81 ; % accleration due to gravity in m/sec^2
% time period in vector form for response spectrum plot
Time_Period_Spectrum = [0.01:0.01:10];
for j=1:2
[PSA(:,:,j), PSV(:,:,j), SD(:,:,j), SA(:,:,j), SV(:,:,j), OUT(:,:,j),vel(:,:,j)] = responseSpectra(xi, Time_Period_Spectrum, ground_acc(:,:,j), dt);
end

PSV_Size = size(PSV)

PSV_Size = 1×3

1 1000 2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Different Intensity Measure Calculation for EQ DATA %%%%%%%%%%%%%%%%
%%% Reference : Ground motion intensity measures for seismic probabilistic risk analysis
%% Peak Based Intensity Measure Calculation
for p=1:2
Peak_Ground_Velocity(p) = max(abs(vel(:,:,p)));
Peak_Ground_Accleration(p) = max(abs(ground_acc(:,:,p)));
%%% Pseudo-spectral accleration at fundamental period
T_f_structure = 1.42; % Fundamental Period of structure  in sec 
PSA_TF(p) = interp1 (Time_Period_Spectrum,PSA(:,:,p),T_f_structure); 
%%% Spectral Intensity (related to kinetic energy stored in the structure During Earthquake)
S_v(:,p) = PSV(1,:,p);                                              % <<== FIXED SUBSCRIPT REFERENCE
T = Time_Period_Spectrum';
Trange_SI = (T >= 0.1) & (T <= 2.5); %-period interval of Structures 
% Qp = p
% QTrange_SI = Trange_SI
% QsizeT = size(T)
% QsizeS_v = size(S_v)
% QT = T(Trange_SI)
% % QS_v = S_v(:,p,(Trange_SI))
% QS_v = S_v((Trange_SI))
% % I_H(p) = trapz(T(Trange_SI), S_v(:,p,(Trange_SI))); %Spectral_Intensity 
I_H(p) = trapz(T(Trange_SI), S_v(Trange_SI,p)); %Spectral_Intensity % <<== FIXED SUBSCRIPT REFERENCE
QI_H = sprintf('I_H(%d) = %23.15E',p,I_H(p))
end 

QI_H = 'I_H(1) =   4.013137327351848E-01'

QI_H = 'I_H(2) =   3.220060876654621E-01'

function [PSA, PSV, SD, SA, SV, OUT,vel] = responseSpectra(xi, Time_Period_Spectrum, ground_acc, dt)
vel = cumtrapz(ground_acc)*dt;
disp = cumtrapz(vel)*dt;

% Spectral solution
for i = 1:length(Time_Period_Spectrum)
    omega_w = 2*pi/Time_Period_Spectrum(i);
    C = 2*xi*omega_w;
    K = omega_w^2;
    y(:,1) = [0;0];
    A = [0 1; -K -C]; 
    Ae = expm(A*dt); AeB = A\(Ae-eye(2))*[0;1];
    
    for k = 2:numel(ground_acc)
        y(:,k) = Ae*y(:,k-1) + AeB*ground_acc(k);
    end
    
    displ = (y(1,:))';                          % Relative displacement vector (cm)
    veloc = (y(2,:))';                          % Relative velocity (cm/s)
    foverm = omega_w^2*displ;                    % Lateral resisting force over mass (cm/s^2)
    absacc = -2*xi*omega_w*veloc-foverm;         % Absolute acceleration from equilibrium (cm/s^2)
    
    % Extract peak values
    displ_max(i) = max(abs(displ));             % Spectral relative displacement (cm)
    veloc_max(i) = max(abs(veloc));             % Spectral relative velocity (cm/s)
    absacc_max(i) = max(abs(absacc));           % Spectral absolute acceleration (cm/s^2)
    
    foverm_max(i) = max(abs(foverm));           % Spectral value of lateral resisting force over mass (cm/s^2)
    pseudo_acc_max(i) = displ_max(i)*omega_w^2;  % Pseudo spectral acceleration (cm/s)
    pseudo_veloc_max(i) = displ_max(i)*omega_w;  % Pseudo spectral velocity (cm/s)
    
    PSA(i) = pseudo_acc_max(i);                 % PSA (cm/s^2)
    SA(i)  = absacc_max(i);                     % SA (cm/s^2)
    PSV(i) = pseudo_veloc_max(i);               % PSV (cm/s)
    SV(i)  = veloc_max(i);                      % SV (cm/s)
    SD(i)  = displ_max(i);                      % SD  (cm)
   
    % Time series of acceleration, velocity and displacement response of SDF system
    OUT.acc(:,i) = absacc;
    OUT.vel(:,i) = veloc;
    OUT.disp(:,i) = displ;
end
end

.

[Sumit Saha]

@Star Strider thanks a lot.

[Star Strider]

As always, my pleasure!