Project Instruction¶
• Please rename this file so that you know which copy you have been working in. Keep a copy safe (especially if you are working in the online Jupyter service). You can download a copy by choosing -File- then -Download as- Notebook from the menu above.
• Complete all of the tasks.
• Make sure your code is readable, organised, and commented appropriately.
Task 1 – Code review¶
This task is to write a code review, not to write python code to solve the problem brief.
A colleague has been asked to write a program to calculate a root of a continuous function using the bisection method as described in the following brief:
Brief¶
The bisection method is a simple numerical technique to find a root of a continuous function in an interval where this function changes sign. According to Intermediate value theorem, the continuous function $f=f(x)$ has at least one root in the interval $[a, b]$, if $f(a)f(b) \leq 0$. Utilising this fact, the following 3-step algorithm will find the root of such function within a desired accuracy:
1. Calculate the midpoint $c = (a+b)/2$ and evaluate the function at this point $f(c)$.
2. If $ (b-a)/2 < \epsilon_1 $ or $|f(c)| < \epsilon_2 $, then return $c$ as the root and stop.
3. Otherwise, depending on the sign of $f(c)$ replace either $a$ or $b$ with $c$, such that $f(a)f(b) < 0$ for the new $a$ and $b$. Then go to step 1.
The criteria introduced in step 2 ensures the difference between the answer and the real root to be less than $\epsilon_1$ or that the function value at the answer is smaller than $\epsilon_2$ (so it can be considered an approximate root). You are supposed to write a function that takes $a$, $b$, $\epsilon_1$ and $\epsilon_2$ as inputs and returns a root using the bisection method.
Note that the bisection method does not guarantee or detect a change of sign in an interval if the two end points have similar signs. Hence, your program should first search for some subinterval where the function changes sign. This can be achieved by halving the intervals consecutively until in one of the smaller subinterval the function changes sign. If no such an interval is found while the size of the smallest interval is still bigger than $\epsilon_1$, the function returns 'None' and prints out "failed to find a root".
Test your code for a continuous function on a given interval.
Your task:¶
You have been asked to write a review of their code. Here is the code they wrote:
In [ ]:
You should write your review here. Things you could choose to discuss:
• Code structure
• Code style
• Does it answer the brief?
• Does it work? If not, could it be fixed?
• Can you explain what it does?
Keep your answer relatively brief (approx. 500 words).
Answer:¶
Add your answer here...
Task 2 - Traffic Modelling¶
Task 2a - update rule¶
For this task, you will work on a model of road traffic using discretised cells, which is an example of a larger class of models in computer science called cellular automaton. Imagine a road that is divided into a number of cells which can contain only one car. For now, let's assume this road is one-way from left to right. We also discretise time into steps. At each step a car moves to its adjacent right cell if it is empty; otherwise, it stays where it is (see the figure below).

We use a periodic boundary condition such that a car that moves off the right-most cell enters into the left-most cell as shown by the red arrow in the schematic above (if you like you can think of this problem as representing a roundabout rather than a straight section of a road).
You should write a function that finds the position of cars in the next step given their positions in the current step. To do this systematically, we can define the numpy array R[i,t] that is 0 if there are no cars in the "i" section of the road, and 1 if a car is present at that section for the time step "t". "i" varies from 1 to $N$ (the total number of cells), and "t" varies from 0 to $T$. The new value R[i,t+1] depends on its old value at time $t$ (i.e. R[i,t]), and also on the old values of the neighbours (R[i-1,t] and R[i+1,t]). You might think about how you would fill out the tables below (on paper - no need to edit the table in your file) to use them to get an explicit form of the update rules (note we use notation here that R[i,t] is the same as $R^t(i)$ below):

To test your update function you can use the following example of a road with 9 cells. If we label the cars by their initial cell position, this demonstrates how they move in the next three steps.

