In[]:=
ClearAll["`*"];Remove["`*"];
In[]:=
b=μ-g;​​g=0.001037052;​​μ=0.00616/12;​​
μ
1
=0.544/12;​​
η
1
=0.0002476904;​​
η
2
=0.088256814;​​
β
1
=0.005026665;​​
β
2
=0.045262083;​​
γ
1
=0;​​
γ
2
=0.91924;​​
γ
3
=0.456/12;​​ρ=0.9384;​​
ρ
1
=0;​​ω=0.080716;​​l=1;(*Iassumed*)​​h=1;(*Iassumed*)​​TIME=12;​​
In[]:=
eq={S'[t]b-(
η
1
S[t]A[t]+
β
1
S[t]c[t])-(μ+g)S[t],Y'[t]gS[t]-(
η
2
Y[t]A[t]+
β
2
Y[t]c[t])-μY[t],L'[t](
η
1
S[t]+
η
2
Y[t])A[t]+(
β
1
S[t]+
β
2
Y[t])c[t]-(μ+
μ
1
+hρ)L[t],A'[t]hρL[t]+
ρ
1
c[t]-(μ+
γ
2
)A[t]-(1-l)ωA[t],c'[t](1-l)ωA[t]-(μ+
μ
1
+
ρ
1
+
γ
3
)c[t],R'[t](1-h)
γ
1
L[t]+l
γ
2
A[t]+
γ
3
c[t]-μR[t],S[0]657824,Y[0]32331634,L[0]0,A[0]20143,c[0]1745,R[0]0}
Out[]=
{
′
S
[t]-0.000523719-0.00155039S[t]-0.00024769A[t]S[t]-0.00502667c[t]S[t],
′
Y
[t]0.00103705S[t]-0.000513333Y[t]-0.0882568A[t]Y[t]-0.0452621c[t]Y[t],
′
L
[t]-0.984247L[t]+c[t](0.00502667S[t]+0.0452621Y[t])+A[t](0.00024769S[t]+0.0882568Y[t]),
′
A
[t]0.-0.919753A[t]+0.9384L[t],
′
c
[t]0.-0.0838467c[t],
′
R
[t]0.91924A[t]+0.038c[t]-0.000513333R[t],S[0]657824,Y[0]32331634,L[0]0,A[0]20143,c[0]1745,R[0]0}
In[]:=
sol=NDSolve[eq,{S[t],Y[t],L[t],A[t],c[t],R[t]},{t,0,TIME}]
Out[]=
S[t]InterpolatingFunction
Domain: {{0.,12.}}
Output: scalar
[t],Y[t]InterpolatingFunction
Domain: {{0.,12.}}
Output: scalar
[t],L[t]InterpolatingFunction
Domain: {{0.,12.}}
Output: scalar
[t],A[t]InterpolatingFunction
Domain: {{0.,12.}}
Output: scalar
[t],c[t]InterpolatingFunction
Domain: {{0.,12.}}
Output: scalar
[t],R[t]InterpolatingFunction
Domain: {{0.,12.}}
Output: scalar
[t]
In[]:=
Plot[Evaluate[{S[t],Y[t],L[t],A[t],c[t],R[t]}/.sol],{t,0,TIME},PlotRange->Full]
Out[]=
2
4
6
8
10
12
5.0×
6
10
1.0×
7
10
1.5×
7
10
2.0×
7
10
2.5×
7
10
3.0×
7
10
(*----------------------------------------------------*)
In[]:=
Plot[Evaluate[{S[t]}/.sol],{t,0,TIME}]
Out[]=
2
4
6
8
10
12
-0.00010
-0.00005
0.00005
In[]:=
Plot[Evaluate[{Y[t]}/.sol],{t,0,TIME}]
Out[]=
2
4
6
8
10
12
-4×
-11
10
-3×
-11
10
-2×
-11
10
-1×
-11
10
1×
-11
10
2×
-11
10
In[]:=
Plot[Evaluate[{L[t]}/.sol],{t,0,TIME},PlotRange->All]
Out[]=
2
4
6
8
10
12
5.0×
6
10
1.0×
7
10
1.5×
7
10
2.0×
7
10
2.5×
7
10
3.0×
7
10
In[]:=
Plot[Evaluate[{A[t]}/.sol],{t,0,TIME},PlotRange->All]
Out[]=
2
4
6
8
10
12
2.0×
6
10
4.0×
6
10
6.0×
6
10
8.0×
6
10
1.0×
7
10
1.2×
7
10
In[]:=
Plot[Evaluate[{c[t]}/.sol],{t,0,TIME},PlotRange->All]
Out[]=
2
4
6
8
10
12
600
800
1000
1200
1400
1600
1800
In[]:=
Plot[Evaluate[{R[t]}/.sol],{t,0,TIME},PlotRange->All]
Out[]=
2
4
6
8
10
12
5.0×
6
10
1.0×
7
10
1.5×
7
10
2.0×
7
10
2.5×
7
10
3.0×
7
10