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.

Question

Ron_S 2022-12-22

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1882917-how-to-write-a-for-loop-to-generate-a-new-set-of-initial-conditions-based-on-a-input-value-that-chan

评论： William Rose 2022-12-26

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

William Rose 2022-12-22

kineticModel1.m

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

2 个评论
显示无隐藏无

Ron_S 2022-12-22

@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 2022-12-22

@Ron_S, you're welcome. You too!

请先登录，再进行评论。

Answer 2

VBBV 2022-12-22

编辑：VBBV 2022-12-22

在 MATLAB Online 中打开

% 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

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

VBBV 2022-12-22

编辑：VBBV 2022-12-22

在 MATLAB Online 中打开

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

William Rose 2022-12-23

在 MATLAB Online 中打开

kineticModel1.m

@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 2022-12-23

在 MATLAB Online 中打开

% 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 2022-12-26

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

William Rose 2022-12-26

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

请先登录，再进行评论。

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

2 个评论
显示无隐藏无

更多回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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.

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

2 个评论 显示 无隐藏 无

更多回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

2 个评论
显示无隐藏无

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