Exam practice

A model is a simplified description of a real system: what it keeps, its rules and its assumptions. Example: the SIR model, three groups of people and two rules for how they move between them.
A simulation is a model run forward in time by a program, step by step, to see what the system does. Example: the program that steps the SIR model forward a day at a time for Bhaktapur district and prints how many are ill on each 30th day.

Leaves out (any three): traffic jams; road works; breakdowns; the weather; the driver; how full the bus is; police checks.
Assumes (any three): the road is open all the way; the speed on each stretch is the same every day; the meal stops always take the same time; the bus doesn’t break down.
In the monsoon the assumption that the road is open is most likely to fail: landslides often block the highway for hours. A good model for the monsoon would include the chance of a blockage and how long it lasts, and give a range of times, not one.

The update rule moves the ball for a whole time step at the speed it had at the start of the step. In reality the speed changes all through the step (the ball slows as it rises), so each step is a little wrong, and the errors add up. With a smaller step the speed changes less within each step, so each step’s error is smaller; for Euler’s method the total error shrinks in proportion to the step (0.5 s: 15.8 m off; 0.01 s: 0.2 m).
It isn’t free because a smaller step needs more steps for the same flight, so more work: 10 times smaller step, 10 times as many steps. Very small steps also let rounding errors build up over millions of steps.

class Ball:    """A ball in the air, moved one time step at a time (no air resistance)."""     def __init__(self, vx, vy):        self.x, self.y = 0.0, 0.0        self.vx, self.vy = vx, vy     def step(self, dt):        self.x = self.x + self.vx * dt        self.y = self.y + self.vy * dt        self.vy = self.vy - 9.8 * dt  ball = Ball(0, 10)t = 0.0while ball.y >= 0:    ball.step(0.001)    t = t + 0.001print(f"In the air for about {t:.2f} s")

In the air for about 2.04 s

The exact answer is 2 × 10 ÷ 9.8 = 2.04 s. A class keeps the ball’s state and its update rule together; a simulation with many balls (or planets, or people) makes one object for each and calls step on all of them.

She is wrong: the drift isn’t from the computer’s rounding, which is tiny. It comes from Euler’s method itself. Each step moves the planet in a straight line along the way it was moving at the start of the step, and a straight line touching a circle lies outside it, so every step lands a little outside the orbit. All these errors push the same way, outwards, so they add up over thousands of steps.
To fix it, change the order of the update: change the velocity first, then move with the new velocity (the Euler–Cromer method), whose errors balance round each orbit; or use a much smaller step (which slows the drift but doesn’t stop it), or a better method such as Runge–Kutta.

On the way up, most people are still susceptible, so each infected person infects more than one other (R₀ × the share susceptible is above 1): new cases outnumber recoveries, and the number ill grows faster and faster.
As people catch it and recover, the susceptible become fewer. At the peak, each case is making exactly one new case on average: the new cases in a day just equal the recoveries in a day, so the number ill stops rising.
After the peak, the susceptible are too few: each case makes fewer than one new case, recoveries outnumber new cases, and the number ill falls until the epidemic ends, with some people never having caught it.
The number ill changes by (new cases − recoveries) each day. It is highest at the moment it stops rising and starts falling, which is when that difference is zero.

1 − 1 ÷ R₀:

• R₀ = 1.5: 1 − 0.67 = 0.33, about 33%.
• R₀ = 3: 1 − 0.33 = 0.67, about 67%.
• R₀ = 6: 1 − 0.17 = 0.83, about 83%.

With R₀ below 1, each case makes less than one new case on average even when everyone can catch it, so the number ill shrinks from the start: 10 cases make fewer than 10, which make fewer still, and it dies out.

1. When there are many prey and few predators, the prey are born faster than they are eaten (births × prey is bigger than eaten × prey × predators), so the prey grow.
2. With many prey, the predators eat well: fed × prey × predators is bigger than deaths × predators, so the predators grow too.
3. As the predators grow, they eat more prey: eaten × prey × predators becomes bigger than the births, and the prey start to fall. That is the prey’s peak.
4. The predators keep growing for a while after that, because there are still enough prey to feed them; only when the prey have fallen below deaths ÷ fed do the predators’ deaths outnumber their young. That is the predators’ peak, after the prey’s.
5. Then both fall, the prey recover once the predators are few, and the cycle begins again.

