storing variables after a loop (easy)

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

[Mathieu NOE]

hello

I see you have also started a similar post there

ode in a for loop updating initial conditions - MATLAB Answers - MATLAB Central (mathworks.com)

now looking at this specific question, what do you mean by " I only get the last part of the loop" ?

I understand that you are solving the ODE in the loop, but for each time the start an end times of the ODE are different, so I do not understand how you can have only identical t vector and y vector dimensions for each simulation;

do you intend to overlay the plots or stack / concatenate the results ?

could you provide the entire code / data to test it ?

[Cris LaPierre]

Conveniently enough, there is an example of how to use indexing to save values from each loop on the for loop documentation page.

You might also consider going through a couple of the relevant chapters in MATLAB Onramp as well (Ch5 and Ch 13 specifically).

[Paul Rogers]

@Mathieu NOE

here's the files. I tryied as i was told in the latter post, but nothing came out

[Mathieu NOE]

hello Paul

I had started to work on your files from the other post but I got some errors

I see you share a new version of the code

First, to speed up a bit the process , I have just changed these two lines , but you can easily come back to your own tuning :

% time_rpm=[25:0.01:40]; % original
time_rpm=[25:1:40]; % MN

....
    
sample_frequency = compressor_frequency*8; %sample frequency [Hz] (choose at least a frequency  

then, for what I understand, you would like to concatenate the results of all iterations of the ode solver - and that's what I did (nothing fancy) and plotted in the last plot (look for the lines with % MN)

here is the entire code :

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]

thanks a lot,

% MN : concat t, y for all iterations
t_all = [t_all; t]; % MN
y_all = [y_all; y]; % MN

this is the line I was looking foor and didn't know about.

For the other tips I don't agree.

% time_rpm=[25:0.01:40]; % original
time_rpm=[25:1:40]; % MN

when I have an higher

rpm_frequency = 20

I need more points in the vector

time_rpm=[25:0.01:40];

Anyway, thanks a lot on how to get all the datas from the iterations

[Paul Rogers]

and I am sorry I can't accept the answer since it's a comment to my first post.

If you repost the answer below I will accept.

[Mathieu NOE]

sure

Glad this little mods solved your problem

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]

noope, post it a new answer, not as a comment to the first post, or i won't get the accept button

Answer 1

[Ryan Fedeli]

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

[Paul Rogers]

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.