How to solve an ODE calculating the value of a variable inside a loop?

Hello, I am having some issues including my dummy variable Xv inside the loop. I want to predict Xv at every point of Tspan.

My script:

%% Initialization

clear; close all; clc

%% User-defined decisions

% Initial conditions

y(1) = 3 ; %Bioreactor volume, L

y(2) = 1000 ; %Dummy variable to include in loop

x0 = [y] ;

%Feeding Profile/Time spam

tspan = [0 10 10.01 20 20.01 30 30.1 40 ] ; %Sampling/Feeding Events

F_glc = [0 2 0 5 0 1 0 0 ] ; %Glucose feed (L/h)

F_out = [0 0.04 0 0.04 0 0.04 0 0 ] ; %Sampling (L/h)

system.F_glc =F_glc' ;

system.F_out =F_out' ;

% Define initial empty arrays for storing solutions and time points

Xs = [];

Ts = [];

%Loop at each sampling/feeding point

for i = 1:(length(tspan)-1)

% Define parameters for current timespan

system.F_glc = F_glc(i);

system.F_out = F_out(i);

% Compute solutions

[T,X] = ode45(@(t,x) bioreactor(t,x,system), [tspan(i) tspan(i+1)], x0);

% Set initial point for next iteration

x0 = X(end);

% Save solutions

Ts = [Ts; T];

Xs = [Xs; X];

end

% Create figure

figure

plot(Ts,Xs(:,1),'.')

------

%my function

function dy=bioreactor(t,y,system)

%% State Variables

% Unwrap variables

F_glc = system.F_glc;

F_out = system.F_out;

V = y(1);

Xv = y(2);

%Volume of Bio

dVdt = F_glc - F_out;

%Dummny variable Test

dXvdt = Xv*F_glc ;

dy =[dVdt; dXvdt] ;

end

I am getting the folliwing error:Index exceeds the number of array elements (1)

Error in Testfeeding>bioreactor (line 52)

Xv = y(2);

I want to calculate Xv at the same points of dVdt

Answer 1

[Star Strider]

The problem is:

x0 = X(end);

so the error only occurs after the first iteration.

The solution is:

x0 = X(end,:)

and thin it works!

% Initial conditions
y(1)         =   3                                                ;        %Bioreactor volume, L
y(2)         =   1000                                             ;        %Dummy variable to include in loop
x0           =   [y]                                                       ;        
%Feeding Profile/Time spam
tspan   =  [0    10   10.01   20    20.01  30     30.1  40 ] ;        %Sampling/Feeding Events 
F_glc   =  [0     2    0      5     0      1      0     0  ] ;        %Glucose feed (L/h)
F_out   =  [0   0.04   0      0.04  0      0.04   0     0  ] ;        %Sampling     (L/h)
system.F_glc  =F_glc'                                          ;
system.F_out  =F_out'                                          ;
% Define initial empty arrays for storing solutions and time points
Xs = [];
Ts = [];
%Loop at each sampling/feeding point
for i = 1:(length(tspan)-1)
    % Define parameters for current timespan
    system.F_glc = F_glc(i);
    system.F_out = F_out(i);
    % Compute solutions
    [T,X] = ode45(@(t,x) bioreactor(t,x,system), [tspan(i) tspan(i+1)], x0);
    
    % Set initial point for next iteration
    x0 = X(end,:);
    
    % Save solutions 
    Ts = [Ts; T];
    Xs = [Xs; X];
end
% Create figure
figure
plot(Ts,Xs(:,1),'.')

% ------
%my function
function dy=bioreactor(t,y,system)
%% State Variables
% Unwrap variables
F_glc = system.F_glc;
F_out = system.F_out;
V     = y(1);
Xv    = y(2);
%Volume of Bio
dVdt  = F_glc - F_out; 
%Dummny variable Test
dXvdt = Xv*F_glc     ;
dy    =[dVdt; dXvdt] ;
end

.

[Rodrigo Barbosa]

Thanks a lot! It solved. But once I try to add components and parameters it gets stuck on "Not enough input arguments."

Error in Testfeeding>profile (line 81)

Y_SX= Param(1);

Would you please kindly help me on this error? thank you!

%% Initialization

clear; close all; clc

%% User-defined decisions

% Initial conditions

y(1) = 3 ; %Bioreactor volume, L

y(2) = 153.6 ;

y(3) = 1.761 ;

y(4) = 0 ;

x0 = [y] ;

%

%Parameters

Y_SX = 0.84;

Y_PX = 0.0064;

k_DW= 3.3;

mu_max = 0.1466;

K= 1.5;

Ki = 8.86;

%

%Parameters

Param(1)= 0.84;

Param(2)= 0.0064;

Param(3)= 3.3;

Param(4)= 0.1466;

Param(5)= 1.5;

Param(6)= 8.86;

p=[Param];

%

%Feeding Profile/Time spam

tspan = [0 10 10.01 20 20.01 30 30.1 40 ] ; %Sampling/Feeding Events

F_glc = [0 2 0 5 0 1 0 0 ] ; %Glucose feed (L/h)

