フィルターのクリア

How do I store an ODE in a vector?

3 ビュー (過去 30 日間)
Chiara
Chiara 2023 年 3 月 21 日
編集済み: Torsten 2023 年 3 月 21 日
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

回答 (1 件)

Torsten
Torsten 2023 年 3 月 21 日
編集済み: Torsten 2023 年 3 月 21 日
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

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

製品


リリース

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by