trouble in plotting file differential equation contain function

9 Feb 2023
1 Answer
Answer Accepted
Updated 9 Feb 2023
11 Views (30 days)
Ran in:

0 votes

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 Comment
Show -1 older comments Hide -1 older comments

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

Sign in to comment.

Sign in to answer this question.

 Accepted Answer

Latifah
Latifah on 9 Feb 2023

0 votes

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 Comments
Show 3 older comments Hide 3 older comments

Torsten
Torsten on 9 Feb 2023
@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 on 9 Feb 2023
Edited: Torsten on 9 Feb 2023
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 on 9 Feb 2023
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 on 9 Feb 2023
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Products

Release

R2021a

Tags

Asked:

Latifah
Latifah
on 9 Feb 2023

Commented:

Latifah
Latifah
on 9 Feb 2023

Community Treasure Hunt

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

Start Hunting!