Filter löschen
Filter löschen

Why won't Matlab plot or finish running my code due to the value of 1 variable?

1 Ansicht (letzte 30 Tage)
I have a value known as h0 that I set as 10^-2 but whenever I run this with this value thge code never runs it just doesn'ty produce any values nor graphs. I am inputting this value as it represents mutants in this case y values that have a strong control over the population. Why wont this run?
this is what the code is based off of with tmp 1 being x0 and so on.
Code:
HematopoieticSystemModel
function [t,dydt1] = HematopoieticSystemModel
time = 0:3000;
%------%
x_0 = 1;
y_0 = 0;
x_1 = 0;
y_1 = 0;
x_2 = 0;
y_2 = 0;
% p0_0=100;
%----------%
p = [x_0 y_0 x_1 y_1 x_2 y_2];
%--------------------------------------------------%
[t,dydt1] = ode45(@myODE,time,p);
%---Matrices-----%
% [m1, m2, m3, m4, m5, m6] = equillibriatest(dydt1(:,1), dydt1(:,2), dydt1(:,3), dydt1(:,4), dydt1(:,5), dydt1(:,6));
% save('m1.mat','m1')
% save('m2.mat','m2')
% save('m3.mat','m3')
% save('m4.mat','m4')
% save('m5.mat','m5')
% save('m6.mat','m6')
%% -----Individual Plots------ %
% figure (1)
% hold on
% subplot (1,3,1);
% plot(t,dydt1(:, 1), 'r');
% title ('HSC');
% xlabel('time')
% ylabel('cells')
% hold off
%
% hold on
% subplot (1,3,2);
% plot (t, dydt1(:,3),'g');
% title ('ST-HSCs');
% xlabel('time')
% ylabel('cells')
% hold off
%
% hold on
% subplot (1, 3, 3);
% plot (t, dydt1(:,5),'b');
% title ('MPPs');
% xlabel('time')
% ylabel('cells')
% hold off
%------------------%
figure (2)
hold on
plot(t,dydt1(:, 1), 'g');
plot (t, dydt1(:,3),'b');
plot (t, dydt1(:,5),'m');
legend ('HSC' , 'ST-HSC', 'MPP');
xlabel('time (days)');
ylabel('cells');
title ('Hematopoietic System Growth');
set(gca, 'YScale', 'log');
hold off
%------------------%
figure (3)
hold on
plot(t,dydt1(:, 2), 'g');
plot (t, dydt1(:,4),'b');
plot (t, dydt1(:,6),'m');
legend ('HSC' , 'ST-HSC', 'MPP');
xlabel('time (days)');
ylabel('mutant cells');
title ('Hematopoietic System: Mutant Growth');
set(gca, 'YScale', 'log');
hold off
figure(4)
hold on
plot (t, dydt1(:,1),'g');
plot (t, dydt1(:,2), 'k');
title ('HSCs');
xlabel('time');
ylabel('cells');
legend ('HSC', 'mutant');
set(gca, 'YScale', 'log');
hold off
figure (5)
hold on
plot (t, dydt1(:,3),'b')
plot (t, dydt1(:,4),'k')
title ('ST-HSCs');
xlabel('time');
ylabel('cells');
legend ('ST-HSC', 'mutant');
set(gca, 'YScale', 'log');
hold off
figure (6)
hold on
plot (t, dydt1(:,5), 'm')
plot (t, dydt1(:,6),'k')
title ('MPPs');
xlabel('time');
ylabel('cells');
legend ('MPPs', 'mutants');
set(gca, 'YScale', 'log');
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);
global P0 P1 P2 h0 r0 r1 r2 d3 u1 s x0 u0 u2 h2 h1
d3=0.04;
s=0.03;
x0=17000;
r0=0.01818;
r1=0.8712;
r2=8.0217;
P0=0.5;
P1=0.4783;
P2=0.4987;
u0=10^-5;
u1=10^-5;
u2=10^-5;
%%strong control
h0= 10^-2;
h1= 10^-2;
h2= 10^-2;
% y(1) = x_0;
% y(2) = y_0;
% y(3) = x_1;
% y(4) = y_1;
% y(5) = x_2;
% y(6) = y_2;
c0=0.5+((h0*x0)/2);
c1=P1+((r0*x0*h1*P1)/(r1*(1-2*P1)));
c2=P2+((2*r0*x0*h2*P2*(1-P1))/(r2*(1-2*P1)*(1-2*P2)));
p0=c0./(1+h0.*(y(1)+y(2)));
p1=c1./(1+h1.*(y(3)+y(4)));
p2=c2./(1+h2.*(y(5)+y(6)));
p_0=P0.*(1+s);
p_1=P1.*(1+s);
p_2=P2.*(1+s);
tmp1 = r0.*y(1).*p0.*(1-u0)-r0.*y(1).*(1-p0);
tmp2 = r0.*y(1).*p0.*u0-r0.*y(2).*(2.*p_0-1);
%%
tmp3 = r0.*y(1).*(1-p0).*(2-u0)+r1.*y(3).*p1.*(1-u1)-r1.*y(3).*(1-p1);
tmp4 = r0.*y(1).*(1-p0).*u0+2.*r0.*y(2).*(1-p_0)+r1.*y(3).*p1.*u1+r1.*y(4).*(2*p_1-1);
%%
tmp5 =r1.*y(3).*(1-p1).*(2-u1)+r2.*y(5).*p2.*(1-u2)-r2.*y(5).*(1-p2);
tmp6 =r1.*y(3).*(1-p1).*u1+2.*r1.*y(4).*(1-p_1)+r2.*y(5).*p2.*u2+r2.*y(6).*(2*p_2-1);
dydt = [tmp1; tmp2; tmp3; tmp4; tmp5; tmp6];
end
function [matrix1, matrix2, matrix3, matrix4, matrix5, matrix6] = equillibriatest (in1, in2, in3, in4, in5, in6)
matrix1 = in1;
matrix2 = in2;
matrix3 = in3;
matrix4 = in4;
matrix5 = in5;
matrix6 = in6;
end

