How to use 'trapz' command in a foor loop

Question

Sumit Saha 2021-5-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/825780-how-to-use-trapz-command-in-a-foor-loop

评论： Star Strider 2021-5-11

采纳的回答： Star Strider

responseSpectra.m

在 MATLAB Online 中打开

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 个评论
显示 11更早的评论隐藏 11更早的评论

Star Strider 2021-5-10

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 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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2021-5-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/825780-how-to-use-trapz-command-in-a-foor-loop#answer_696465

在 MATLAB Online 中打开

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 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2021-5-10

在 MATLAB Online 中打开

‘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 2021-5-10

@Star Strider thanks a lot.

Star Strider 2021-5-11

As always, my pleasure!

请先登录，再进行评论。

How to use 'trapz' command in a foor loop

13 个评论
显示 11更早的评论隐藏 11更早的评论

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to use 'trapz' command in a foor loop

13 个评论 显示 11更早的评论隐藏 11更早的评论

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

13 个评论
显示 11更早的评论隐藏 11更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论