Why my plot is wrong if I decrease dt

1 次查看（过去 30 天）
cindyawati cindyawati
评论： cindyawati cindyawati 2023-11-20
Ran in:
I have a code, while I running the code is nothing error in the result. But if dt decrease (ex: 0.001) the plot is wrong.
I know the plot is wrong because I do that in excel. In excel is nothing wrong. I don't know what's wrong with my plot code
clc;clear;
%parameter
delta=50;
gamma=75;
K1=10^-4;
K2=5*10^-4;
K3=10^-3;
K4=5*10^-3;
K5=10^-2;
K6=5*10^-2;
Ko=0.1;
n=6;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_4=10^-3;
mu_5=10^-3;
mu_6=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
%input for time
t(1)=0;
dt=0.001; %time interval
t=0:dt:180; %time span
%component for microglia
Mg = 0.047;
lambdaMf = 5*10^-2;
Fo = 3.36*10^-11;
KFo = 2.58*10^-11;
lambdaMa = 2.23*10^-2;
Ao = 0.14;
KAo = 10^-7;
beta = 10;
epsilon1 = 0.3333;
epsilon2 = 0.8;
lambdaM1Tb = 6*10^-3;
Tb=10^-6;
KTb = 2.5*10^-7;
makro1 = 0.02;
makro2 = 0.02;
dM1= 0.015;
dM2= 0.015;
tmakro1 = Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon1*beta/beta*(epsilon1+epsilon2))-(lambdaM1Tb*(Tb/(Tb+KTb)*makro1))-(dM1*makro1);
tmakro2 = Mg*(lambdaMf*(Fo/(Fo+KFo))+lambdaMa*(Ao/(Ao+KAo)))*(epsilon2/beta*(epsilon1+epsilon2))+(lambdaM1Tb*(Tb/(Tb+KTb))*makro1)-(dM2*makro2);
%component drugs
Abi = 10^-6;
lambdaN = 8*10^-4;
lambdaA = 8*10^-4;
daboM = 2*10^-3;
dabom = 10^-2; %clearance rate by macrophages
makrofag1 = 0;
makrofag2 = 0;
mikro1 = 0.02;
mikro2 = 0.02;
teta = 0.9;
h = 10;
Kaob = 7*10^-3;
AoB = 10^-8;
dtdab = Abi*(dt*(0.14*1.9*10^-4*exp(-1.9*10^-4)))+ (lambdaN*1)+(lambdaA*1)-(dabom*(makrofag1+(teta*makrofag2))+daboM*(mikro1+(teta*mikro2))*(1+h))*(AoB/(AoB+Kaob));
%initial condition with microglia
M1(1)=10;
M2(1)=0;
M3(1)=0;
M4(1)=0;
M5(1)=0;
M6(1)=0;
O(1)=0;
P(1)=0;
%initial condition without microglia
M12(1)=10;
M22(1)=0;
M32(1)=0;
M42(1)=0;
M52(1)=0;
M62(1)=0;
O2(1)=0;
P2(1)=0;
%initial condition with drugs depend time
M14(1)=10;
M24(1)=0;
M34(1)=0;
M44(1)=0;
M54(1)=0;
M64(1)=0;
O4(1)=0;
P4(1)=0;
%empty array with microglia
T=zeros(length(t)+1,1); %empty array for t
M1=zeros(length(t)+1,1); %empty array for M1
M2=zeros(length(t)+1,1); %empty array for M2
M3=zeros(length(t)+1,1); %empty array for M3
M4=zeros(length(t)+1,1); %empty array for M4
M5=zeros(length(t)+1,1); %empty array for M5
M6=zeros(length(t)+1,1); %empty array for M6
O=zeros(length(t)+1,1); %empty array for O
P=zeros(length(t)+1,1); %empty array for P
sumter=K2*M2+K3*M3+K4*M4+K5*M5;
%empty array without microglia
M12=zeros(length(t)+1,1); %empty array for M1
M22=zeros(length(t)+1,1); %empty array for M2
M32=zeros(length(t)+1,1); %empty array for M3
M42=zeros(length(t)+1,1); %empty array for M4
M52=zeros(length(t)+1,1); %empty array for M5
M62=zeros(length(t)+1,1); %empty array for M6
O2=zeros(length(t)+1,1); %empty array for O
P2=zeros(length(t)+1,1); %empty array for P
sumter2=K2*M22+K3*M32+K4*M42+K5*M52;
%empty array with drugs depend time
M14=zeros(length(t)+1,1); %empty array for M1
M24=zeros(length(t)+1,1); %empty array for M2
M34=zeros(length(t)+1,1); %empty array for M3
M44=zeros(length(t)+1,1); %empty array for M4
M54=zeros(length(t)+1,1); %empty array for M5
M64=zeros(length(t)+1,1); %empty array for M6
O4=zeros(length(t)+1,1); %empty array for O
P4=zeros(length(t)+1,1); %empty array for P
sumter2=K2*M24+K4*M34+K4*M44+K5*M54;
for j = 1:length(t)
%with microglia
T(j+1)=T(j)+dt;
M1(j+1)=M1(j)+1./(1+exp(-T(j)));
M1(j+1) = M1(j)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*sumter(j+1))-((Oa-n)*K6*M1(j+1)*M6(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1)));
M2(j+1) = M2(j)+(dt*(K1*M1(j)*M1(j)-K2*M1(j)*M2(j))-(mu_2*M2(j+1))-tmakro1-tmakro2);
M3(j+1) = M3(j)+(dt*(K2*M1(j)*M2(j)-K3*M1(j)*M3(j))-mu_3*M3(j));
M4(j+1) = M4(j)+(dt*(K3*M1(j)*M3(j)-K4*M1(j)*M4(j))-mu_4*M4(j));
M5(j+1) = M5(j)+(dt*(K4*M1(j)*M4(j)-K5*M1(j)*M5(j))-mu_5*M5(j));
M6(j+1) = M6(j)+(dt*(K5*M1(j)*M5(j)-K6*M1(j)*M6(j))-mu_6*M6(j));
O(j+1) = O(j)+(dt*(K6*M1(j)*M6(j)-Ko*M1(j)*O(j)-mu_o*O(j)));
P(j+1) = P(j)+(dt*(Ko*M1(j)*O(j)-mu_p*P(j)));
%without microglia
T(j+1)=T(j)+dt;
M12(j+1) = M12(j)+1./(1+exp(-T(j)));
M12(j+1) = M12(j)+(dt*(delta*M12(j+1)*(1-(M12(j+1)/gamma))-2*K1*M12(j+1)*M12(j+1)-M12(j+1)*sumter2(j+1))-((Oa-n)*K6*M12(j+1)*M62(j+1))-((Pa-Oa)*Ko*M12(j+1)*O2(j+1))-(mu_1*M12(j+1)));
M22(j+1) = M22(j)+(dt*(K1*M12(j)*M12(j)-K2*M12(j)*M22(j))-(mu_2*M22(j+1)));
M32(j+1) = M32(j)+(dt*(K2*M12(j)*M22(j)-K3*M12(j)*M32(j))-mu_3*M32(j));
M42(j+1) = M42(j)+(dt*(K3*M12(j)*M32(j)-K4*M12(j)*M42(j))-mu_4*M42(j));
M52(j+1) = M52(j)+(dt*(K4*M12(j)*M42(j)-K5*M12(j)*M52(j))-mu_5*M52(j));
M62(j+1) = M62(j)+(dt*(K5*M12(j)*M52(j)-K6*M12(j)*M62(j))-mu_6*M62(j));
O2(j+1) = O2(j)+(dt*(K6*M12(j)*M62(j)-Ko*M12(j)*O2(j)-mu_o*O2(j)));
P2(j+1) = P2(j)+(dt*(Ko*M12(j)*O2(j)-mu_p*P2(j)));
%with drugs depend by time
T(j+1)=T(j)+dt;
M14(j+1) = M14(j)+1./(1+exp(-T(j)));
M14(j+1) = M14(j)+(dt*(delta*M14(j+1)*(1-(M14(j+1)/gamma))-2*K1*M14(j+1)*M14(j+1)-M14(j+1)*sumter2(j+1))-((Oa-n)*K6*M14(j+1)*M64(j+1))-((Pa-Oa)*Ko*M14(j+1)*O4(j+1))-(mu_1*M14(j+1)));
M24(j+1) = M24(j)+(dt*(K1*M14(j)*M14(j)-K2*M14(j)*M24(j))-(mu_2*M24(j+1))-tmakro1-tmakro2-dtdab);
M34(j+1) = M34(j)+(dt*(K2*M14(j)*M24(j)-K3*M14(j)*M34(j))-mu_3*M34(j));
M44(j+1) = M44(j)+(dt*(K3*M14(j)*M34(j)-K4*M14(j)*M44(j))-mu_4*M44(j));
M54(j+1) = M54(j)+(dt*(K4*M14(j)*M44(j)-K5*M14(j)*M54(j))-mu_5*M54(j));
M64(j+1) = M64(j)+(dt*(K5*M14(j)*M54(j)-K6*M14(j)*M64(j))-mu_6*M64(j));
O4(j+1) = O4(j)+(dt*(K6*M14(j)*M64(j)-Ko*M14(j)*O4(j)-mu_o*O4(j)));
P4(j+1) = P4(j)+(dt*(Ko*M14(j)*O4(j)-mu_p*P4(j)));
end
subplot (2,2,1)
plot(T,M1,'k', T,M12,'r',T,M14,'b','Linewidth',2)
legend ('M1 dengan with immune', 'M1 without immune', 'M1 with drugs');
%xticks ([60 120 180 240 300 360 420 480 540 600 660 720 780 840 900 960 1020 1080])
xticks ([30 60 90 120 150 180 ])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M1')
subplot (2,2,2)
plot(T,M2,'k',T,M22,'r',T,M24,'b','Linewidth',2)
legend ('M2 with immune', 'M2 without immune','M2 with drug');
xticks ([30 60 90 120 150 180 ])
yticks ([ 0 2 4 6 8 10 12 14 16])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M2')
subplot (2,2,3)
plot(T,M3,'k',T,M32,'r',T,M34,'b','Linewidth',2)
legend ('M3 with immune', 'M3 without immune','M3 with drug');
xticks ([30 60 90 120 150 180 ])
yticks ([0 0.5 1 1.5 2 2.5 3 3.5])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M3')
subplot (2,2,4)
plot(T,M4,'k',T,M42,'r',T,M44,'b','Linewidth',2)
legend ('M4 with immune', 'M4 without immune','M4 with drug');
xticks ([30 60 90 120 150 180 ])
yticks ([0 0.1 0.2 0.3 0.4 0.5])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ('M4')
subplot (2,2,1)
plot(T,M5,'k',T,M52,'r',T,M54,'b','Linewidth',2)
legend ('M5 with immune', 'M5 without immune','M5 with drug');
xticks ([30 60 90 120 150 180 ])
yticks ([0 0.05 0.1 0.15 0.2])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("M5")
subplot (2,2,2)
plot(T,M6,'k',T,M62,'r',T,M64,'b','Linewidth',2)
legend ('M6 with immune', 'M6 without imune','M6 with drug');
xticks ([30 60 90 120 150 180 ])
yticks ([0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("M6")
subplot (2,2,3)
plot(T,O,'k',T,O2,'r',T,O4,'b','Linewidth',2)
legend ('O with immune', 'O without immune','O with drug');
xticks ([30 60 90 120 150 180 ])
yticks ([0 0.005 0.01 0.015 0.02])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("O")
subplot(2,2,4)
plot(T,P,'k',T,P2,'r',T,P4,'b','Linewidth',2)
legend ('P with immune', 'P without immune','P with drug');
xticks ([30 60 90 120 150 180 ])
%yticks ([0 2 4 6 8 10 12])
xlabel('time (day)','Fontsize',14,'FontWeight','bold')
ylabel('concentration (gr/mL)','Fontsize',14,'FontWeight','bold')
title ("P")
  8 个评论
显示 6更早的评论
Sam Chak
Sam Chak 2023-11-19
Hi @cindyawati cindyawati
I wish to see the amyloid β math equations like this (an excerpt from the Mathematical Model on Alzheimer’s disease). I'm not a coder by profession, so it is rather difficult for me to decipher the math equations from the code alone.
Furthermore, seeing the equations allows us to countercheck if the equations in your code are described correctly or not.
cindyawati cindyawati
cindyawati cindyawati 2023-11-19
This is my equation for Amyloid beta that I use. and for the drugs code you are right, that is my code that i use (but i only use ODE not ODP) @Sam Chak

请先登录，再进行评论。

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

采纳的回答

Sam Chak
Sam Chak 2023-11-19
Ran in:
Hi @cindyawati cindyawati
Please check if the differential equations and the initial values are enterer correctly or not.
% Call ode45 solver
tspan = [0 180]; % integration time interval
x0 = [10 0 0 0 0 0 0 0]; % initial values
[t, x] = ode45(@drugODE, tspan, x0);
% Plot results
subplot(221)
plot(t, x(:,5)), grid on, title('M_{5}')
subplot(222)
plot(t, x(:,6)), grid on, title('M_{6}')
subplot(223)
plot(t, x(:,7)), grid on, title('O')
subplot(224)
plot(t, x(:,8)), grid on, title('P')
% Differential equations
function dxdt = drugODE(t, x)
% definitions
dxdt = zeros(8, 1);
M1 = x(1);
M2 = x(2);
M3 = x(3);
M4 = x(4);
M5 = x(5);
M6 = x(6);
O = x(7);
P = x(8);
% Parameters
delta = 50;
gamma = 75;
K1 = 1e-4;
K2 = 5e-4;
K3 = 1e-3;
K4 = 5e-3;
K5 = 1e-2;
K6 = 5e-2;
Ko = 0.1;
n = 6;
Oa = 10;
Pa = 100;
mu1 = 1e-3;
mu2 = 1e-3;
mu3 = 1e-3;
mu4 = 1e-3;
mu5 = 1e-3;
mu6 = 1e-3;
muo = 1e-4;
mup = 1e-5;
% ODEs
dxdt(1) = delta*M1*(1 - M1/gamma) - 2*K1*M1^2 - M1*(K2*M2 + K3*M3 + K4*M4 + K5*M5) - (Oa - n)*K6*M1*M6 - (Pa - Oa)*Ko*M1*O - mu1*M1;
dxdt(2) = K1*M1*M1 - K2*M1*M2 - mu2*M2;
dxdt(3) = K2*M1*M2 - K3*M1*M3 - mu3*M3;
dxdt(4) = K3*M1*M3 - K4*M1*M4 - mu4*M4;
dxdt(5) = K4*M1*M4 - K5*M1*M5 - mu5*M5;
dxdt(6) = K5*M1*M5 - K6*M1*M6 - mu6*M6;
dxdt(7) = K6*M1*M6 - Ko*M1*O - muo*O;
dxdt(8) = Ko*M1*O - mup*P;
end
  3 个评论
显示 1更早的评论
Sam Chak
Sam Chak 2023-11-19
Hi @cindyawati
It is possible to use the user-defined Euler (ode1) solver if the step size is made arbitrarily small. As mentioned earlier, the Euler method is known for propagating errors at every time step. With the built-in ode45 solver, my only concern is ensuring that I've entered the equations and parameters correctly. In contrast, with Euler, one needs to ensure the correctness of the algorithm, not to mention the absence of tools to verify the accuracy of the results.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by