F_out = [0 0.04 0 0.04 0 0.04 0 0 ] ; %Sampling (L/h)

system.F_glc =F_glc' ;

system.F_out =F_out' ;

% Define initial empty arrays for storing solutions and time points

Xs = [];

Ts = [];

%Loop at each sampling/feeding point

for i = 1:(length(tspan)-1)

% Define parameters for current timespan

system.F_glc = F_glc(i);

system.F_out = F_out(i);

% Compute solutions

[T,X] = ode45(@(t,x) profile(t,x,system), [tspan(i) tspan(i+1)], x0);

% Set initial point for next iteration

x0 = X(end,:);

% Save solutions

Ts = [Ts; T];

Xs = [Xs; X];

end

% Create figure

figure

plot(Ts,Xs(:,1),'.')

% function dy=profile(t,y,system,Param)

function dy=profile(t,y,system, Param)

% Unwrap variables

F_glc = system.F_glc;

F_out = system.F_out;

%State variables

V = y(1);

Xv = y(2);

S = y(3);

P =y(4);

%

% %Parameters

% %

Y_SX= Param(1);

Y_PX= Param(2);

k_DW= Param(3);

mu_max= Param(4);

K= Param(5);

Ki= Param(6);

%% Algebraic Equations

%

% mu=mu_max*S/(K+S+(S^2/Ki));

% S_in = 10;

%% ODE

dVdt= F_glc - F_out ;

% dXvdt= -F_glc/V*Xv + mu*Xv ;

dXvdt= -F_glc/V*Xv ;

dSdt= F_glc/V*(S_in-S) - Y_SX*mu*Xv ;

dPdt= - F_glc/V*P + k_DW*Y_PX*mu*Xv ;

%

dy =[dVdt; dXvdt;dSdt;dPdt ] ;

end

[Star Strider]

My pleasure!

Add ‘Param’ to the function call —

[T,X] = ode45(@(t,x) profile(t,x,system,Param), [tspan(i) tspan(i+1)], x0);

With that (and a few other) corrections, it now runs without error and gives what I consider to be appropriate results —

% Initial conditions
y(1)         =   3                                                ;        %Bioreactor volume, L
y(2)         = 153.6                                              ;
y(3)         = 1.761                                              ;
y(4)         = 0                                                  ;
x0           =   [y]                                              ;        
% 
%Parameters
Y_SX = 0.84;
Y_PX = 0.0064;
k_DW=  3.3;
mu_max = 0.1466;
K= 1.5;
Ki = 8.86;
% 
%Parameters
Param(1)= 0.84;
Param(2)= 0.0064;
Param(3)=  3.3;
Param(4)= 0.1466;
Param(5)= 1.5;
Param(6)= 8.86;
p=[Param];
% 
%Feeding Profile/Time spam
tspan   =  [0    10   10.01   20    20.01  30     30.1  40 ] ;        %Sampling/Feeding Events 
F_glc   =  [0     2    0      5     0      1      0     0  ] ;        %Glucose feed (L/h)
F_out   =  [0   0.04   0      0.04  0      0.04   0     0  ] ;        %Sampling     (L/h)
system.F_glc  =F_glc'                                        ;
system.F_out  =F_out'                                        ;
% Define initial empty arrays for storing solutions and time points
Xs = [];
Ts = [];
%Loop at each sampling/feeding point
for i = 1:(length(tspan)-1)
    % Define parameters for current timespan
    system.F_glc = F_glc(i);
    system.F_out = F_out(i);
    % Compute solutions
    [T,X] = ode45(@(t,x) profile(t,x,system,Param), [tspan(i) tspan(i+1)], x0);
    
    % Set initial point for next iteration
    x0 = X(end,:);
    
    % Save solutions 
    Ts = [Ts; T];
    Xs = [Xs; X];
end
% Create figure
figure
plot(Ts,Xs(:,1),'.')

% function dy=profile(t,y,system,Param)
function dy=profile(t,y,system, Param)
% Unwrap variables
F_glc = system.F_glc;
F_out = system.F_out;
%State variables
V     = y(1); 
Xv    = y(2);
S     = y(3);
P      =y(4);
% 
% %Parameters
% % 
Y_SX=   Param(1);
Y_PX=   Param(2);
k_DW=   Param(3);
mu_max= Param(4);
K=      Param(5);
Ki=     Param(6);
%% Algebraic Equations
% 
 mu=mu_max*S/(K+S+(S^2/Ki));
 S_in    = 10;
%% ODE
dVdt=   F_glc - F_out                 ;
%  dXvdt= -F_glc/V*Xv + mu*Xv            ;
dXvdt= -F_glc/V*Xv          ;
dSdt=   F_glc/V*(S_in-S) - Y_SX*mu*Xv ;
dPdt= - F_glc/V*P + k_DW*Y_PX*mu*Xv   ;
% 
dy    =[dVdt; dXvdt;dSdt;dPdt ]        ;

end

.

[Star Strider]

If my Answer helped you solve your problem, please Accept it!

.

[Rodrigo Barbosa]

You are an star!!

[Star Strider]

Thank you!

.