AMCS 252 A4
Zein Elrez
March 24, 2025
1 Question 1
Please refer to the scanned written solution at the end of this report.
2 Question 2
Please refer to the scanned written solution at the end of this report.
3 Question 3
Please refer to the scanned written solution at the end of this report.
Next, we implement our scheme on python, in order to approximate the
exact solution. The following is the code:
1 import numpy as np
2 import matplotlib . pyplot as plt
3
4 timeStep = 0.1
5 lamb = -0.6
6 T = 10
7 nu = 0.5
8
9 tVals = np . arange (0 ,10.1 , timeStep )
10 z = timeStep * lamb
11
12 x1 = z + np . sqrt ( z **2 + 1)
13 x2 = z - np . sqrt ( z **2 + 1)
14
15
16 c1 = nu *( 1 + np . sqrt ( z **2+1) ) /(2* np . sqrt ( z **2 + 1) )
1
17 c2 = nu *( np . sqrt ( z **2+1) - 1 ) /(2* np . sqrt ( z **2 + 1) )
18
19 nVals = np . arange ( T / timeStep +1)
20
21 approxSol = c1 * x1 ** nVals + c2 * x2 ** nVals
22
23 exactSol = nu * np . exp ( lamb * tVals )
24
25
26 plt . xlabel ( " t " )
27 plt . ylabel ( " u ( t ) " )
28 plt . title ( " Approx . Sol . vs . Theo . Sol . " )
29 plt . plot ( tVals , approxSol , label = " Approximate Solution " )
30 plt . plot ( tVals , exactSol , label = " Theoretical Solution " )
31 plt . legend ()
32 plt . show ()
Looking at the plot of the theoretical solution vs. approximate solution
below, we notice that the approximate solution begins to diverge for larger
values of T.
Figure 1: Approximate solution vs. Theoretical solution
Notice the root x2 is causing a divergence in the approximate solution. The
reason is that x2 may be negative, so, the term c2 x2 will fluctuate about the
time axis. More specifically, when n is odd, xn2 is negative and when n is even
xn2 is positive (or vice versa). Let’s plot c1 xn1 , and c2 xn2 to observe the behavior
2
of each fundamental solution. The following is the python code that implements
it:
1 c1x1 = lambda n : c1 * x1 ** n
2 c2x2 = lambda n : c2 * x2 ** n
3
4 plt . xlabel ( " n " )
5 plt . ylabel ( " U ^ n " )
6 plt . title ( " c1x1 ^ n vs . c2x2 ^ n " )
7 plt . plot ( nVals , [ c1x1 ( n ) for n in nVals ] , label = " c1x1 ^ n " )
8 plt . plot ( nVals , [ c2x2 ( n ) for n in nVals ] , label = " c2x2 ^ n " )
9
10 plt . legend ()
11 plt . show ()
And we get the following plot as an output :
Figure 2: c1 xn1 vs. c2 xn2
Looking at c1 xn1 vs. c2 xn2 , we notice that the plot for c1 xn1 is smooth (con-
vergent), but c2 xn2 fluctuates frequently √(divergent). Hence, we should try to
2
eliminate that term. Looking at c2 = η( 2√z z+1−1)2 +1
, as z →
− 0, then c2 → − 0. So,
n
c2 x 2 →
− 0. So, let us pick a small z; for example, z = -0.06, and let us check the
approximate solution vs. theoretical solution for this case. We implement the
aforementioned with the following code :
1 timeStep = 0.01
2 lamb = -0.6
3
3 T = 10
4 nu = 0.5
5
6 tVals = np . arange (0 ,10.01 , timeStep )
7 z_star = timeStep * lamb
8
9 x1 = z_star + np . sqrt ( z_star **2 + 1)
10 x2 = z_star - np . sqrt ( z_star **2 + 1)
11
12 nVals = np . arange ( T / timeStep +1)
13
14 c1 = nu *( 1 + np . sqrt ( z_star **2+1) ) /(2* np . sqrt ( z_star **2 + 1) )
15 c2 = nu *( np . sqrt ( z_star **2+1) - 1 ) /(2* np . sqrt ( z_star **2 + 1) )
16
17 approxSol = c1 * x1 ** nVals + c2 * x2 ** nVals
18
19 exactSol = nu * np . exp ( lamb * tVals )
20
21
22 plt . xlabel ( " t " )
23 plt . ylabel ( " u ( t ) " )
24 plt . title ( " Approx . Sol . vs . Theo . Sol . " )
25 plt . plot ( tVals , approxSol , label = " Approximate Solution " )
26 plt . plot ( tVals , exactSol , label = " Theoretical Solution " )
27 plt . legend ()
28 plt . show ()
And we get the following plot as an output:
Figure 3: Approximate solution vs. Theoretical solution
4
Notice from the figure above, that the approximate solution converges to
that of the theoretical solution. Hence, picking a very small value for z solves
the paradox.
4 Question 4
4.1 Part (i)
Let x(t) : Susceptible and y(t) : infections. Then the SIR model is as follow :
x′ (t) = −βxy, and y ′ (t) = βxy − γy, where β is the contract/transmission rate
and γ1 is the time of infectiousness. Using Euler’s Method: xn+1 = xn − βxn yn
and yn+1 = βxn yn − γyn , we solve the IVP with the following code:
1 import numpy as np
2 import matplotlib . pyplot as plt
3
4 def EulerSolve ( x0 , y0 , timeStep , Beta , Gamma ) :
5 x = x0
6 y = y0
7 h = timeStep
8
9 x_vals = [ x0 ]
10 y_vals = [ y0 ]
11 t = [0]
12 i =1
13
14 while True :
15 x = x - Beta * h * x * y
16 y = y + h *( Beta * x * y - Gamma * y )
17
18 x_vals . append ( x )
19 y_vals . append ( y )
20
21 t . append ( i * h )
22
23 if ( abs ( x_vals [ len ( x_vals ) -1] - x_vals [ len ( x_vals ) -2]) < 10** -6
and abs ( y_vals [ len ( y_vals ) -1] - y_vals [ len ( y_vals ) -2]) < 10** -6)
:
24 break
25 i += 1
26
27 return x_vals , y_vals , t
28
29 Beta = 0.3
30 Gamma = 0.1
31 timeStep = 0.01
32
33 x0 = 0.99
34 y0 = 0.01
35
36 xVals , yVals , tVals = EulerSolve ( x0 , y0 , timeStep , Beta , Gamma )
37 newxVals = [ i * 100 for i in xVals ]
5
38 newyVals = [ i * 100 for i in yVals ]
39
40 plt . xlabel ( " Time ( Days ) " )
41 plt . ylabel ( " % of population " )
42
43 plt . title ( " Susceptible and Infectious given Beta = 0.3 " )
44 plt . plot ( tVals , newxVals , label = " Susceptible " )
45 plt . plot ( tVals , newyVals , label = " Infectious " )
46 plt . legend ()
47 plt . show ()
48
49 maxInfected = print ( f " Given Beta = { Beta } , Max infected is : { round
( max ( yVals ) *100 ,2) }%. " )
The following is the plot of the infectiousness and susceptibility generated
by the code (these are the indicators of the proportions of the population that
are either susceptible or infected at any point in time over the duration of the
pandemic):
Figure 4: Susceptibility vs. Infectiousness given β = 0.3
Notice that the values of x(t), y(t) are nearly constant after about 90 days,
so the pandemic lasts roughly 90 days. In addition, following the output of
the code, the maximum infected at any time instance in the duration of the
pandemic is 30.36% of the population. (this can be checked by running the
above code).
6
4.2 Part (ii)
To simulate social distancing, quarantining, and so forth, we drop the trans-
mission rate β to 0.15. We re-implement what we did in Part(i) with the same
code, but changing the value of β to the aforementioned.
Figure 5: Susceptibility vs. Infectiousness given β = 0.15
The values of x(t), y(t) are nearly constant after about 180 days, so, the
pandemic last for roughly 180 days. In addition, the maximum infected at any
time instance in the duration of the pandemic is 6.9% of the population (this
can be checked by running the above code with β = 0.15). Notice that the
curves in this case are much flatter than β = 0.3.
Let’s continue by dropping the transmission rate β to 0.05. Subsequently,
the above code (in Part (i)) generates the following plot :
7
Figure 6: Susceptibility vs. Infectiousness given β = 0.05
Notice that the values of x(t), y(t) become constant almost instantaneously
(after 1-2 days). In this case, the pandemic would have not had much of an
impact on the population. However, the maximum infected at any time instance
in the duration of the pandemic is 1% of the population (this can be checked
by running the above code with β = 0.05).
In conclusion, decreasing the rate of transmission leads to a decreased impact
of the pandemic on the population. This suggests that we are using a correct
mathematical model to simulate the spreading of infectious diseases across a
population.
8
►
' /6.),
L) r tot 7
,. _ -= e
/t90
.e - .
I,; J(
( 1 I I
I
j •
\ I
A
•
/1!/J~ .
I ' '
) '.
r \ ,
' \
I •
I \, •
l • ' . .
\ I \ 1' l t· \\ I t
'· ) I
' I } '\ . )
I
I ': t
·( :.
I I'
I '. I
I
f , .. (
\
.' \ i ) ,, I
!': 'I
'· :
j " '
:f )
,.
j (
I
) '·
\I; (
\ \ i'( I ,! ' t)
I
V
..,.
t ()
I I () 4 II
L '2.-
l ~ r + r _, 0 2 •l
-12.
~·
e, 0 0 --3a -',D
, 0 \ \l
f) 1 -1 ... ,'l,
0 I 1
-
l- 'l.1 " - - J
1
:2Jc~t1
( I
,'
I ( r'
I
f r
I I ,I '
I t ' ,.,
f
/
/ /
I
('
'i
)
/' (
( . ' ,,
t _,
'
, (,I
I
I
1 /
I
(
I !
·/✓ 1 /,c, .i
"r1
