MATLAB Answers

How to use 'trapz' command in a foor loop

조회 수: 6(최근 30일)
Sumit Saha
Sumit Saha 2021년 5월 9일
댓글: Star Strider 2021년 5월 11일
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
Sumit Saha
Sumit Saha 2021년 5월 10일
@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

댓글을 달려면 로그인하십시오.

채택된 답변

Star Strider
Star Strider 2021년 5월 10일
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.
  댓글 수: 6
Star Strider
Star Strider 2021년 5월 11일
As always, my pleasure!

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by