find equilirium point with all parameter known

4 次查看（过去 30 天）
Latifah
Latifah 2023-4-2
回答： Torsten 2023-4-3
Hello
I have a mathematical model for degree 7 non-linear system. I have calculated the equilibrium points for a disease-free system (x2=0 and x3=0).
However, it is not possible to find endemic equilibrium points manually, so it is necessary to use numerical methods. but I don't know how to use matlab's feature to get the endemic equilibrium point (x2≠0 , x3≠0) .
can you help me change this script, here i already have the parameter value.
Pg=0.0068;
Ps=0.012;
Pr=0.006;%0.012 0.002 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,y0]=ode45(@glioma0,[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 = glioma0(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.0;
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));
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 个）

Torsten
Torsten 2023-4-3
I don't know if this solves your problem, but of course you have to reverse the order of the commands here:
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));
You cannot use the if-statement using hdot(4) before it is computed.

类别

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

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by