trouble in plotting file differential equation contain function

9 Feb. 2023
1 Antwort
Antwort akzeptiert
Aktualisiert 9 Feb. 2023
11 Ansichten (30 Tage)
Ran in:

0 Stimmen

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 Kommentar
-1 ältere Kommentare anzeigen -1 ältere Kommentare ausblenden

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

Melden Sie sich an, um zu kommentieren.

Melden Sie sich an, um diese Frage zu beantworten.

 Akzeptierte Antwort

Latifah
Latifah am 9 Feb. 2023

0 Stimmen

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 Kommentare
3 ältere Kommentare anzeigen 3 ältere Kommentare ausblenden

Torsten
Torsten am 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 am 9 Feb. 2023
Bearbeitet: Torsten am 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 am 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 am 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.
Latifah
Latifah am 9 Feb. 2023
okay thank you so much

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Loops and Conditional Statements finden Sie in Hilfe-Center und File Exchange

Produkte

Version

R2021a

Tags

Gefragt:

Latifah
Latifah
am 9 Feb. 2023

Kommentiert:

Latifah
Latifah
am 9 Feb. 2023

Community Treasure Hunt

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

Start Hunting!

Translated by