The array R, for this example, becomes
$$ \left(\begin{array}{cccc} 1 & 0 & 1 & 0\\ 0 & 1 & 0 & 1\\ 0 & 0 & 1 & 1\\ 1 & 1 & 1 & 0\\ 1 & 1 & 0 & 1\\ 1 & 0 & 1 & 0\\ 0 & 1 & 0 & 1\\ 0 & 0 & 1 & 0\\ 1 & 1 & 0 & 1\\ \end{array}\right) $$
Throughout this project, let's denote the total number of cells by $N$, the total number of cars by $M$ and the final time step by $T$.
For the first part of this task you should write code to perform each of the following sub-tasks or markdown text to answer any discussion questions:
1) Write a function updateR(R[i,t]), which takes R[i,t] at the current time step as its argument and returns the updated value R[i,t+1] for the next step. Make sure you properly implement the periodic boundary condition in this function. Test your function with the example above. Your function should work for any choice of positive integer $N>3$.
2) Numbering your cells from $1$ to $N=100$, put the cars initially (at $t=0$) in those cells with a prime number index and also the cells in $(40,55]$. Calculate and print out the position of the cars at the final time $T = 400$.
3) Plot the position of the cars from problem 2) in the following three time intervals: $t=[0, 20]$, $[300, 320]$ and $[380, 400]$. Discuss whether you think that the traffic reaches some steady state by the end of the simulation. What is the average velocity of all the cars at the end of the simulation? Note you can either point markers in plot, or use the plotting function imshow to create a checkerboard-like plot.
In [ ]:
In [ ]:
Task 2b – Trajectory of the cars and average velocity¶
Storing R[i,t] at each time step is unnecessarily memory-consuming. Moreover, it is hard to track individual cars through the zeros and ones of R[i,t]. Instead, it is more efficient for many purposes to store the trajectories of each car in an array, Traj[car,t], where each row records the trajectory for a different car, with the value in each successive column recording the location at each successive time step. Our next purpose is to use this array to calculate the average velocities of each car. So it is now better not to use the periodic boundary condition immediately, but to store the value of $N+1$ for the car that moves off the last cell (similarly storing the value of $N+2$ instead of 2 etc). For instance, for the 9-cell road shown in the previous figure Traj is constructed as below:
$$ \left(\begin{array}{ccc} 1 & 2 & 3 & 3\\ 4 & 4 & 4 & 5\\ 5 & 5 & 6 & 7\\ 6 & 7 & 8 & 9\\ 9 & 9 & 10& 11\\ \end{array}\right) $$
To derive Traj[:,t+1] from Traj[:,t], we only need the current column of R and not the entire matrix R[0:N,0:t+1]. Hence, we are going to only keep the vector r[0:N], which is the current state of the road with empty cells represented by 0 and cells with a car by 1 (i.e. r=R[:,t]). Now write a new funtion new_position(r,Traj) that takes the vector r and the current position vector Traj[:,t], and returns the next column Traj[:,t+1].
Using Traj[:,t+1] and modular arithmetic you can update r more easily. To update r from the state at t to t+1 you will need to:
• start a new array for r at t+1 by filling r[:] values as zeros initially
• take each element of Traj[:,t+1]
• find the remainder when you divide that element by $N$, call it j
• set r[j-1]=1.
(Note what is happening here – you are storing the cars in cells 1 up to N in array positions indexed 0 up to N-1).
Write another function new_updateR(Traj[:,t+1]) that takes the trajectory at t+1 and returns a new vector r without using any if statements using this procedure outlined above. (Note here the size of r is $N$ and the size of Traj is $M \times T$).
The speed of car i at time t can be readily calculated as Traj[i,t]-Traj[i,t-1], which is either 0 or 1 cell per step. Write the function ave_vel(Traj,t) that takes the array of trajectories and the time step and returns the average velocity of all cars at time t.
Using the functions you have written; new_position(), new_updateR() and ave_vel(), complete the following list of tasks:
1) Test your function new_position() by using the example above.
2) Put $M = 15$ cars in the first fifteen cells of a road that has $N = 50$ cells. Plot the average speed of all cars as a function of time up to $T=50$. How long does it take to reach a steady average speed?
3) Repeat question 2) for $M = 25$ and $M = 35$ and discuss the results. Is the final average speed the same for all $M$? Which is the largest value of $M$ for which you can reach the maximum average speed?
In [ ]:
In [ ]:
4) Does the average speed of the cars depend on their initial configuration? To answer this question, initally fill $M = 20$ cars in $N=50$ cells in three different ways:
a) Fill the first 20 cells with 20 cars.
b) Use random.sample(range(0,N),M) to randomly distribute the cars (you will need to import random to use this function).
c) Put them in pairs with one cell space between the pairs. That is, fill the cells 1,2,4,5,7,8,10,11,…
Plot the average speed of the cars as a function of time up to $T=50$ for each case and discuss the results.
5) Repeat the previous question for $M = 30$ cars. Does the average speed reach the same limit?
In [ ]:
In [ ]:
6) How does the number of cars affect the final average speed? In other words, we want to see how traffic fluidity depends on how busy the road is. Use scipy.stats.bernoulli.rvs(p, size=N) to fill each cell in the inital configuration (see workshop 5 if you are not familiar with this function). This means the chance of having a car in a cell is p. Therefore, if N is large enough, you expect to have N*p cars on the road. Set $N=800$, vary p from $0.2$ to $0.8$, and plot the average speed of cars at T=200 as a function of p. Comment on the results. How does the number of cars affect their final average speed?
7) Repeat the last question for $N=50$ and $N = 2000$. Do the results change? Why?
In [ ]:
In [ ]:
Task 2c – Fast and slow cars¶
Now imagine there are two types of cars on our road: fast cars who want to move with the speed of 2 cells per step if there are enough spaces in front of them, and the slow cars who always move with the speed of 1 cell per step as in the previous task. If there are less than two free cells in front of the fast cars, they either move one cell per step or don’t move at all (like slow cars).
1) Numbering your cells from $1$ to $N=100$, put all fast cars initially in the prime indexed cells and the cells in the interval $(60,85]$. Now plot the position of the cars in the time interval $t=[0, 60]$. Do you think a part of road will stay congested forever? In which direction is the congested part of the traffic flow moving?
(You will find it helpful to construct a new function to help new_position2(r,Traj,Vel), where the new vector Vel stores the velocity of each car).
In [ ]:
2) Now put fast cars initially in the prime indexed cells only like those in the last question, but this time make the last car (in cell 97) slow (i.e. moving only one cell per step). How does this car affect the whole traffic flow?
3) By running the simulation up to $T = 100$, calculate the average velocity of the cars at this final time.
In [ ]:
In [ ]: