How to write a for loop to generate a new set of initial conditions based on a input value that changes over different time intervals.

I have a ODE system with four equations and a input variable "Input".

I would like to show the output of species B, C and Input on seperate plots over time, similar to the plots attached in the screenshot (generated in Python).

The Input changes over time as follows:

Input = 0.5 for t < 50,

Input = 1 for t < 100,

Input = 1.5 for t < 150,

Input = 1 for t < 200,

Input = 0.5 for t < 250.

The Matlab code I have so far is:

% Set the initial values

A = 1;

B = 1;

C = 1;

D = 1;

Input = 1;

% Set the model parameters

k1 = 1;

k2 = 3;

k3 = 2;

k4 = 1;

k5 = 50;

k6 = 1;

k = [k1,k2,k3,k4,k5,k6];

init = [A B C D];

tspan = [0 200];

%t = linspace (0,200,100);

Looping Process

%Add code here..

% Perform the numerical integration

[t,u] = ode45(@(t,u) gene(t,u,k,Input), tspan, init);

plot(t,u(:,1),'--', t,u(:,2),'-', t,u(:,3),'--',t,u(:,4),'--','LineWidth',2.0)

title('');

xlabel('Time t');

ylabel('Solution y');legend('u_1','u_2','u_3','u_4')

ODE System

function eqns = gene(t,u,k,Input)

eqns = zeros(4,1); % To start with we havefour empty equations

% Using u = [A B C D]

% Using k = [k1,k2,k3,k4,k5,k6]

eqns(1) = k(1)*u(4) - k(2)*u(1);

eqns(2) = k(3)*Input*u(1) - k(4)*u(2);

eqns(3) = k(4)*u(2) - k(5)*u(3)*u(4);

eqns(4) = k(6) - k(5)*u(3)*u(4);

end

Thanks in advance.

Answer 1

[William Rose]

@Ron_S,

See attached code. I think you will be able to uunderstand what I have done by comparing the code to your version. It generates the plot below. I define Input inside function gene(). Good luck!

[Ron_S]

@William Rose,

Thank you so much! This is exactly the output I was after and the code is clear to understand. Merry Christmas!

[William Rose]

@Ron_S, you're welcome. You too!

Answer 2

[VBBV]

% Set the initial values
A = 1;
B = 1; 
C = 1;
D = 1;
% % Input = 1;
% Set the model parameters
k1 = 1;
k2 = 3;
k3 = 2;
k4 = 1;
k5 = 50;
k6 = 1;
k = [k1,k2,k3,k4,k5,k6];
init = [A B C D];
tspan = [0 250];
%t = linspace (0,200,100);
Input = [0.5 1 1.5 1 0.5];

 
    [t1,u1] = ode45(@(t,u) gene(t,u,k,Input(1)), tspan, init);
    idx1 = find(t1 < 50);
     
    
    [t2,u2] = ode45(@(t,u) gene(t,u,k,Input(2)), tspan, init);
    idx2 = find(t2 > 50 & t2 < 100 );
     
     
    [t3,u3] = ode45(@(t,u) gene(t,u,k,Input(3)), tspan, init);
    idx3 = find(t3 > 100 & t3 < 150 );
      
    [t4,u4] = ode45(@(t,u) gene(t,u,k,Input(4)), tspan, init);
    idx4 = find(t4 > 150 & t4 < 200 );
     
     
    [t5,u5] = ode45(@(t,u) gene(t,u,k,Input(5)), tspan, init);
    idx5 = find(t5 > 200 & t5 < 250 );
      
    
    t = [t1(idx1); t2(idx2); t3(idx3); t4(idx4); t5(idx5)];
    U1 = [u1(idx1,1); u2(idx2,1); u3(idx3,1); u4(idx4,1); u5(idx5,1)];
    U2 = [u1(idx1,2); u2(idx2,2); u3(idx3,2); u4(idx4,2); u5(idx5,2)];
    U3 = [u1(idx1,3); u2(idx2,3); u3(idx3,3); u4(idx4,3); u5(idx5,3)];
    U4 = [u1(idx1,4); u2(idx2,4); u3(idx3,4); u4(idx4,4); u5(idx5,4)];
    
    figure(1)
    plot(t,U1,t,U2,'--',t,U3,'-.',t,U4,'k-')
   
title('');
xlabel('Time t');
ylabel('Solution y');
legend('u_1','u_2','u_3','u_4')

% ODE System
function eqns = gene(t,u,k,Input)
eqns = zeros(4,1); % To start with we havefour empty equations
% Using u = [A B C D] 
% Using k = [k1,k2,k3,k4,k5,k6]
eqns(1) = k(1)*u(4) - k(2)*u(1);
eqns(2) = k(3)*Input*u(1) - k(4)*u(2);
eqns(3) = k(4)*u(2) - k(5)*u(3)*u(4);
eqns(4) = k(6) - k(5)*u(3)*u(4);
end