import random random.seed(5)rolls = 100000sevens = 0for i in range(rolls):    if random.randint(1, 6) + random.randint(1, 6) == 7:        sevens += 1print(f"Chance of 7: {sevens / rolls:.4f}")print(f"Exact: {6 / 36:.4f}")

Chance of 7: 0.1663Exact: 0.1667

Exactly: of the 36 equally likely pairs, 6 add up to 7 (1+6, 2+5, 3+4, 4+3, 5+2, 6+1), so the chance is 6 ÷ 36 = 1/6 ≈ 0.1667. With 100,000 rolls the estimate is usually within about 0.002 of it.

A Monte Carlo answer depends on the random numbers used, so another run gives a different answer. With few trials the answers can differ a lot (the birthday chance from 10 classes might be 60% or 100%), so one number alone could mislead someone into thinking it is exact.
To find how much it varies: run the whole simulation several times (or with different seeds) and report the range or the spread of the answers; and remember that the spread shrinks with 1 ÷ √trials, so to halve it, use four times the trials.

Reasons:

• Fewer births: women in Nepal had far fewer children by the 2010s (2.1 each by 2022, more than 6 in the 1970s). The model assumed the growth rate would stay the same, but the share of births fell.
• People abroad: about 2.19 million Nepalis were living abroad in 2021, almost three times as many as in 2001, and the census counts only the people living in Nepal. The model assumed the same share of people would leave each year.

Both show the model’s one assumption, a steady growth rate, was false: the things that make up growth (births, deaths, people leaving and coming back) changed.
To test it, the planner should have fitted it on the earlier censuses only (say 1961 to 1981) and checked its prediction against 1991 and 2001, data it hadn’t seen; and given a range, not one number, for example by fitting rates over different periods. Better still, a model with births, deaths and migration as separate parts, each fitted to its own data.

Any four, such as:

1. What does the model assume? Doubling every week until the end of the year means the cases grow without limit; does it allow for people becoming immune, for the mosquito season ending as the weather cools, for people protecting themselves?
2. Was it tested on data it wasn’t fitted to? Did its forecasts for past weeks match the real cases?
3. How sure is it? Does it give a range, and what happens to the forecast if its parameters are a little different?
4. How far ahead? A forecast for next week can be trusted far more than one months ahead.
5. Where does the data come from? Reported cases depend on how many people are tested; a rise could be partly more testing.
6. Who made it, and is it published so others can check it?

(i) The sensor is the float (or the level switch) in the tank, which measures the water level; the controller is the box with the rule “on below 250 liters, off at 950”; the actuator is the pump. The set point is a full tank; the gap between switching on at 250 liters and off at 950 is the dead band.

(ii) and (iii):

def use(minute):    """Liters a minute the house uses, by the time of day."""    hour = minute // 60    if 6 <= hour < 8 or 18 <= hour < 20:        return 8                     # mornings and evenings: cooking, washing    if 8 <= hour < 22:        return 2    return 0                         # night  def day(low, high):    """One day with the pump on below low liters and off at high. Gives back how many times it starts."""    level = 600.0                    # liters in the tank at midnight    pump = False    starts = 0    for minute in range(24 * 60):        if level < low and not pump:            pump = True            starts += 1        elif level >= high:            pump = False        level = level + (20 if pump else 0) - use(minute)   # the pump gives 20 liters a minute    return starts  for low, high in [(250, 950), (940, 950)]:    print(f"on below {low} L, off at {high} L: the pump starts {day(low, high)} times")

on below 250 L, off at 950 L: the pump starts 3 timeson below 940 L, off at 950 L: the pump starts 120 times

With a dead band of only 10 liters, the pump starts 120 times in a day: in the daytime the house uses those 10 liters in about 5 minutes, so the pump starts every few minutes, runs for half a minute and stops. Each start draws a large current and wears the motor and the switch, so the pump would fail much sooner. The wide band lets the pump run a few long times a day instead, which is the same reason a room thermostat has a dead band. The tank is a little less full on average, which doesn’t matter while it never runs dry.