# storing variables after a loop (easy)

6 ビュー (過去 30 日間)
Paul Rogers 2021 年 1 月 8 日
コメント済み: Paul Rogers 2021 年 1 月 11 日
Hi, I managed to make a loop in wich I run a function.
The expected outputs are:
t,y(:,1),y(:,2)
The problem is that at the end of the iteration I only get the last part of the loop, but
I would llike to have all the datas, for those vectors, computed during the looop.
Here's the part of the code:
for i=1:(length(giri)-2)
U1(i) = (D1*pi.*giri(i))/60;
U2(i) = (D2*pi.*giri(i))/60;
B(i) = U1(i)./(2.*wh.*Lc);
[t,y]=ode113(@greitzer,[t0(i),tf(i)],[y1(i),y2(i)],options,U1(i),U2(i),B(i)); %I found ode113 is way more efficient
t0(i+1) = time_rpm(i+1); %simulation's start [s]
tf(i+1) = time_rpm(i+2); %simulation's finish [s]
y1(i+1) = y(end,1); %initial condition 1
y2(i+1) = y(end,2); %initial condition 2
plot(t,y(:,2))
grid on
grid minor
xlabel('t [s]')
ylabel('\Psi')
hold on
end
I'd like to call:
t=time
y(:,1)=phi
y(:,2)=psi
##### 8 件のコメント6 件の古いコメントを表示6 件の古いコメントを非表示
Mathieu NOE 2021 年 1 月 11 日
sure
I put the same code again below so you can accept it
good luck for the future
clear
clc
close all
%% Gas properties
a01 = 340; %speed of sound in gas(m/s)
p01 = 1e5; %Compressor inlet pressure(Pa)
T01 = 303.35; %Compressor inlet temperature(K)
cp = 1005; %specific heat capacity for constant pressure, property of gas(J/(kg*K))
ro1 = 1.15; %Gas density(kg/m^3)
Re = 100000; %Reynolds number()
k = 1.4; %Gas constant ratio, cp/cv()
%% Compressor characteristics
Lc = 0.53; %Length of compressor and duct(m)
Vp = 0.01; %Plenum volume(m^3)
Ac = 0.0026;
% kl = 0.0008; %Throttle gain, proportional to the throttle opening
Dt1 = 0.074; %Impeller diameter at inducer tip(m)
Dh1 = 0.032; %Impeller diameter at hub casing(m)
D1 = 1/sqrt(2)*sqrt((Dt1)^2+(Dh1)^2); %Average diameter(m)
Di = 0.02; %Mean hydraulic diameter of impeller(m)
Dd = 0.02; %Mean hydraulic diameter of impeller(m)
I = 0.001; %Moment of inertia for the compressor spool(kg*m^2)
alfa1 = pi/2; %Flow angle at inducer, for alpha1 = pi/2 rad, assumed no pre-whirl
A = ((pi*D1^2)/4)*4;
Ai = ((pi*D1^2)/4); % Cross section area impeller(m^2)
Ad = ((pi*D1^2)/4); % Cross section area diffuser(m^2)
li = 0.053; %Mean channel length(m)
ld = 0.053; %Mean channel length(m)
beta1b = 0.61; %Impeller blade inlet angle()
Ch = 4*0.3164*(Re)^-0.25; %Friction loss coefficient()
kf = (Ch*li)/(2*Di*ro1^2*Ai^2*(sin(beta1b))^2)*4; %Friction loss impeller
kfd = (Ch*ld)/(2*Dd*ro1^2*Ad^2*(sin(beta1b))^2)*4; %Friction loss diffuser
sigma = 0.9; %Slip factor()
D2 = 0.128; %Diameter of impeller tip(m)
alfa2b = atan((D1*tan(beta1b))/(sigma*D2));
m = [-0.2:0.001:0.8]; %[kg/s] Mass Flow
cn=10;
%% Speed
giri_medium = 55000; %throttle valve aperture
delta_giri = 400; %range between your gamma_T_max e gamma_T_min
giri_max = giri_medium+(delta_giri/2); %maximum aperture
giri_min = giri_medium-(delta_giri/2); %minimum aperture
A_rpm = (giri_max - giri_min)/2; %amplitude
b_rpm = (giri_max + giri_min)/2;
rpm_frequency = 20; %engine's frequency [Hz]
w_rpm = rpm_frequency*2*pi; %engine's frequency [radiants/s]
% time_rpm=[25:0.01:40]; % original
time_rpm=[25:1:40]; % MN
giri = A_rpm*sin(w_rpm*time_rpm)+b_rpm;
wh = a01.*(Ac./Vp.*Lc).^0.5; %Helmotz's frequency
%% Efficency Losses
deltanbf = 0.02; %Back flow loss() 0.02 annular diffuser - 0,03 medium volume - 0.05 vaned diffuser
deltanv = 0.035; %Volute loss()
deltanc = 0.00; %Clearance loss()
deltand = 0.00; %Diffusion loss()
deltan = deltanbf+deltanv+deltanc+deltand;
%% Throttle valve's parameters
gamma_T = .17; %throttle valve aperture
Delta_gamma_T = .02; %range between your gamma_T_max e gamma_T_min
gamma_T_max = gamma_T+(Delta_gamma_T/2); %maximum aperture
gamma_T_min = gamma_T-(Delta_gamma_T/2); %minimum aperture
A = (gamma_T_max - gamma_T_min)/2; %amplitude
b = (gamma_T_max + gamma_T_min)/2;
%% Compressor's frequency (only frequency)
compressor_frequency = 10; %engine's frequency [Hz]
w = compressor_frequency*2*pi; %engine's frequency [radiants/s]
% save main_parameters.mat
%% Ode's parameters
t0(1) = 0; %simulation's start [s]
tf(1) = time_rpm(2); %simulation's finish [s]
y1(1) = 0; %initial condition 1
y2(1) = 0; %initial condition 2
init = [y1 y2]'; %initial conditions vector
sample_frequency = compressor_frequency*8; %sample frequency [Hz] (choose at least a frequency
% 10-12 times the one being analyzed
% for a good definition. Just 2-3 times for a quick try)
maximum_time_step = 1/sample_frequency;
options= odeset('MaxStep',maximum_time_step); %maximum time-step size
save main_parameters.mat
t_all = []; % MN
y_all = []; % MN
for i=1:(length(giri)-2)
U1(i) = (D1*pi.*giri(i))/60;
U2(i) = (D2*pi.*giri(i))/60;
B(i) = U1(i)./(2.*wh.*Lc);
[t,y]=ode113(@greitzer,[t0(i),tf(i)],[y1(i),y2(i)],options,U1(i),U2(i),B(i)); %I found ode113 is way more efficient
t0(i+1) = time_rpm(i+1); %simulation's start [s]
tf(i+1) = time_rpm(i+2); %simulation's finish [s]
y1(i+1) = y(end,1); %initial condition 1
y2(i+1) = y(end,2); %initial condition 2
% MN : concat t, y for all iterations
t_all = [t_all; t]; % MN
y_all = [y_all; y]; % MN
% array_val = zeros(length(giri)-2) ;
% plot(t,y(:,1))
% grid on
% grid minor
% xlabel('t [s]')
% ylabel('\Phi')
% hold on
% figure()
% plot(t,y(:,2))
% grid on
% grid minor
% xlabel('t [s]')
% ylabel('\Psi')
% hold on
% savefig('10Hz')
end
plot(t0,y1)
grid on
grid minor
xlabel('t [s]')
ylabel('\Phi')
% hold on
figure()
plot(t0,y2)
grid on
grid minor
xlabel('t [s]')
ylabel('\Psi')
% hold on
figure()
plot(t_all,y_all)
grid on
grid minor
xlabel('t [s]')
ylabel('y all')
Paul Rogers 2021 年 1 月 11 日
noope, post it a new answer, not as a comment to the first post, or i won't get the accept button

