# How do I store an ODE in a vector?

2 Ansichten (letzte 30 Tage)
Chiara am 21 Mär. 2023
Bearbeitet: Torsten am 21 Mär. 2023
I have a fully coded ODE modelling population growth of cells and I want to find the homeostasis levels or when the graph flattens out. in order to do this I need to store the ODEs, (tmp1, tmp2, tmp3, tmp4) which are under dydt, in vectors. How can I store each ODE into a vector and if this is confusing how is an easier way to find the point when the curve flattens and there is no change (slope=0).
Also note I cannot use the symbolic set as my license for matlab doesn't allow it.
function [t,dydt1] = HematopoieticSystemModel
time = 0:4000;
%------%
x_0 = 1;
x_1 = 3;
x_2 = 9;
x_3 = 50;
p0_0=100;
%----------%
p = [x_0 x_1 x_2 x_3 p0_0];
%--------------------------------------------------%
[t,dydt1] = ode45(@myODE,time,p);
%% -----Individual Plots------ %
figure (1)
hold on
subplot (2,2,1);
plot(t,dydt1(:, 1), 'r');
title ('HSC');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2,2,2);
plot (t, dydt1(:,2),'g');
title ('ST-HSCs');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2, 2, 3);
plot (t, dydt1(:,3),'b');
title ('MPPs');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2, 2, 4);
plot (t, dydt1(:,4), 'm');
title ('CLPs and CMPs');
xlabel('time')
ylabel('cells')
hold off
%------------------%
figure (2)
hold on
plot(t,dydt1(:, 1), 'r');
plot (t, dydt1(:,2),'g');
plot (t, dydt1(:,3),'b');
plot (t, dydt1(:,4), 'm');
legend ('HSC' , 'ST-HSC', 'MPP', 'CLP and CMP');
xlabel('time (days)');
ylabel('cells');
title ('Hematopoietic System Growth');
set(gca, 'YScale', 'log');
hold off
%------------------%
figure (3)
hold on
plot(t,dydt1(:, 1), 'r');
title ('HSC');
xlabel('time')
ylabel('cells')
hold off
figure(4)
hold on
plot (t, dydt1(:,2),'g');
title ('ST-HSCs');
xlabel('time')
ylabel('cells')
hold off
figure (5)
hold on
plot (t, dydt1(:,3),'b')
title ('MPPs');
xlabel('time')
ylabel('cells')
hold off
figure (6)
hold on
plot (t, dydt1(:,4), 'm')
title ('CLPs and CMPs');
xlabel('time')
ylabel('cells')
hold off
% figure(7)
% hold on
% plot (t, dydt1(:,5), 'k')
% title ('Self Renewal for HSCs');
% xlabel('time')
% ylabel('cells')
% hold off
end
function dydt = myODE(t,y, r);
P0=0.7;
h0=0.0000235;
r0=0.01818;
r1=0.08712;
r2=8.0217;
d3=0.1;
p1=0.4783;
p2=0.4987;
%%
% y(1) = x_0;
% y(2) = x_1;
% y(3) = x_2;
% y(4) = x_3;
p0=P0./(1+h0.*y(1));
tmp1 = r0.*y(1).*(2.*p0-1);
tmp2 = 2.*r0.*y(1).*(1-p0) + r1.*y(2).*(2.*p1-1);
tmp3 = 2.*r1.*y(2).*(1-p1) + r2.*y(3).*(2.*p2-1);
tmp4 =2.*r2.*y(3).*(1-p2) - d3.*y(4);
dydt = [tmp1; tmp2; tmp3; tmp4; p0];
end
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Antworten (1)

Torsten am 21 Mär. 2023
Bearbeitet: Torsten am 21 Mär. 2023
Either you use the code below and work with the time derivatives after the integration to decide when the curves have flattened or you use the Event facility of the ODE integrators to stop integration when the time derivatives of the equations become small enough for your purpose.
%function [t,dydt1] = HematopoieticSystemModel
time = 0:4000;
%------%
x_0 = 1;
x_1 = 3;
x_2 = 9;
x_3 = 50;
p0_0=100;
%----------%
p = [x_0 x_1 x_2 x_3 p0_0];
%--------------------------------------------------%
[t,y] = ode15s(@myODE,time,p);
dydt = zeros(size(y));
for i = 1:numel(t)
dydt(i,:) = myODE(t(i),y(i,:));
end
plot(t,dydt)
%% -----Individual Plots------ %
figure (1)
hold on
subplot (2,2,1);
plot(t,y(:, 1), 'r');
title ('HSC');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2,2,2);
plot (t, y(:,2),'g');
title ('ST-HSCs');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2, 2, 3);
plot (t, y(:,3),'b');
title ('MPPs');
xlabel('time')
ylabel('cells')
hold off
hold on
subplot (2, 2, 4);
plot (t, y(:,4), 'm');
title ('CLPs and CMPs');
xlabel('time')
ylabel('cells')
hold off
%------------------%
figure (2)
hold on
plot(t,y(:, 1), 'r');
plot (t, y(:,2),'g');
plot (t, y(:,3),'b');
plot (t, y(:,4), 'm');
legend ('HSC' , 'ST-HSC', 'MPP', 'CLP and CMP');
xlabel('time (days)');
ylabel('cells');
title ('Hematopoietic System Growth');
set(gca, 'YScale', 'log');
hold off
%------------------%
figure (3)
hold on
plot(t,y(:, 1), 'r');
title ('HSC');
xlabel('time')
ylabel('cells')
hold off
figure(4)
hold on
plot (t, y(:,2),'g');
title ('ST-HSCs');
xlabel('time')
ylabel('cells')
hold off
figure (5)
hold on
plot (t, y(:,3),'b')
title ('MPPs');
xlabel('time')
ylabel('cells')
hold off
figure (6)
hold on
plot (t, y(:,4), 'm')
title ('CLPs and CMPs');
xlabel('time')
ylabel('cells')
hold off
% figure(7)
% hold on
% plot (t, y(:,5), 'k')
% title ('Self Renewal for HSCs');
% xlabel('time')
% ylabel('cells')
% hold off
%end
function dydt = myODE(t,y)
P0=0.7;
h0=0.0000235;
r0=0.01818;
r1=0.08712;
r2=8.0217;
d3=0.1;
p1=0.4783;
p2=0.4987;
%%
% y(1) = x_0;
% y(2) = x_1;
% y(3) = x_2;
% y(4) = x_3;
p0=P0./(1+h0.*y(1));
tmp1 = r0.*y(1).*(2.*p0-1);
tmp2 = 2.*r0.*y(1).*(1-p0) + r1.*y(2).*(2.*p1-1);
tmp3 = 2.*r1.*y(2).*(1-p1) + r2.*y(3).*(2.*p2-1);
tmp4 =2.*r2.*y(3).*(1-p2) - d3.*y(4);
dydt = [tmp1; tmp2; tmp3; tmp4; p0];
end
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!

Translated by