COSC1003/1903 Assignment
COSC1003/1903 Computational Science
Assignment — due Monday October 15th, 11:59pm
The assignment consists of two questions (COSC1003) or three questions (COSC1903). It is worth
10% of your total assessment for this Unit of Study. You will be marked on the correctness and
quality of the code, as well as your written responses to the questions. You may submit your solutions
independently or as a pair. If you submit as a pair you will both receive the same mark. Both members
of the pair must be in the same stream (1003 or 1903).
Download and unzip the assignment zip file (codes.zip) from the course Canvas site. You
should write your solutions in these template files. All programs should include a one-paragraph
description of how they work as a comment at the top of the file. Your written solutions can be in
.txt, .docx or .pdf format.
Your code (and all other files) must be submitted as a single zip file to the assignment tab on
Canvas. The zip file should be named using your SID (e.g. cosc 012345678.zip), or both SIDs
if you are submitting as a pair (e.g. cosc 012345678-111222333.zip).
Ensure you have included the following in your zip file:
1. The programs zombiesim.m, montyhall.m, and (COSC1903 only) montecarlo.m;
2. Written responses to Q1.1, Q1.2, Q1.3, Q1.4, Q2.2, Q2.3, and (COSC1903 only) Q3.2, Q3.3,
and Q3.4; and
3. Plots for Q1.2 and (COSC1903 only) Q3.3.
You must submit your zip file by 11:59pm on the due date.
We encourage cooperation between students in completing assignments, and either one or two
students may hand in a single assignment; we recommend you complete the assignment with your
usual COSC lab partner. If you submit this assignment with another student, you certify that you have
both made a fair contribution to it, and are happy to both receive the same mark.
We will NOT accept assignments that are simply copied. Copying the work of another person
without acknowledgement is plagiarism and contrary to University policies. By uploading your sub-
mission to Canvas, you are certifying that you have read and understood the University Academic
Dishonesty and Plagiarism in Coursework Policy.
1
COSC1003/1903 Assignment
Question 1 (15 marks)
Once you know how to run simulations, the obvious thing to do is see what would happen in a zombie
apocalypse.
Conceptually, the zombie apocalypse works as follows. We begin with a 2-dimensional grid of
squares (the “world”). A number of healthy people are randomly distributed throughout the world, as
well as one zombie. Each day, each person either stays still or moves one square up, down, right, or
left. If a zombie spends a day on the same square as a healthy person, the healthy person becomes
infected (i.e., becomes a zombie).
Your program zombiesim.m should have four parts. First, specify the parameters of the simula-
tion. Second, initialise the people who will live (and die) in the simulation. Third, run the simulation
(i.e., at each timestep, update the locations of each person and whether or not they are a zombie).
Finally, plot and output the results.
Parameters
The simulation should be initialized with four parameters: the size of the 2-dimensional world people
inhabit (sidelength, default 40 squares per side), the maximum time to run the simulation for
(maxtime, default 1000 days), the number of people (npeople, default 100), and the rate at which
zombies spontaneously enter remission and return to being healthy humans (remission, default 0).
People
The best way to represent people in a simulation is using a structure array. Create a structure array
to hold each person, with fields for (i) their x coordinate (a random integer between 1 and the size of
the world), (ii) their y coordinate (same), and (iii) whether or not they’re a zombie (no, except for one
person).
To create the structure array, use something like this:
1 people = struct();
2 for p = 1:npeople
3 people(p).x = ...
4 people(p).y = ...
5 people(p).zombie = ...
6 end
Simulation
The “simulation” consists of the calculations that occur inside the loop over time (e.g. for t =
1:maxtime).
First, we need to figure out where people move to. For this, we’ll need an inner loop over people
(e.g. for p = 1:npeople). Create a random variable and map it onto each of the five possible
movements (none, left, right, up, down), then update people(p).x and people(p).y accord-
ingly. If people move outside the world (x;y sidelength), reverse the direction of
their step (e.g., a person at x = 1 who tries to move to x = 0 should move to x = 2 instead).
Second, we need to figure out which people are zombies and which are not zombies. To get the
indices of the zombies, we can use
1 zombies = find([people.zombie]);
2
COSC1003/1903 Assignment
(Note the square brackets – we need this since people.zombie returns a different output for each
element rather than an array.) Once you know the indices of the zombies, there are various ways you
can calculate the indices of the non-zombies.
Third comes the tricky part: we need to handle infections. To do this, we can loop over just the
zombies using e.g. for z = zombies (remember that zombies is already an array, so we can
use it directly in our for loop). Next, we need to loop over all the non-zombies using e.g. for nz =
nonzombies. Finally, for each zombie/non-zombie pair, we need to check whether the non-zombie
is occupying the same square as the zombie (i..e., their x and y coordinates are identical), and if so,
infect the non-zombie (i.e., set people(nz).zombie = 1).
Finally, for each zombie, we need to check if they spontaneously return to being a human: i.e.,
with probability remission, set people(z).zombie = 0.
Plotting
To plot healthy people in blue and zombies in green, include code like this inside the time loop (again
note the square brackets!). This is just an example, so feel free to customise:
1 clf % Clear the current figure
2 hold on % Allow us to plot multiple things simultaneously
3 scatter([people.x], [people.y],'filled') % Plot the locations of all people
4 scatter([people(zombies).x], [people(zombies).y],'filled','g') % Plot the zombies, in green
5 title(sprintf('t=%i zombies=%i', t, length(zombies))) % Put useful information in the title
6 xlim([1,sidelength]) % Set the axis limits to include the entire world
7 ylim([1,sidelength])
8 pause(0.01) % Pause for a short moment so the plot will update
(Note: since the array zombies isn’t updated until the next timestep, including this code immediately
after the new zombies are calculated will technically plot the zombies from the previous step. That’s
fine; consider that it takes a day for someone to become fully infected or recovered.)
If you’ve done everything correctly, the simulation should take about 1 second of actual time for
every 10 days of simulated time, and you should get an output like this:
3
COSC1003/1903 Assignment
Questions
(1.1) What is the average length of time until the last human gets infected? (Run just enough times
to get a decent estimate of the uncertainty, e.g. 10 times.)
(1.2) Calculate the theoretical number of infections per day as a function of zombie prevalence (i.e,
the fraction of people who are zombies). You can do this by calculating (a) the probability that
any given square is occupied by a zombie (you can ignore the fact that there could be more
than one zombie per square), and (b) the total number of healthy people remaining. Next, for
a single run of the simulation, plot the number of zombies as a function of time. (Hint: at
each timestep, simply store the current number of zombies in a vector.) Discuss how the curve
showing number of zombies as a function of time is related to the curve showing the number of
infections as a function of prevalence.
(1.3) If you halve the size of the world (i.e., set sidelength=20), how does this change the time it
takes until the last human is infected? Briefly discuss the real-world implications of this result,
e.g. in terms of how diseases such as tuberculosis spread in crowded hospitals and prisons.
(1.4) Using the original size of the world (sidelength=40), change the spontaneous zombie re-
mission rate from 0 to 0.01 (i.e., every 100 days on average a given zombie will spontaneously
return to being human). Discuss the new dynamics of the model, with reference to steady states
and/or equilibrium points (if you wish, you can also include a graph of the number of zombies
over time to illustrate your point).
4
COSC1003/1903 Assignment
Question 2 (5 marks)
(2.1) In the famous Monty Hall game show problem, contestants pick a door from a set of three doors.
One door hides a car, while the other two have goats. The host, Monty, then opens one of the
doors containing a goat. He then asks if the contestent would like to stay with their initially
selected door, or change to the other unopened door. Then their door is opened. This can be
summarised as
(a) Pick a door where the car is, and pick a door for the contestent
(b) Open a door that was not picked, and contains a goat. If the contestant picked the door
with the car, choose randomly between the two possible doors
(c) Contestent can choose whether to switch or not
(d) See if the contestent has selected the car
Write a function montyhallgame() that performs a single game, and returns 1 if the contes-
tant won (which, following tradition, means choosing the car), and 0 if they lost. The function
should take an argument specifying if the contestant switches doors or not. [Note: you may
use the montyhall.m example covered in lectures for inspiration, but your code should re-
flect your own understanding of the problem. Marks will be deducted for not following the
instructions above, such as not writing a function that performs a single game.]
(2.2) Run your program 10,000 times, and determine how many times the contestant wins if the
contestant never switches, or the contestant always switches. What is the optimal strategy for
the contestant to follow?
(2.3) Change the number of doors from three to 100 (such that there are now 99 goats and one car).
Comment on how the results change and why.
5
COSC1003/1903 Assignment
Question 3 (COSC1903 only) (10 marks)
(3.1) An alternative to traditional methods of numerical integration such as the trapezoidal rule is
Monte Carlo integration. In this method, the estimate of the integral of a function f(x) over an
interval [a;b] on x is given by Z
b
a
fdx (b a)hfi (1)
where
hfi= 1N
NX
i=1
f(xi) (2)
It can be calculated using the following steps:
(a) Choose N random points on the interval [a;b]
(b) Calculate the value of f(x) at each of these points
(c) Use the average value of f(x) in the estimate of the integral
Write a function integratemonte() to implement Monte Carlo integration for a given
number of points (N). Write a program montecarlo.m that calls this function.
(3.2) Run your program to calculate the integral of
f(x) = x3 sinx+x2 + 50 (3)
using 50 points over the range 0 x 5. Compare your answer to the analytic solution.
(3.3) The error term for Monte Carlo integration is given by
E = (b a)rhf2i hfi2N (4)where
hf2i= 1N
NX
i=1f2(xi) (5)
Write a function to calculate the error E. Plot the error for a varying number of points between
1 and 10;000. Use your plot to discuss how the error changes with number of points N. Plot
your results using a log scale and extrapolate (or extend your plotting range) to work out how
many points are needed to calculate the value of the integral to one decimal place.
(3.4) Briefly discuss (with examples) an advantage and disadvantage of Monte Carlo integration in
comparison to a traditional technique such as the trapezoidal rule.