I have mathematical model with 7 equation . and its difficult to find endemic equilibrium point with analitical , i want to solve with numerical method
6 次查看(过去 30 天)
显示 更早的评论
can you help me to find the endemic equilibrium point , with all the parameter known
close all; clear all; clc;
Pg=0.0068;
Ps=0.012;
Pr=0.006;
pe=0.002;
beta1=0.00702;
beta2=0.00702;
beta3=0.00702;
d4=1.7143;
u=0.05;
a1=1.10;
a2=4.6205;
a4=4.6666;
a5=1.10;
c1=0.0002;
c2=0.032;
c4=0.032;
c5=0.0002;
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.01813;
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 = 1000;
options = [];
%options=odeset('Abstol', 1e-6, 'Reltol', 1e-6);
[t,y1]=ode45(@glioma1,[t0:0.1:tf],X0,options,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,u,d10,d11,d12,d20,d21,d22,d50,d51,d52,d4,tau,mu, beta3,psi, delta, gamma,alpha, Pr);
function hdot = glioma1(t,x,Pg,Ps,pe, beta1, beta2,a1, a2, a4, a5, c1,c2,c4,c5,u,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));
v=0.6;
Phi=3.3*10^-3;
H=0;
if x(6)>0 || x(7)>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));
hdot(3)=Pr*x(3)*(1-(x(2)+x(3))/(1+tau*x(4)))-beta3*x(1)*x(3)+u*x(2)*H-v*H*x(3);
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 个评论
Walter Roberson
2023-4-12
The mathematics of ode45 is not valid when you have discontinuities in the first or second deriviatives of the equations. You need to construct event functions to detect each condition under which you are changing your equations -- detect x(6) or x(7) changing between negative and positive, detect hdot(1) changing between negative and positive.
Question: what should your hdot(5) be if hdot(1) is exactly 0?
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!