Practical 4: Growth and Decay Model
(Exponential case)
Question 1: Plot the solution of the
differential equation
dY
dt
=0.05 Y, where Y(0) = 1000 over the
domain [0,10].
In[42]:= sol = DSolve[{D[Y[t], t] ⩵ 0.05 * Y[t],
Y[0] ⩵ 1000}, Y[t], t]
Out[42]= Y[t] → 1000. ⅇ0.05 t
In[43]:= Plot[Y[t] /. sol, {t, 0, 10},
PlotStyle → {Thick, Red},
PlotLabel → "Growth Model"]
Growth Model
1600
1500
1400
Out[43]=
1300
1200
1100
2 4 6 8 10
Question 2: Plot the solution of the
differential equation
2 Practical 4 (Growth and Decay Model).nb
dY
dt
= -0.05 Y, where Y(0) = 5000 over the
domain [0,10].
In[66]:= sol1 = DSolve[{D[Y[t], t] ⩵ -0.05 * Y[t], Y[0] ⩵ 5000}, Y[t], t]
Out[66]= Y[t] → 5000. ⅇ-0.05 t
In[52]:= Plot[Y[t] /. sol1, {t, 0, 10},
PlotStyle → {Magenta, Dashed}, PlotLabel → "Decay Model"]
Decay Model
5000
4500
Out[52]=
4000
3500
2 4 6 8 10
Question 3: Suppose that the population of a certain
country grows at annual rate of 2% . If current
population is 3 million what will be the population in
(i) 10 years
(ii) 20 years
(iii)50 years
In[53]:= sol2 = DSolve[{D[Y[t], t] ⩵ 0.02 * Y[t], Y[0] ⩵ 3 000 000}, Y[t], t]
Out[53]= Y[t] → 3. × 106 ⅇ0.02 t
Practical 4 (Growth and Decay Model).nb 3
In[55]:= Plot[Y[t] /. sol2, {t, 0, 100}, PlotStyle → {Red, Thick},
PlotLabel → "Population of a Country"]
Population of a Country
2.0 × 107
1.5 × 107
Out[55]=
1.0 × 107
5.0 × 106
20 40 60 80 100
In[56]:= Y[t] /. sol2 /. t → 10
Out[56]= 3.66421 × 106
The population of the country in 10 year would be 3.66421 x
106.
In[57]:= Y[t] /. sol2 /. t → 20
Out[57]= 4.47547 × 106
The population of the country in 10 year would be 4.47547 x
106.
In[59]:= Y[t] /. sol2 /. t → 50
Out[59]= 8.15485 × 106
The population of the country in 10 year would be 8.15485 x 106.
Question 4: A city has population of 25, 000 in 1960 and
30, 000 in 1970. Assuming the population to grow
exponentially, what is the expected population in 2020.
In[60]:= sol3 = DSolve[{D[Y[t], t] ⩵ r * Y[t], Y[0] ⩵ 25 000}, Y[t], t]
Out[60]= Y[t] → 25 000 ⅇr t
4 Practical 4 (Growth and Decay Model).nb
In[61]:= Y[t] /. sol3 /. t → 10
Out[61]= 25 000 ⅇ10 r
In[63]:= NSolve25 000 * ⅇ10*r ⩵ 30 000, r, Reals
Out[63]= {{r → 0.0182322}}
In[64]:= X[t_] = Y[t] /. sol3 /. r → 0.0182322
Out[64]= 25 000 ⅇ0.0182322 t
In[65]:= X[60]
Out[65]= {74 649.8}
The Population of the city will be 74649.8 in 2020.
Practical 5
Lake Pollution Model
General Compartmental Model
This leads us to the word equation, for the mass of pollutant in
the lake,
Lake Pollution Model.nb ��� 5
6 ��� Lake Pollution Model.nb
Question3. Solve Example 2.8
F F
��� ��� sol2 = DSolveC '[t] ⩵ cin - C[t], C[0] ⩵ c0 , C[t], t
V V
Ft Ft
c0 - cin + ⅇ cin
-
���� �� t → ⅇ V V
��� ��� C[t] /. sol2 /. cin → 0
Ft
c0
-
���� �� ⅇ V
t
Solveⅇ c0 ⩵ 0.05 * c0 , t, Reals
-F*
��� ��� V
������ ����� ��� ������ �� ����� ��� ������ ���� ������� ������������ ��� ������ ��� �������� �� ������� �
������������� ����� ������ ��� ������������ ��� �������
2.99573 V
���� �� t →
F
Lake Pollution Model.nb ��� 7
Lake Erie
V=458 x 109 m3 , F = 1.75 x 1011 m3/year
F F
��� ��� sol3 = DSolveC '[t] ⩵ cin - C[t], C[0] ⩵ c0 , C[t], t
V V
Ft Ft
c0 - cin + ⅇ cin
-
���� �� t → ⅇ V V
��� ��� C[t] /. sol3 /. cin → 0, F → 1.75 * 1011 , V → 458 * 109
���� �� ⅇ-0.382096 t c0
��� ��� Solveⅇ-0.382096 t c0 ⩵ 0.05 * c0 , t, Reals
������ ����� ��� ������ �� ����� ��� ������ ���� ������� ������������ ��� ������ ��� �������� �� ������� �
������������� ����� ������ ��� ������������ ��� �������
���� �� {{t → 7.84026}}
��� ��� Plot[{Evaluate[Exp[- 0.382096 * t] * c0 /.
{{c0 → 150 * 10 ^ 9}, {c0 → 140 * 10 ^ 9}, {c0 → 130 * 10 ^ 9}, {c0 → 120 * 10 ^ 9}}],
0.05 * 150 * 10 ^ 9}, {t, 6, 9}, PlotLegends → "Expressions"]
1.4 × 1010
1.2 × 1010
150 000 000 000 ⅇ-0.382096 t
140 000 000 000 ⅇ-0.382096 t
10
1.0 × 10
���� �� 130 000 000 000 ⅇ-0.382096 t
8.0 × 109 120 000 000 000 ⅇ-0.382096 t
0.05 × 150 × 109
6.0 × 109
6.5 7.0 7.5 8.0 8.5 9.0
It takes pollution of Lake Erie to reduce to 5% of its current pollution in 7.84
years.
Lake Ontario
8 ��� Lake Pollution Model.nb
V=1636 x 109 , F = 2.089 x 1011
F F
��� ��� sol4 = DSolveC '[t] ⩵ cin - C[t], C[0] ⩵ c0 , C[t], t
V V
Ft Ft
c0 - cin + ⅇ cin
-
���� �� t → ⅇ V V
��� ��� C[t] /. sol4 /. cin → 0, F → 2.089 * 1011 , V → 1636 * 109
���� �� ⅇ-0.127689 t c0
��� ��� Solveⅇ-0.127689 t c0 ⩵ 0.05 * c0 , t, Reals
������ ����� ��� ������ �� ����� ��� ������ ���� ������� ������������ ��� ������ ��� �������� �� ������� �
������������� ����� ������ ��� ������������ ��� �������
���� �� {{t → 23.4612}}
It takes pollution of Lake Ontario to reduce to 5% of its current pollution in
23.46 years.
Case Study: Lake Burley Griffin
F F
��� ��� sol5 = DSolveC '[t] ⩵ cin - C[t], C[0] ⩵ c0 , C[t], t
V V
Ft Ft
c0 - cin + ⅇ cin
-
���� �� t → ⅇ V V
��� ��� a = C[t] /. sol5 /. {cin → 3, F → 48, V → 28}
���� �� ⅇ-12 t/7 - 3 + 3 ⅇ12 t/7 + c0
��� ��� a = Solveⅇ-12 t/7 - 3 + 3 ⅇ12 t/7 + c0 ⩵ 0.05 * c0 , t, Reals
������ ����� ��� ������ �� ����� ��� ������ ���� ������� ������������ ��� ������ ��� �������� �� ������� �
������������� ����� ������ ��� ������������ ��� �������
20. - 3. + c0
���� �� t → ConditionalExpression0.583333 Log , c0 > 60. || c0 < 3.
- 60. + c0
��� ��� a /. c0 → 6
���� �� {{t → Undefined}}
Lake Pollution Model.nb ��� 9
��� ��� Plot[Evaluate[a /. {{c0 → 6}, {c0 → 5}, {c0 → 4}, {c0 → 3}, {c0 → 2}, {c0 → 1}}],
{t, 0, 8}, PlotLegends → "Expressions", PlotRange → All,
PlotStyle → {Red, Blue, Orange, Dashed, Magenta, Black}]
5 ⅇ-12 t/7 3 + 3 ⅇ12 t/7
ⅇ-12 t/7 2 + 3 ⅇ12 t/7
4
ⅇ-12 t/7 1 + 3 ⅇ12 t/7
���� ��
3
3
ⅇ-12 t/7 -1 + 3 ⅇ12 t/7
2 ⅇ-12 t/7 -2 + 3 ⅇ12 t/7
2 4 6 8
��� ��� Plot[{Evaluate[a /. {c0 → 6}], a /. {c0 → 5}, a /. {c0 → 4}, a /. {c0 → 3},
a /. {c0 → 2}, a /. {c0 → 1}}, {t, 0, 8}, PlotLegends → "Expressions",
PlotRange → All, PlotStyle → {Red, Blue, Orange, Dashed, Magenta, Black}]
5
ⅇ-12 t/7 3 + 3 ⅇ12 t/7
a /. {c0 → 5}
4
a /. {c0 → 4}
���� ��
a /. {c0 → 3}
3
a /. {c0 → 2}
2 a /. {c0 → 1}
2 4 6 8
��� ��� Plot[Evaluate[a /. {c0 → {6, 5, 4, 3, 2, 1}}],
{t, 0, 8}, PlotRange → All, PlotLegends → "Expressions"]
5 ⅇ-12 t/7 3 + 3 ⅇ12 t/7
ⅇ-12 t/7 2 + 3 ⅇ12 t/7
4
ⅇ-12 t/7 1 + 3 ⅇ12 t/7
���� ��
3
3
ⅇ-12 t/7 -1 + 3 ⅇ12 t/7
2 ⅇ-12 t/7 -2 + 3 ⅇ12 t/7
2 4 6 8
10 ��� Lake Pollution Model.nb
��� ��� sol6 =
1 + 6 * Sin[2 π t]
NDSolveC '[t] ⩵ * 10 + 10 * Cos[2 * π * t] - C[t], C[0] ⩵ 6,
28
C, {t, 0, 20}
1 + 6 * Sin[2 π t]
sol7 = NDSolveC '[t] ⩵ * 10 + 10 * Cos[2 * π * t] - C[t],
28
C[0] ⩵ 5, C, {t, 0, 20}
1 + 6 * Sin[2 π t]
sol8 = NDSolveC '[t] ⩵ * 10 + 10 * Cos[2 * π * t] - C[t],
28
C[0] ⩵ 4, C, {t, 0, 20}
1 + 6 * Sin[2 π t]
sol9 = NDSolveC '[t] ⩵ * 10 + 10 * Cos[2 * π * t] - C[t],
28
C[0] ⩵ 3, C, {t, 0, 20}
������� {{��� ���}}
���� �� C → InterpolatingFunction ������� ������
������� {{��� ���}}
���� �� C → InterpolatingFunction ������� ������
������� {{��� ���}}
���� �� C → InterpolatingFunction ������� ������
������� {{��� ���}}
���� �� C → InterpolatingFunction ������� ������
��� ��� Plot[{Evaluate[C[t] /. sol6], Evaluate[C[t] /. sol7],
Evaluate[C[t] /. sol8], Evaluate[C[t] /. sol9]}, {t, 0, 8}]
���� ��
5
2 4 6 8
Limited Growth Model without
harvesting
or
Logistic Equation
dX X
dt
= r X 1 - K , X (0) = x0
Question. Solve and plot the limit growth model with r
=1 and K=1000 by taking different initial Values.
��� ��� sol = DSolveX '[t] ⩵ X[t] * 1 - X[t] 1000, X[0] ⩵ a, X[t], t
1000 a ⅇt
���� �� X[t] →
1000 - a + a ⅇt
��� ��� Plot[
{X[t] /. sol /. {a → {50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650,
700, 750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200}},
1000}, {t, 0, 8}, PlotRange → All]
1200
1000
800
���� �� 600
400
200
2 4 6 8
��� ��� sol1 = DSolveX '[t] ⩵ r * X[t] * 1 - X[t] K, X[0] ⩵ a, X[t], t
a ⅇr t K
���� �� X[t] →
- a + a ⅇr t + K
2 ��� Limit growth model.nb
��� ��� Plot[
{X[t] /. sol1 /. {a → {50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650,
700, 750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200},
r → 1, K → 1000}, 1000}, {t, 0, 8}, PlotRange → All]
1200
1000
800
���� �� 600
400
200
2 4 6 8
��� ��� sol2 = DSolveX '[t] ⩵ r * X[t] * 1 - X[t] K, X[0] ⩵ a, X[t], t
a ⅇr t K
���� �� X[t] →
- a + a ⅇr t + K
��� ��� y[t] = X[t] /. sol2 /. {r → 1, K → 1000}
1000 a ⅇt
���� ��
1000 - a + a ⅇt
��� ��� Plot[
{y[t] /. {a → {50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700,
750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200}},
1000}, {t, 0, 8}, PlotRange → All]
1200
1000
800
���� �� 600
400
200
2 4 6 8
Limited Growth Model with
harvesting
Limit growth model.nb ��� 3
or
Logistic Equation with harvesting
dX X
dt
= r X 1 - K
- h, X (0) = x0
Question. Solve and plot the limit growth model with
harvesting h=100 , r =1 and K=1000 by taking different
initial Values.
X[t]
��� ��� sol4 = DSolveX '[t] ⩵ X[t] * 1 - - 100, X[0] ⩵ a, X[t], t
1000
3 3 3
t t t
100 1000 - 5 a + 15 a - 1000 ⅇ 5
+5aⅇ 5
+ 15 a ⅇ 5
���� �� X[t] →
3 3 3
t t t
500 + 100 15 - a - 500 ⅇ 5
+ 100 15 ⅇ 5
+aⅇ 5
��� ��� Plot[
{X[t] /. sol4 /. {a → {150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700,
750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200}}, 1000},
{t, 0, 8}, PlotRange → All, AxesOrigin → {0, 0}]
1200
1000
800
���� �� 600
400
200
2 4 6 8
X[t]
��� ��� sol5 = DSolveX '[t] ⩵ X[t] * 1 - - 500, X[0] ⩵ a, X[t], t
1000
������ ������� ��������� ��� ����� ���� �� ������ �� ���� ��������� ��� ��� �� ������ ��� ������ ��� �������� ��������
������������
1 500 - a
���� �� X[t] → 500 1 + Tan - t - 2 ArcTan ,
2 500
1 1
X[t] → 500 1 + Tan - t + 2 ArcTan - 500 + a
2 500
4 ��� Limit growth model.nb
��� ��� Plot[
{Evaluate[X[t] /. sol5[[1]] /. {a → {266, 290, 300, 350, 400, 450, 500, 550, 600,
650, 700, 750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200}}], 1000},
{t, 0, 3.5}, PlotRange → {{0, 8}, {0, 1200}}, AspectRatio → 0.9]
1200
1000
800
���� �� 600
400
200
0
0 2 4 6 8
��� ��� X[t] /. sol5[[1]] /. {a → 1150}
1 13
���� �� 500 1 + Tan - t + 2 ArcTan
2 10
1 4 1 9
��� ��� Plot500 1 + Tan - t - 2 ArcTan , 500 1 + Tan - t - 2 ArcTan ,
2 5 2 10
t 1 1
500 1 - Tan , 500 1 + Tan - t - 2 ArcTan ,
2 2 5
1 13
500 1 + Tan - t + 2 ArcTan , {t, 0, 8}, PlotRange → {{0, 8}, {0, 1300}}
2 10
1200
1000
800
���� �� 600
400
200
0
0 2 4 6 8
Limit growth model.nb ��� 5
��� ��� Plot[Tan[x], {x, 0, 10}]
6
���� ��
2 4 6 8 10
-2
-4
-6
��� ��� Plot[{X[t] /. sol5[[1]] /.
{{a → 400}, {a → 595}, {a → 600}, {a → 650}, {a → 700}, {a → 750}, {a → 800}},
1000}, {t, 0, 8}, PlotRange → {{0, 8}, {0, 1300}}]
1200
1000
800
���� �� 600
400
200
0
0 2 4 6 8
dX r
dt = - K X 2-K X + K hr )
9
r = 1, K = 10, h =
10
X[t] 9
��� ��� sol6 = DSolveX '[t] ⩵ X[t] * 1 - - , X[0] ⩵ a, X[t], t
10 10
9 - a - 9 ⅇ4 t/5 + 9 a ⅇ4 t/5
���� �� X[t] →
9 - a - ⅇ4 t/5 + a ⅇ4 t/5
6 ��� Limit growth model.nb
��� ��� Plot[Evaluate[X[t] /. sol6 /.
{a → {0.96, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 7, 8, 9, 10, 10.2, 10.5}}],
{t, 0, 8}, PlotRange → {{0, 8}, {0, 12}}, AspectRatio → 0.9]
12
10
���� �� 6
0
0 2 4 6 8
Predator-Prey Model (Lotka–Volterra predator-prey System)
��� ��� α2 = 0.5;
β = 1;
c1 = 0.01;
c2 = 0.005;
sol = NDSolve[{X '[t] ⩵ β * X[t] - c1 * X[t] * Y[t],
Y '[t] ⩵ c2 * X[t] * Y[t] - α2 * Y[t], X[0] ⩵ 200, Y[0] ⩵ 80}, {X[t], Y[t]}, {t, 20}]
������� {{��� ���}}
���� �� X[t] → InterpolatingFunction ������� ������
[t],
������� {{��� ���}}
Y[t] → InterpolatingFunction ������� ������
[t]
��� ��� Plot[{X[t] /. sol, Y[t] /. sol}, {t, 0, 20},
PlotStyle → {Black, Gray}, AxesOrigin → {0, 0}]
200
150
���� �� 100
50
5 10 15 20
Predator and Prey Model (DDT)
2 ��� Predator Prey Model.nb
sol1 = NDSolve[{X '[t] ⩵ X[t] - 0.01 * X[t] * Y[t] - 0.1 * X[t],
Y '[t] ⩵ 0.005 * X[t] * Y[t] - 0.5 * Y[t] - 0.1 * Y[t],
X[0] ⩵ 200, Y[0] ⩵ 80}, {X[t], Y[t]}, {t, 20}]
������� {{��� ���}}
���� �� X[t] → InterpolatingFunction ������� ������
[t],
������� {{��� ���}}
Y[t] → InterpolatingFunction ������� ������
[t]
��� ��� Plot[{X[t] /. sol1, Y[t] /. sol1}, {t, 0, 20},
PlotStyle → {Black, Gray}, AxesOrigin → {0, 0}]
200
150
���� �� 100
50
5 10 15 20
Predator-Prey Model (Density Dependent)
X[t]
��� ��� sol2 = NDSolveX '[t] ⩵ X[t] * 1 - - 0.01 * X[t] * Y[t], Y '[t] ⩵
1000
0.005 * X[t] * Y[t] - 0.5 * Y[t], X[0] ⩵ 200, Y[0] ⩵ 80, {X[t], Y[t]}, {t, 40}
������� {{��� ���}}
���� �� X[t] → InterpolatingFunction ������� ������
[t],
������� {{��� ���}}
Y[t] → InterpolatingFunction ������� ������
[t]
Predator Prey Model.nb ��� 3
��� ��� Plot[{X[t] /. sol2, Y[t] /. sol2}, {t, 0, 40},
PlotStyle → {Black, Gray}, AxesOrigin → {0, 0}]
200
150
���� �� 100
50
10 20 30 40
Epidemic Model of Influenza
������� beta := 2.18 * 10-3 ;
gamma := 0.44;
sol = NDSolve[{S '[t] ⩵ - beta * S[t] * T[t], T '[t] ⩵ beta * S[t] * T[t] - gamma * T[t],
S[0] ⩵ 762, T[0] ⩵ 1}, {S[t], T[t]}, {t, 0, 20}]
������� {{��� ���}}
������� S[t] → InterpolatingFunction ������� ������
[t],
������� {{��� ���}}
T[t] → InterpolatingFunction ������� ������
[t]
������� Plot[{S[t] /. sol, T[t] /. sol}, {t, 0, 20},
PlotStyle → {Black, Gray}, PlotLegends → {"Susceptible", "Infected"}]
800
600
Susceptible
������� 400
Infected
200
5 10 15 20
2 ��� Epidemic copy.nb
��� ��� a = 1 80;
b = 1 80;
beta = 5.2 * 10-4 ;
gamma = 52;
sol = NDSolveS '[t] ⩵ b * S[t] + R[t] + T[t] - beta * S[t] * T[t] - a * S[t],
T '[t] ⩵ beta * S[t] * T[t] - gamma * T[t] - a * T[t], R '[t] ⩵ gamma * T[t] - a * T[t],
S[0] ⩵ 105 , T[0] == 10, R[0] ⩵ 106 - 10 - 105 , {S[t], T[t], R[t]}, {t, 0, 20}
������� {{��� ���}}
���� �� S[t] → InterpolatingFunction ������� ������
[t],
������� {{��� ���}}
T[t] → InterpolatingFunction ������� ������
[t],
������� {{��� ���}}
R[t] → InterpolatingFunction ������� ������
[t]
��� ��� Plot[T[t] /. sol, {t, 0, 20}, PlotStyle → Black, PlotRange → {{0, 20}, {0, 1000}}]
1000
800
600
���� ��
400
200
0
0 5 10 15 20