サインインしてコメントする。

### 回答 (1 件)

Ryan Fedeli 2021 年 1 月 8 日
This example is simplified, but to save a value from each loop iteration, all you need to do is save it into an array. The lines with asterisks are the important parts:
z = linspace(1,50,10); % arbitrary vector
b = z; % arbitrary values for example function
array_val = zeros(1,length(z)); % preallocate array for loop speed ***
for i = 1:length(z)
a = b(i)*1.1; % example function to create data
array_val(i) = a; % save the value from this iteration into an array ***
end
disp(array_val)
But if "a" in this example is not a single value, but an array itself, you'd need to save it into a matrix. In this case, it would look something like:
array_val = zeros(length(z)); % 1) preallocate *matrix* of zeros, length(z) by length(z) in size.
% 2) you'll need to figure out what these dimensions should be.
% 3) This line is helpful so the loop runs faster, but it's not necessary
for i = 1:length(z)
a = b(i).*z; % example function to create a row vector of data (array of data)
array_val(:,i) = a; % save the data in the row vector into the matrix in the i'th position
end
disp(array_val)
array_val is now a matrix, where each row is the data from the i'th iteration of the loop
I hope this helps
##### 1 件のコメント-1 件の古いコメントを表示-1 件の古いコメントを非表示
Paul Rogers 2021 年 1 月 8 日
I am trying to follow your heads up, but in your case you know that a is going too have the same length a z, infacct I get a Subscripted assignment dimension mismatch.
In my case, I don't know how many points the ode is gooing to use for that step.

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeSpectral Measurements についてさらに検索

R2014b

### Community Treasure Hunt

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

Start Hunting!

Translated by