Antworten (1)

Walter Roberson
Walter Roberson am 12 Jun. 2023
You have to reduce the upper bound on the time by quite a bit. With you current code, with ode45, after running for about 45 minutes, it would not have proceeded much past 42 simulated seconds.
If you switch to ode15s, then it quits the integration after a small fraction of a second.
If you switch to ode23s, then after roughly 0.097 seconds it really starts to slow down; 0.098 as an upper bound is tolerable but beyond that it gets very slow.
With ode45 an upper bound of 20 seconds is tolerable time.
With both ode23s and ode45, dydt1(:,3) and dydt1(:,4) go negative early on, and so are mostly not plotted because you set log axes.
You also have a little bit of a settings problem. You hold on before drawing anything at all in the plots, and that locks in a ylim of [0 1]. You should use code closer to
figure (2)
plot(t,dydt1(:, 1), 'g');
hold on
plot (t, dydt1(:,3),'b');
plot (t, dydt1(:,5),'m');
legend ('HSC' , 'ST-HSC', 'MPP');
xlabel('time (days)');
ylabel('cells');
title ('Hematopoietic System Growth');
set(gca, 'YScale', 'log');
hold off
ylim auto
  1 Kommentar
Walter Roberson
Walter Roberson am 12 Jun. 2023
I ran with ode45 up to time t = 20. It took a while. I used [0, 20] rather than 0:20 for the time, so it emitted an output (four outputs actually) each time it took a successful step. It emitted over 12 million output rows up to t = 20.
dydt1(:, 3) and dydt1(:, 4) both go negative pretty much immediately, and they are both quite unstable -- nearly all of the hard work is in correctly tracing :3 and :4 . The :3 is pretty much pure noise; the :4 is noise overlayed on a distinct arc.
I would suggest that you re-check that your implementation matches your equations.

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu MATLAB finden Sie in Help Center und File Exchange

Produkte


Version

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by