[Ron_S]

Thank you for providing this code VBBV, it's good to learn the different ways of getting to a solution!

[VBBV]

Yes, the previous code captures transients only at beginning of time, but does not capture transients for different inputs at later times, The below one does

% Set the initial values
A = 1;
B = 1; 
C = 1;
D = 1;
% % Input = 1;
% Set the model parameters
k1 = 1;
k2 = 3;
k3 = 2;
k4 = 1;
k5 = 50;
k6 = 1;
k = [k1,k2,k3,k4,k5,k6];
init = [A B C D];
tspan = [0 250];
%t = linspace (0,200,100);
Input = [0.5 1 1.5 1 0.5];

 
    [t1,u1] = ode45(@(t,u) gene(t,u,k,Input(1)),[0 50], init);
    idx1 = find(t1 <= 50);
     
    
    [t2,u2] = ode45(@(t,u) gene(t,u,k,Input(2)),[51 100], init);
    idx2 = find(t2 > 50 & t2 <= 100 );
     
     
    [t3,u3] = ode45(@(t,u) gene(t,u,k,Input(3)),[101 150], init);
    idx3 = find(t3 > 100 & t3 <= 150 );
      
    [t4,u4] = ode45(@(t,u) gene(t,u,k,Input(4)),[151 200], init);
    idx4 = find(t4 > 150 & t4 <= 200 );
     
     
    [t5,u5] = ode45(@(t,u) gene(t,u,k,Input(5)),[201 250], init);
    idx5 = find(t5 > 200 & t5 <= 250 );
      
    
    t = [t1(idx1); t2(idx2); t3(idx3); t4(idx4); t5(idx5)];
    U1 = [u1(idx1,1); u2(idx2,1); u3(idx3,1); u4(idx4,1); u5(idx5,1)];
    U2 = [u1(idx1,2); u2(idx2,2); u3(idx3,2); u4(idx4,2); u5(idx5,2)];
    U3 = [u1(idx1,3); u2(idx2,3); u3(idx3,3); u4(idx4,3); u5(idx5,3)];
    U4 = [u1(idx1,4); u2(idx2,4); u3(idx3,4); u4(idx4,4); u5(idx5,4)];
    
    figure(1)
    plot(t,U1,t,U2,'--',t,U3,'-.',t,U4,'m-','linewidth',2)
   
title('');
xlabel('Time t');
ylabel('Solution y');
legend('u_1','u_2','u_3','u_4','location','best')

% ODE System
function eqns = gene(t,u,k,Input)
eqns = zeros(4,1); % To start with we havefour empty equations
% Using u = [A B C D] 
% Using k = [k1,k2,k3,k4,k5,k6]
eqns(1) = k(1)*u(4) - k(2)*u(1);
eqns(2) = k(3)*Input*u(1) - k(4)*u(2);
eqns(3) = k(4)*u(2) - k(5)*u(3)*u(4);
eqns(4) = k(6) - k(5)*u(3)*u(4);
end

[Ron_S]

@VBBV Thank you so much! Is there a way to add a plot of Input over the different time intervals or perhaps as a subplot?

[William Rose]

@Ron_S,

There is no simple elegant way to plot the input versus time, using the code as it is currently structured.

There is an inelegant way: Copy the line below, which is in genes(), to the main script body:

Input = 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);

and, in the main script body, modify it as follows:

t1=0:.5:200;        %vector of times
Input1 = ones(1,length(t1))*0.5 + 0.5*(50<=t1 & t1<100) + 1*(100<=t1 & t1<150) + 0.5*(150<=t1 & t1<200);

Input1 will be a vector with same size at t1. Then you can modify the plot command:

plot(t,u(:,1),'-r', t,u(:,2),'-g', t,u(:,3),'-b',t,u(:,4),'-m',t1,Input1,'-k',LineWidth',2.0)
title(''); legend('u_1','u_2','u_3','u_4','Input')

Then you will see Input1 on the plot, as shown below:

The solution above is inelegant and suboptimal because, if you change Input in function genes(), it will not automatically be changed in the plot, and vice versa, so the user has to remember to do both, if they do either.

A more elegant solution is to define an anonymous function in the main body, and pass it to genes, and use it for plotting. This means there will only be one "Input", which is used by genes() and for plotting. To implement this solution, instead of doing the above, do the following:

Add line

input1 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);

to the main script, somewhere before the call to ode45(). This creates anonymous function input1(). Change ode45() to

[t,u] = ode45(@(t,u) gene(t,u,k,input1), tspan, init);

Change function genes() to

%% ODE System
function eqns = gene(t,u,k,in1)
% Using u = [A B C D] 
% Using k = [k1,k2,k3,k4,k5,k6]
eqns = zeros(4,1); % To start with we have four empty equations
eqns(1) = k(1)*u(4) - k(2)*u(1);
eqns(2) = k(3)*in1(t)*u(1) - k(4)*u(2);
eqns(3) = k(4)*u(2) - k(5)*u(3)*u(4);
eqns(4) = k(6) - k(5)*u(3)*u(4);
end

