trouble in plotting file differential equation contain function

2023 2 9
1 个回答
回答已采纳
更新时间：2023 2 9
8 次查看（30 天）
Ran in:

0 个投票

Ran in:
I want to plot my script, glioma_if but it doesn't work, can you help me to plot?
close all; clear all; clc;
Pg=0.0068;
Ps=0.012;
Pr=0.002;%0.012 0.002 0.006
pe=0.002;
beta1=0.018;
beta2=0.0018;
beta3=0.0018;
d4=1.7143;
a1=1;
a2=1;
a4=4.6666;
a5=1;
c1=0.0002;
c2=0.032;
c4=0.032;
c5=0.0002;
v=2.0;
d10=1.2*10^-7;
d11=4.2*10^-8;
d12=1.0*10^-7;
d20=0.2051;
d21=0.00431;
d22=19.4872;
d50=1.2*10^-7;
d51=4.2*10^-8;
d52=1.0*10^-7;
tau=0.1615;
mu=0.00371;
Phi=2.0;
psi=0.2;
delta=2.4*10^-4;
gamma=0.136;
%t = linspace(0,0.1,100)';
%tspan = [0 1000];
X0 = [1.0 0.01 0.0 1.0 0.0];
t0 = 0;
tf = 800;
options = [];
%options=odeset('Abstol', 1e-6, 'Reltol', 1e-6);
[t,y0]=ode45(@glioma0,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
Index exceeds the number of array elements. Index must not exceed 5.

Error in solution>glioma0 (line 95)
if (x(6)>0)

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
[t,y01]=ode45(@glioma01,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
[t,y1]=ode45(@glioma1,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
figure
subplot(1,1,1)
%plot(t,y0(:,1),'b','LineWidth',2);
plot(t,y0(:,1),'b',t,y01(:,1),'r',t,y1(:,1),'k');
title('g(t)');
%axis([0 750 0 1]);
%subplot(4,1,2)
%plot(t,y0(:,2),'r','LineWidth',2);
%plot(t,y0(:,2),'b',t,y01(:,2),'r',t,y1(:,2),'k');
%title('s(t)');
%axis([0 750 0 1]);
%subplot(4,1,3)
%plot(t,y0(:,3),'g','LineWidth',2);
%plot(t,y0(:,3),'b',t,y01(:,3),'r',t,y1(:,3),'k');
%title('r(t)');
%axis([0 750 0 1]);
%subplot(4,1,4)
%plot(t,y0(:,4),'b',t,y01(:,4),'r',t,y1(:,4),'k');%'LineWidth',2);
%title('n(t)');
%axis([0 750 0 1]);
%subplot(5,1,5)
%plot(t,y(:,4),'k','LineWidth',2);
%plot(t,y(:,5),'k','LineWidth',2);
%title('q(t)');
%axis([0 750 0 1]);
%figure
%plot(t,y0(:,5),'k','LineWidth',2);
%title('q(t)');
function hdot = glioma0(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.0;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end
function hdot = glioma01(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.001;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end
function hdot = glioma1(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.01;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Torsten
Torsten 2023-2-9
If you have 7 differential equations, you have to supply 7 initial conditions, not 5.

请先登录，再进行评论。

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

 采纳的回答

Latifah
Latifah 2023-2-9

0 个投票

I've added the initial condition, but it still error , can you fix it?
close all; clear all; clc;
Pg=0.0068;
Ps=0.012;
Pr=0.002;%0.012 0.002 0.006
pe=0.002;
beta1=0.018;
beta2=0.0018;
beta3=0.0018;
d4=1.7143;
a1=1;
a2=1;
a4=4.6666;
a5=1;
c1=0.0002;
c2=0.032;
c4=0.032;
c5=0.0002;
v=2.0;
d10=1.2*10^-7;
d11=4.2*10^-8;
d12=1.0*10^-7;
d20=0.2051;
d21=0.00431;
d22=19.4872;
d50=1.2*10^-7;
d51=4.2*10^-8;
d52=1.0*10^-7;
tau=0.1615;
mu=0.00371;
Phi=2.0;
psi=0.2;
delta=2.4*10^-4;
gamma=0.136;
%t = linspace(0,0.1,100)';
%tspan = [0 1000];
X0 = [1.0 0.01 0.0 1.0 0.0 0.01 0.01];
t0 = 0;
tf = 800;
options = [];
%options=odeset('Abstol', 1e-6, 'Reltol', 1e-6);
[t,y0]=ode45(@glioma0,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
[t,y01]=ode45(@glioma01,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
[t,y1]=ode45(@glioma1,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
figure
subplot(1,1,1)
%plot(t,y0(:,1),'b','LineWidth',2);
plot(t,y0(:,1),'b',t,y01(:,1),'r',t,y1(:,1),'k');
title('g(t)');
%axis([0 750 0 1]);
%subplot(4,1,2)
%plot(t,y0(:,2),'r','LineWidth',2);
%plot(t,y0(:,2),'b',t,y01(:,2),'r',t,y1(:,2),'k');
%title('s(t)');
%axis([0 750 0 1]);
%subplot(4,1,3)
%plot(t,y0(:,3),'g','LineWidth',2);
%plot(t,y0(:,3),'b',t,y01(:,3),'r',t,y1(:,3),'k');
%title('r(t)');
%axis([0 750 0 1]);
%subplot(4,1,4)
%plot(t,y0(:,4),'b',t,y01(:,4),'r',t,y1(:,4),'k');%'LineWidth',2);
%title('n(t)');
%axis([0 750 0 1]);
%subplot(5,1,5)
%plot(t,y(:,4),'k','LineWidth',2);
%plot(t,y(:,5),'k','LineWidth',2);
%title('q(t)');
%axis([0 750 0 1]);
%figure
%plot(t,y0(:,5),'k','LineWidth',2);
%title('q(t)');
function hdot = glioma0(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.0;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end
function hdot = glioma01(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.001;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end
function hdot = glioma1(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.01;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Torsten
Torsten 2023-2-9
@Latifah comment moved here:
Error using alpha
Too many output arguments.
Error in glioma_if>glioma0 (line 115)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in glioma_if (line 45)
[t,y0]=ode45(@glioma0,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5,
c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, Phi, psi, delta, gamma, Pr);
Torsten
Torsten 2023-2-9
编辑：Torsten 2023-2-9
Do you see a variable with name "alpha" somewhere in your code ? I cannot find it.
And this cannot be correct:
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
You use hdot(4) in an if-statement to compute hdot(3) although hdot(4) is not yet computed.
Latifah
Latifah 2023-2-9
finally got it to work thankyou, but I think it's still wrong. What if I just swap positions? can it be true? can you help me so that it can be calculated correctly?
close all; clear all; clc;
Pg=0.0068;
Ps=0.012;
Pr=0.002;%0.012 0.002 0.006
pe=0.002;
beta1=0.018;
beta2=0.0018;
beta3=0.0018;
d4=1.7143;
a1=1;
a2=1;
a4=4.6666;
a5=1;
c1=0.0002;
c2=0.032;
c4=0.032;
c5=0.0002;
v=2.0;
d10=1.2*10^-7;
d11=4.2*10^-8;
d12=1.0*10^-7;
d20=0.2051;
d21=0.00431;
d22=19.4872;
d50=1.2*10^-7;
d51=4.2*10^-8;
d52=1.0*10^-7;
tau=0.1615;
mu=0.00371;
psi=0.2;
delta=2.4*10^-4;
gamma=0.136;
alpha=2.0;
%t = linspace(0,0.1,100)';
%tspan = [0 1000];
X0 = [1.0 0.01 0.0 1.0 0.0 0.01 0.01];
t0 = 0;
tf = 800;
options = [];
%options=odeset('Abstol', 1e-6, 'Reltol', 1e-6);
[t,y0]=ode45(@glioma0,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma, alpha, Pr);
[t,y01]=ode45(@glioma01,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3, psi, delta, gamma, alpha, Pr);
[t,y1]=ode45(@glioma1,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);
figure
subplot(1,1,1)
%plot(t,y0(:,1),'b','LineWidth',2);
plot(t,y0(:,1),'b',t,y01(:,1),'r',t,y1(:,1),'k');
title('g(t)');
%axis([0 750 0 1]);
%subplot(4,1,2)
%plot(t,y0(:,2),'r','LineWidth',2);
%plot(t,y0(:,2),'b',t,y01(:,2),'r',t,y1(:,2),'k');
%title('s(t)');
%axis([0 750 0 1]);
%subplot(4,1,3)
%plot(t,y0(:,3),'g','LineWidth',2);
%plot(t,y0(:,3),'b',t,y01(:,3),'r',t,y1(:,3),'k');
%title('r(t)');
%axis([0 750 0 1]);
%subplot(4,1,4)
%plot(t,y0(:,4),'b',t,y01(:,4),'r',t,y1(:,4),'k');%'LineWidth',2);
%title('n(t)');
%axis([0 750 0 1]);
%subplot(5,1,5)
%plot(t,y(:,4),'k','LineWidth',2);
%plot(t,y(:,5),'k','LineWidth',2);
%title('q(t)');
%axis([0 750 0 1]);
%figure
%plot(t,y0(:,5),'k','LineWidth',2);
%title('q(t)');
function hdot = glioma0(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.0;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end
function hdot = glioma01(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.001;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end
function hdot = glioma1(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,v,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);
% The time-dependent term is A * sin(w0 * t - theta)
%H = 1*sin(-1*t);
%G = x(1)
hdot=zeros(7,1);
%H= heaviside(x(5));
u=0.01;
Phi=200.0;
H=0;
if (x(6)>0)
H = 1;
%else%if (x(5)<=0);
% H = 0;
end
hdot(1)=Pg*x(1)*(1-x(1))-beta1*x(1)*(x(2)+x(3))-(d10+d11*x(4)+d12*x(7))*(x(1)*x(6))/(a1+x(1));
hdot(2)=Ps*x(2)*(1-(x(2)+x(3))/1+tau*x(4))-beta2*x(1)*x(2)-u*x(2)*H-(d20+d21*x(4)+d22*x(7))*(x(2)*x(6))/(a2+x(2));
if (hdot(4)<0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H+v*hdot(4)*x(3);
elseif (hdot(4)>0)
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/1+tau*x(4))-beta3*x(1)*x(3)+u*x(2)*H;
end
hdot(4)=mu*(x(2)+x(3))+pe*x(4)*(1-x(4))-d4*((x(4)*x(7))/a4+x(4));
if (hdot(1)<0)
hdot(5)=alpha*hdot(1)*x(5)-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
elseif (hdot(1)>0)
hdot(5)=-(d50+d51*x(4)+d52*x(7))*(x(5)*x(6))/a5+x(5);
end
hdot(6)=Phi-(psi+c1*x(1)/(a1+x(1))+c2*x(2)/(a2+x(2))+c5*x(5)/(a5+x(5)))*x(6);
hdot(7)=delta-(gamma+c4*x(4)/(a4+x(4)))*x(7);
%hdot = hdot';
%hdot=[hdot(:);H,h1];
% To make xdot a column
% End of FUN1.M
end
Torsten
Torsten 2023-2-9
What if I just swap positions? can it be true?
Yes, of course. hdot(4) must be calculated before you use it in the if-statement to derive hdot(3).
But I would really be surprised if the integrator wouldn't come into trouble with all your if-statements in the function routines.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

产品

版本

R2021a

标签

Community Treasure Hunt

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

Start Hunting!

Translated by