How to use 'trapz' command in a foor loop
1 view (last 30 days)
Show older comments
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?
13 Comments
Star Strider
on 10 May 2021
In the context of the code, what are ‘1.txt’ and ‘2.txt’?
They are vectors, however everything else seem to be matrices.
Accepted Answer
Star Strider
on 10 May 2021
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
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.
6 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!