Change the plot command to

t1=0:.5:200;
plot(t,u(:,1),'-r', t,u(:,2),'-g', t,u(:,3),'-b',t,u(:,4),'-m',t1,input1(t1),'-k','LineWidth',2.0)
title(''); legend('u_1','u_2','u_3','u_4','Input')

This produces a plot identical to the plot above.

There is a built-in Matlab function input(), so I would rather not use "Input" as a variable name or function name, even though it is distinct (due to Matlab's case sensitivity).

See attached.

[VBBV]

% Set the initial values
A = 1;
B = 1; 
C = 1;
D = 1;
% % Input = 1;
% Set the model parameters
k1 = 1;
k2 = 3;
k3 = 2;
k4 = 1;
k5 = 50;
k6 = 1;
k = [k1,k2,k3,k4,k5,k6];
init = [A B C D];
tspan = [0 250];
%t = linspace (0,200,100);
Input = [0.5 1 1.5 1 0.5];


[t1,u1] = ode45(@(t,u) gene(t,u,k,Input(1)),[0 50], init);
idx1 = find(t1 <= 50);


[t2,u2] = ode45(@(t,u) gene(t,u,k,Input(2)),[51 100], init);
idx2 = find(t2 > 50 & t2 <= 100 );


[t3,u3] = ode45(@(t,u) gene(t,u,k,Input(3)),[101 150], init);
idx3 = find(t3 > 100 & t3 <= 150 );

[t4,u4] = ode45(@(t,u) gene(t,u,k,Input(4)),[151 200], init);
idx4 = find(t4 > 150 & t4 <= 200 );


[t5,u5] = ode45(@(t,u) gene(t,u,k,Input(5)),[201 250], init);
idx5 = find(t5 > 200 & t5 <= 250 );


t = [t1(idx1); t2(idx2); t3(idx3); t4(idx4); t5(idx5)];
U1 = [u1(idx1,1); u2(idx2,1); u3(idx3,1); u4(idx4,1); u5(idx5,1)];
U2 = [u1(idx1,2); u2(idx2,2); u3(idx3,2); u4(idx4,2); u5(idx5,2)];
U3 = [u1(idx1,3); u2(idx2,3); u3(idx3,3); u4(idx4,3); u5(idx5,3)];
U4 = [u1(idx1,4); u2(idx2,4); u3(idx3,4); u4(idx4,4); u5(idx5,4)];

figure(1)
subplot(411)
plot(t,U1,'linewidth',2); ylim([0 1.5]);title('U1'); 
subplot(412)
plot(t,U2,'linewidth',2);ylim([0 1.5]);title('U2'); 
subplot(413)
plot(t,U3,'linewidth',2);ylim([0 1.5]);title('U3'); 
subplot(414)
plot(t,U4,'linewidth',2);ylim([0 3.5]);title('U4');xlabel('time [s]')

    % response plot of signal U1 at different time intervals 
    figure('Name','u1')
    subplot(511)
    plot(t1,u1(idx1,1),'linewidth',2); ylim([0 1.5]); 
    subplot(512)
    plot([t1; t2],[u1(idx1,1); u2(idx2,1)],'linewidth',2);ylim([0 1.5]); 
    subplot(513)
    plot([t1 ;t2 ;t3],[u1(idx1,1); u2(idx2,1); u3(idx3,1)],'linewidth',2);ylim([0 1.5]); 
    subplot(514)
    plot([t1 ;t2 ;t3 ;t4],[u1(idx1,1); u2(idx2,1); u3(idx3,1); u4(idx4,1)],'linewidth',2);ylim([0 1.5]);   
    subplot(515)
    plot([t1; t2; t3 ;t4 ;t5],[u1(idx1,1); u2(idx2,1); u3(idx3,1); u4(idx4,1); u5(idx5,1)],'linewidth',2);ylim([0 1.5]); xlabel('time [s]')    

% xlabel('Time t');
% ylabel('Solution y');
% legend('u_1','u_2','u_3','u_4','location','best')
% ODE System
function eqns = gene(t,u,k,Input)
eqns = zeros(4,1); % To start with we havefour empty equations
% Using u = [A B C D] 
% Using k = [k1,k2,k3,k4,k5,k6]
eqns(1) = k(1)*u(4) - k(2)*u(1);
eqns(2) = k(3)*Input*u(1) - k(4)*u(2);
eqns(3) = k(4)*u(2) - k(5)*u(3)*u(4);
eqns(4) = k(6) - k(5)*u(3)*u(4);
end

[Ron_S]

@William Rose, thank you so much for your considered response and detailed explanations, this is the exact outcome I was looking for!

[Ron_S]

@VBBV Thank you again for providing alternative insights in relation to the question, much appreciated.

[William Rose]

@Ron_S, you're welcome. Good luck with your work.