#### 您當前位置：首頁 >> Python編程Python編程

###### 日期：2019-11-11 11:10

Exercises for IN1900

October 14, 2019

Preface

This document contains a number of programming exercises made for the course

IN1900. The chapter numbers and titles correspond to the chapters of the book

“A primer on Scientific Programming with Python” by Hans Petter Langtangen.

The exercises are meant to be a supplement to the exercise collection in the book,

and most are motivated by applications in science and applied mathematics.

The exercise collection is used for the first time in 2018, and there may be typos

and small errors. If you find any errors, or have other comments or questions

1

Chapter 1

Computing with Formulas

Problem 1.1. Throw a ball

When throwing a ball in the air, the position of the ball can be calculated using

the acceleration of the ball. When neglecting air resistance, the acceleration will

be the negative of the gravitational constant, ?g. The height of the ball relative

to its starting point is

y(t) = v0t ?12gt2,

where v0 is the initial velocity of the ball and t is the time after the throw. The

ball reaches its maximum height at time

tmax =v0g.

Write a program computing the maximum height of the ball, that is y(tmax),

when v0 = 8.2m/s and g = 9.81m/s2. Print the result.

Filename: ball.py

Problem 1.2. Population growth

The growth of a population can often be described by a logistic function

N(t) = B1 + Ce?kt ,

where B is the carrying capacity of the species in the environment, i.e., the

maximum size of the population that the environment can sustain indefinitely.

The constant k tells us something about how fast the population grows, while C

is given by the initial conditions. Let us consider a bacterial colony where we

take the carrying capacity to be B = 50000 and k = 0.2h?1. If the population is

5000 at t = 0, find C and write a code that finds the number of bacteria in the

colony after 24 hours.

Filename: population.py

2

Problem 1.3. Solve the quadratic equation

ax2 + bx + c = 0,

the two roots are.

Make a program evaluating the roots of.

Print out both roots with two decimals.

Filename: find_roots.py

Problem 1.4. Forces in the hydrogen atom

There are two kinds of forces acting between the proton and the electron in the

hydrogen atom; Coulomb force and gravitational force. The Coulomb force can

be expressed as,

where ke is Coulomb’s constant, e is the elementary charge, and r is the distance

between the proton and the electron.

The gravitational force can be expressed as

FG = Gmpmer2,

where G is the gravitational constant, mp is the mass of the proton, me is the

mass of the electron, and r is the distance between the particles.

We can use these expressions for FC and FG to illustrate the difference in

strength of these two forces, i.e., the electromagnetic force and gravitational force.

Use the values ke = 9.0·109 Nm2C

?2

, e = 1.6·10?19 C, G = 6.7·10?11 Nkg?2m2

,

mp = 1.7·10?27 kg and me = 9.1·10?31 kg.You can take the distance between the

proton and electron to be approximately the Bohr radius r = a0 = 5.3 · 10?11 m.

Make a program that computes both the Coulomb force and the gravitational

force between the proton and the electron. Write out the forces in scientific

notation with one decimal in units of Newton (N = kgm/s

2

). Also print the

ratio between the two forces.

Filename: hydrogen.py

3

Chapter 2

Loops and Lists

Problem 2.1. Multiply by five

Write a code printing out 5 · 1, 5 · 2, ..., 5 · 10, using either a for or a while loop.

Filename: multiplication.py

Problem 2.2. Multiplication table

Write a new code based on the one from Problem 2.1. This code should print

the whole miltiplication table from 1 · 1 to 10 · 10.

Hint: You may want to consider using one loop inside another.

Filename: mult_table.py

Problem 2.3. Stirling’s approximation

Stirling’s approximation can be written ln(x!) ≈ x ln x ? x. This is a good

approximation for large x. Write out a nicely formatted table of integer x values,

the actual value of ln(x!), and Stirling’s approximation to ln(x!).

Filename: stirling.py

Problem 2.4. Errors in summation

The program has three errors and therefore does not work. Find the three errors

and write a correct program. Put comments in your program to indicate what

the mistakes were.

Hint: There are two basic ways to find errors in a program:

1. read the program carefully and think about the consequences of each

statement,

2. print out inter mediate results and compare with hand calculations.

4

1. First, try method 1 and find as many errors as you can. Thereafter, try method

2 for M D 3 and compare the evolution of s with your own hand calculations.

Lastly, write a similar code evaluating the same sum using a while loop.

Check that the two loops compute the same answer.

Filename: sum_for.py

5

Problem 2.5. Binomial coefficient

The binomial coefficient is indexed by two integers n and k and is written . Compute the same value using Eq. (2.1) and check that the results are

correct.

Hint: The Q

sign is a product sign.

When checking the result you will need math.factorial.

Filename: binomial.py

Problem 2.6. Table showing population growth

Consider again the bacterial colony from Problem 1.2. Let us study the number

of individuals for n + 1 uniformly spaced t values throughout the interval [0, 48].

Set n = 12. First store the t and N values in two lists t and N. Thereafter, write

out a nicely formatted table of t and N values by traversing the two lists with a

(separate) for loop.

Filename: population_table.py

Problem 2.7. Nested list

a) Compute two lists of t and N values as explained in Problem 2.6. Store the

two lists in a new nested list tN1 such that tN1[0] is the list containing t-values

and tN[1] correspond to the list containing N-values. Write out a table with t

and N values in two columns by looping over the data in the tN1 list. Each t

and N value should be written in the table as integers.

b) Make a nested list tN2 where tN2[i] contains the i-th element of both the

t-list and the N-list. Loop over the tN2 list and write out the t and N values in

the table as integers.

Filename: population_table2.py

Problem 2.8. Calculate Cesaro mean

Let (an)∞n=1 be a sequence of numbers, sk =Pkn=0 an = a0 + . . . , +ak,

is called a Catalan number. Compute and print the first 10 Catalan numbers.

Filename: catalan.py

Problem 2.10. Molar Mass of Alkanes

Alkanes are saturated hydrocarbons with the chemical formula CnH2n+2. If

there are n Carbon atoms in the alkane, there will be m = 2n + 2 Hydrogen

atoms. The molar mass of the hydrocarbon is MCnHm = nMC + mMH, where

MC is the molar mass of Carbon and MH is the molar mass of Hydrogen.

Use a for-loop or a while-loop to compute and print out the molar mass

of the alkanes with two through nine Carbon atoms (n ∈ [2, 9]). The output

should specify the chemical formula of the alkane as well as the molar mass. An

example on how the formatted output should look like for n = 2 is given below.

M(C2H6) = 30.069 g/mol

You can set the molar masses of the atoms to be MC = 12.011 g/mol and

MH = 1.0079 g/mol

Filename: alkane.py

Problem 2.11. Matrix elements

This exercise involves no programming.

The answers should be written in a text file called matrix.txt

Consider a two dimensional 3 × 3 matrix

In Python, the matrix A can be represented as a nested list A, either as a

list of rows or a list of columns. Find the indices i, j of the Python list A such

that A[i][j] = a11 and the indices k, l such that A[k][l] = a32 for the two

following cases:

a) When A is represented as a list of rows. This means that A contains three

lists, where each list corresponds to a row in A.

b) When A is represented as a list of columns. This means that each element

in A contains a list with the elements of a column in A.

Filename: matrix.txt

7

Chapter 3

Functions and Branching

Problem 3.1. Implement a function for population growth

Consider again the function

N(t, k, B, C) = B

1 + Ce?kt .

Implement N as a python function population(t, k, B, C) that returns the

number of individuals in a population after a time t.

Write out a nicely formatted table of t and N values for the time interval

t ∈ [0, 48] using the values from Problem 2.6.

Filename: pop_func.py

Problem 3.2. Sum of integers

We consider the sum Pn

i=1 i = 1 + 2 + · · · + n of positive integers up to n. It

can be shown that the sum is equal to n(n+1)

2

.

a) Write a function sumint(n) that returns the sum of all positive integers

up to n.

b) Write a function implementing n(n+1)

2

.

c) Write test functions for both a) and b) testing for specific known values.

Filename: sumint.py

Problem 3.3. Implement the factorial

a) The factorial can be implemented by a so called recursive function call. Use

a recursive function call to implement a function myfactorial(n) that returns

n!.

b) Write a test function where you call the myfactorial function and check

the value of the returned object for one value of n using math.factorial.

Filename: factorial.py

8

Problem 3.4. Half-wave rectifier

In a half-wave rectifier the positive part of a signal passes, while the negative

part is blocked. Thus, for a signal passing through a half-wave rectifier, the

negative values are set to zero. Let us look at a sine signal that has passed

through a half-wave rectifier:

f(x) = (

sin x if sin x > 0

0 if sin x ≤ 0.

Implement f(x) as a Python function f(x) and make a test function for testing

the implementation of f(x) in both cases.

Filename: half_wave.py

Problem 3.5. Compute the area of an arbitrary triangle An arbitrary triangle

can be described by the coordinates of its three vertices: (x1, y1),(x2, y2),(x3, y3),

numbered in a counterclockwise direction. The area of the triangle is given by

the formula

Write a function triangle_area(vertices) that returns the area of a triangle

whose vertices are specified by the argument vertices, which is a nested list of

the vertex coordinates. Make sure your implementation passes the following test

function, which also illustrates how the triangle_area function works:

def test_triangle_area():

"""

Verify the area of a triangle with vertices

(0,0), (1,0), and (0,2).

"""

v1 = (0,0); v2 = (1,0); v3 = (0,2)

vertices = [v1, v2, v3]

expected = 1

computed = triangle_area(vertices)

tol = 1E-14

success = abs(expected - computed) < tol

msg = f"computed area={computed} != {expected}(expected)"

assert success, msg

Filename: triangle_area.py

Problem 3.6. Primality checker

Recall that a prime number is a number greater than 1 that has exactly 2 divisors.

Said differently, a number greater than one is a prime if it is divisible by only

itself and one. A number that is not prime is called composite. Every number n

can be written as a unique product of primes (e.g. 12 = 2 · 2 · 3), this is called

the prime factorization of n.

a) Make a function that takes a number n, and returns true if it’s prime, and

false if it’s not. Use the program to find all prime numbers up to 100.

p

Hint: You will only need to check divisibility for numbers up to and including

(n), because any greater divisor will imply that there is a divisor less than this.

9

b) Make a function that instead finds the prime factorization of the input

number. It should print “prime” and return nothing if the number is prime,

and both print and return the factorization if it’s composite. Find the prime

factorization of 5525612.

c) Make test functions for the two functions above where you check for small

values of n.

d) Compare the runtime of the two functions with the number 33425626272.

Is the difference big? If so, why do you think one is faster than the other? The

following code returns the mean time it takes for your program to run once:

import timeit

timeit.timeit(’your_func(args)’, \

’from __main__ import your_func’,number=1)

Filename: prime.py

Problem 3.7. Eulers totient function

Two numbers n and m are called relatively prime if they have no common divisors

except for 1. That is, no number greater than one should divide both numbers

with no residue.

a) Make a function that takes two numbers and returns true if they’re relatively

prime and false if they’re not.

b) Euler’s totient function is defined as

φ(d) = #{Numbers less than d which are relatively prime to d}.

Implement Eulers totient function and print φ(d) for d = 10, 50, 100, 200.

c) Make a test function for both a) and b).

Filename: euler.py

Problem 3.8. Simple Statistical Functions

In this problem you will implement two statistical functions and test them by

comparing the results with statistical functions from the numpy module. We will

trust that the functions from the numpy module are correct, and will use them

as benchmark values in the test functions. When you import the numpy module

you should follow the convention of renaming it np, as shown below.

import numpy as np

a) The mean of a set of numbers x1, x2, x3, ..., xN is defined as

where N is the size of the set. Implement a function mean(x_list) that returns

the mean value of a list of numbers.

10

b) Make a test function test_mean() that tests the function from a). Compare

the returned value with the result from numpy.mean. (Such that

expected =np.mean(x_test_values)).

c) The standard deviation of a set of numbers x1, x2, x3, ..., xN is defined as.

Implement a function standard_deviation(x_list) which returns the standard

deviation of a list of numbers. The mean value of the list will be necessary

to calculate the relative deviation. Obtain the mean value inside the

standard_deviation function by calling the function you made in a).

d) Make a test function test_standard_deviation() that tests the function

from c). Compare the returned value with the result from numpy.std. (Such

that expected =np.std(x_test_values)).

You may use the list below as an example for your test functions.

x_test_values = [0.699, 0.703, 0.698, 0.688, 0.701]

Filename: stat.py

Problem 3.9. Münchhausen Numbers

A Münchhausen number is a number such that the sum of every digit to the

power of itself equals the original number. E.g. 1

1 = 1 is a Münchhausen number,

and 5

5 + 33 + 22 = 3156 6= 532, so 532 is not.

Make a function that checks if a number is Münchhausen. Find a Münchhausen

number different from one.

Hint: There is only one such number different from 1 and also under one

million

Filename: m_numbers.py

11

Chapter 4

User Input and Error

Handling

Problem 4.1. Quadratic with user input

Consider the usual formula for computing solutions to the quadratic equation

ax2 + bx + c = 0 given by.

Write a program that asks the user for values ( a =, b = and c = ) to get values

for a, b, and c through the users keyboard. Use input (or raw_input if you are

using Python 2). Print the solutions.

Problem 4.2. Quadratic with command line

Modify the program from 4.1 such that a, b and c are read from the command

line.

Extend the program from 4.2 with exception handling such that missing command

line arguments are detected. In the except IndexError block, use input (or

raw_input if you are using Python 2) to ask the user for missing input data.

Problem 4.4. Quadratic with raising Error

In this exercise, use the sqrt function imported from math.

Consider the program from Problem 4.1. Not all inputs yield real solutions.

Modify the program such that it raises ValueError if the values for a, b and c

yield complex roots. (That is if b

2 ?4ac < 0). Provide a suitable Error-message.

Test that your program prints out real roots and that it raises ValueError when

the roots are complex. (An example of values that provide complex roots could

be a = 1, b = 1, c = 1, while a = 1, b = 0, c = ?1 provide real roots).

12

Problem 4.5. Estimating harmonic series

Let f(x) be the function

Write a program that approximates f(x) (that is, evaluates fN (x) = PN

with values of x and N given as command line arguments. Run the program for

x = 0.9, x = 1, and N = 10000. Print the results.

Remark. For x = 1 this is known as the harmonic series. Despite the low values

for large N, the series does not converge, but diverges very slowly. Try to run

the program for different values of N to see how big you can get the value of

f(1).

Filename: harmonic.py

Problem 4.6. Estimating harmonic series extended

Using the program from Problem 4.5, consider the following values for x and N

in a text file

x: 0.9 1

N: 500 1000 10 100 50000 10000 5000

a) Write a function to read a file containing information in the above format

that returns two lists containing the values of x and N.

b) Write a test function for a) that generates a file in the given format and

checks that the values returned by the function is correct.

c) Use the program from Problem 4.5 to evaluate fN (x) for the different values

of x and N. Create a function that writes the information to a file in a table

format with the first column containing the values of N in increasing order, and

the second and third the values of fN (x) at 0.9 and 1 respectively.

Filename: harmonic_table.py

Isotopes of a chemical element in its ground state have the same number of

protons but differ in the number of neutrons. The weight of isotopes of the same

chemical element will therefore be different.

The molar mass, M, of a chemical element, can be calculated by summing

over all its isotopes M =Pi miwi, where mi

is the weight of the i-th isotope

and wi the corresponding natural abundance.

The file Oxygen.txt, which is given below, contains the information on

Oxygen’s isotopes (16O,

17O and 18O).

Isotope weight [g/mol] Natural abundance

(16)O 15.99491 0.99759

(17)O 16.99913 0.00037

(18)O 17.99916 0.00204

13

Write a script in Python to read the file Oxygen.txt and extract the weights

and the natural abundance of all the isotopes of Oxygen. Use these to calculate

the molar mass of Oxygen. Print out the result with four decimals and provide

the correct units.

Problem 4.8. A result on prime numbers

A famous result concerning prime numbers states that the number of primes

below a natural number n, denoted π(n), is approximately given by

tends to 1 as n → ∞. The following

table contains the exact values of π(n) for some values of n.

n: 10**20 10**4 10**2 10**1 10**12 10**4 10**6 10**15

pi(n): 2220819602560918840 1229 25 4 37607912018 168

78498 29844570422669

a) Write a function that reads the file given above and returns two tuples

containing sorted values of n and π(n). It is important that the correspondence

in the orderings are correct, that is, the same as in the table above.

b) Write a test function that generates a file with the format above and tests

that the returned values are correct. It should test that the order of the elements

are in correspondence as in the file.

Hint: The == operator on tuples will take the order into account. The same

operator on lists will not.

c) Create a function that writes the values of n and p(n) to a file in a table

format in increasing order with the values of n in the first column and the

corresponding values of p(n) in the second column.

Bonus problem There are better approximations to π(n), for example the

function

Approximate the integral for different values of n and modify the program to

write these into a third column.

Hint: Implement an algorithm for approximating the integral (e.g. the

trapezoidal rule) and compute the difference as before.

Filename: primes.py

14

Problem 4.9. Conversion from other bases

Recall that a binary number is a sequence of zeros and ones which converted to

the decimal system becomes P

i

2

i where i is a term in the sequence containing

a 1 (e.g. 100101 = 25 + 22 + 20 = 37).

a) Write a function that takes a binary number and converts it to a decimal

number. If the argument is not a binary number, a message should be printed

and nothing returned.

Hint: Let the number in the argument be of type string to avoid problems

with numbers starting with a zero.

b) Let the binary number from a) be taken as a command line argument.

Use exceptions (IndexError) to handle missing input. Print the conversion of

100111101.

c) Extend the program with a function to also handle numbers written in base

3.

Hint: An example of a ternary number(a number in base 3) converted to a

decimal number: 1201 = 1 · 3

Filename: base_conversion.py

Problem 4.10. Read temperatures from two files

We consider data sets from the Norwegian Meteorological Institute, containing

daily mean temperatures of any month of any year at Blindern (Oslo).[Ins19]

Each file looks typically like this:

Year: 1997. Month: April. Location: Blindern(Oslo)

9.0 12.3 15.8 13.4 11.0 16.2 13.3

12.9 14.0 14.1 12.0 17.3 15.5 15.4

...

The observations are given chronologically, and the temperatures are given in

degrees Celsius. There are no empty lines in the bottom of the file.

a) Write a function extract_data(filename) that reads any such file and

returns a list of the temperatures from the given month.

b) In the two files temp_oct_1945.dat and temp_oct_2014.dat you will

find observations of daily mean temperatures in October 1945 and October

2014, respectively. Store the temperatures in two lists oct_1945 and oct_2014.

Calculate the average, maximum and minimum value of the temperatures of both

months, and print the results. You may use the numpy.mean(), numpy.max()

and numpy.min() methods.

c) Write a function write_formatting() that takes at least filename, list1

and list2 as parameters, and creates a new file with two nicely formatted

columns containing the temperatures of the given months (you can assume that

the months have equal lengths). Finally, call the function such that the file

temp_formatted.txt is created, using the lists oct_1945 and oct_2014.

15

Chapter 5

Array Computing and

Curve Plotting

Problem 5.1. Fill arrays; loop version

We study the function

f(x) = ln(x).

We want to fill two arrays x and y with x and f(x) values, respectively. Use 101

uniformly spaced x values in the interval [1, 10]. Create empty x and y arrays

and compute each element in x and y with a for loop.

Filename: fill_log_arrays_loop.py

Problem 5.2. Fill arrays; vectorized version

Vectorize the code in Problem 5.1 by creating the x values using the linspace

function from the numpy package and evaluating f(x) with an array argument.

Since the calculation should be vectorized, you may not use any form of loop.

Filename: fill_log_arrays_vectorized.py

Problem 5.3. Plot the population growth

Again, we’re considering a population undergoing logistic growth. The number

of individuals in the population is given by

N(t, k, B, C) = B

1 + Ce?kt .

Plot this function for t ∈ [0, 48] with a carrying capacity B = 50000, C = 9 from

the initial condition that we have 5000 individuals at t = 0 and a steepeness of

k = 0.2.

Filename: population_plot.py

Problem 5.4. Oscillating spring

A rock of mass m is hung from a spring, and pulled down a length A. When

released, the rock will oscillate up and down with a vertical position given by

y(t) = Ae?γt cos r

Here, y is the vertical position of the rock, k is the spring constant, and γ is

a friction coefficient representing air resistance. Set k = 4 kg s?2 and γ = 0.15

s

?1

, m = 9 kg, and A = 0.3 m.

16

a) Create arrays t_array and y_array of size 101, both initially filled with

zeros. Use a for loop to fill them with time values in the range from 0 to 25

seconds, and the corresponding y(t) values.

b) Vectorize your program by using the NumPy’s linspace function to generate

the t_array, and send it into a function y(t) to generate the y_array.

Your program should now be free of for loops. paragraphc) Plot the position

of the rock against time in the given time interval. Use the arrays from both

exercise a) and b), and confirm that they give the same result. Put the correct

units on both axes.

Filename: oscillating_spring.py

Problem 5.5. Plot Stirling’s approximation

Stirling’s approximation is

ln(x!) ≈ x ln x ? x.

a) Make two functions stirling(x) and exact(x), returning Stirling’s approximation

and the exact value of ln(x!), respectively. Plot both the approximation

and the exact curve in the same figure.

Hint: To implement a vectorized version of the exact function, you can use

scipy.special.gamma(x). This function is a “generalized factorial” which can find

the “factorial” of float numbers. It works such that n! = gamma(n + 1). You

can also just consider integer values and plot the value of ln(x!) for each integer

x in the interval you’re considering. Keep in mind that math.factorial is not

vectorized.

b) Use a while loop and find the minimal value of x for the relative error to

be less than 0.1%.

Hint: Relative error is given as (a ? a?)/a, where a is the exact value and a?

is the approximation. Also, do not start with x smaller than or equal to 1, why?

Filename: stirling_plot.py

Problem 5.6. Plotting roots of a complex number

The n’th roots of a complex number z = reiθ can be found by,

for k = 0, 1, ..., n ? 1. The roots can be rewritten to separate the real component

xk and the imaginary component yk, such that ωk = xk + iyk. Through the

relation between the exponential function and the sine functions we get,

for k = 0, 1, ..., n ? 1.

17

a) Write a function that takes the angle θ, the radius r and the degree n of

the roots as parameters. The function should calculate and return all of the n’th

roots of a complex number reiθ, as two lists or arrays corresponding to the real

part x = x0, x1, ..., xn?1 and the complex part y = y0, y1, ..., yn?1 of the roots.

An example of a function call on the function you will write is given below.

x, y = roots(r, theta, n)b) Consider the complex number z = 10?4e

i2π

. Use the function from a) to

get all the roots of order n = 6, n = 12 and n = 24. Plot the roots as points,

and plot all the three orders of roots in the same plot. Label the different orders

of roots. And example of code for plotting the roots of order n = 6 is given

below.

plt.plot(x_n_6, y_n_6, "o", label="n = 6")

Filename: roots.py

18

Problem 5.7. Fermi-Dirac distribution

The Fermi-Dirac distribution says something about the probability of an energy

state being occupied by a particle, or more precisely a fermion, e.g. an electron.

It is a function of energy and temperature given by

f(E, T) = 1

1 + e

(E?μ)/kT , (5.1)

where E is energy, T is temperature, k is Boltzmann’s constant and μ is the

so-called chemical potential. Use k = 8.6 · 10?5

eVK?1 and μ = 4.74eV and

make a program that visualizes the Fermi-Dirac distribution on the interval

E ∈ [0, 10]eV when T = 0.1K. (eV is a unit of energy, 1eV = 1.6 · 10?19J.)

Filename: Fermi_Dirac.py

Problem 5.8. Animate the temperature dependence of the FermiDirac

distribution

Make an animation of the Fermi-Dirac distribution f(E, T) from Problem 5.7

We’re interested in studying how the distribution changes when we raise the

temperature. Plot f as a function of E on [0, 10] for a set of temperatures

T ∈ [0.1, 3 · 104

]. Also make an animated GIF file. Remember to label your axes

and include a legend to show the value of the temperature.

Hint: A suitable resolution can be 1000 intervals (1001 points) along the E

axis, 60 intervals (61 points) in temperature, and 6 frames per second in the

animated GIF file. Use the recipe in Section 5.3.4 and remember to remove the

family of old plot files in the beginning of the program.

Filename: Fermi_Dirac_movie.py

Problem 5.9. Bump functions

Consider the function

f(x) = (

ke? 1

1?x2 ?1 < x < 1

0 otherwise.

a) Plot the function with k = 1 on the interval ?2 ≤ x ≤ 2 by implementing a

b) Animate the function on the same interval as above when k decreases from

1 to 0.

Filename: bump.py

Problem 5.10. Band structure of solids

Electrons in solids are waves. These waves have different wave lengths λ. Often,

waves are characterised by their wave number k = 2π/λ, and the wave number

is associated with the energy of the electron. The energies of electrons in solids

have a band structure, i.e., there are different bands of energies separated by a

band gap.

The file bands.txt contains k-values and corresponding energies for the

three first bands of a solid. Have your program read the values for k and the

energies and plot the energy bands as functions of k in the same figure. You will

see that some energies can never be obtained by electrons in the solids. These

areas of non-allowed energies are called the band gaps.

Filename: band_structure.py

19

Problem 5.11. Half-wave rectifier vectorized

In Problem 3.4, we implemented a function illustrating a sine signal after it had

passed through a half-wave rectifier. Vectorize this function and plot f(x) for

x ∈ [0, 10π].

Hint: The numpy.where(condition, x1, x2) function returns an array

of the same length as condition, whose element number i equals x1[i] if

condition is True, and x2[i] otherwise.

Filename: half_wave_vec.py

Problem 5.12. Singularity plot

In this problem we consider the function. Create arrays of r and θ values on the unit

circle centered at the origin with n uniformly spaced values. Fix axes between

-0.5 and 0.5 for x and y and visualize the function for n = 10, 50, 100, 500. You

can use the following to generate the correct values for r and θ:

theta = np.linspace(0,2*np.pi,100)

r = np.linspace(0.01,1,100)

r, theta = np.meshgrid(r,theta)

Remark. If we had an ideal computer that could calculate every value in an

interval and plot it, then the image we have plotted would touch every single

value in the plane, except for at most one! In our program we have 0.01 < r < 1.

The remarkable thing is that the same is true if we replace the inequality with

0 < r <  for any  > 0. Not only that, but all those points are hit an infinite

number of times!

Filename: ess_sing.py

Problem 5.13. Approximate |x|

The absolute value f(x) = |x| can be written as a sum.

Write a program that calculates the first N terms for N = 1, 2, 3, 4 and plots it

against the exact function. Let the x-axis be [?π, π] with a suitable y-axis.

Filename: approx_abs.py

Problem 5.14. Plotting graphs

A graph is a collection of lines and points in the plane such that each line

connects two points. In this exercise we will create functions for plotting graphs

on a set of points.

a) Make a function plot_line(p1,p2) that takes two points as input arguments

and plots the line between them. The two input arguments should be

lists or tuples specifying x- and y-coordinates, i.e. p1 =(x1,y1). Demonstrate

that the function works by plotting a vertical and a horizontal line.

20

0.0 0.2 0.4 0.6 0.8 1.0

ro

1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

ro

Figure 5.1: Examples of complete graphs from Problem 5.14 b). The left shows

the graph on the corners of the unit square, the right is the graph for eight

equally spaced points on a circle.

b) A complete graph is a graph such that any two points has a line that

connects them. Make a function that takes a list of points and plots the

complete graph on those points. To verify that the function works, first choose

the four corners of the square ((0, 0),(1, 0),(0, 1),(1, 1)) and then the points

(1, 0),(α, α),(0, 1),(?α, α),(?1, 0),(?α, ?α),(0, ?1),(α, ?α), with α =

2/2.

The resulting complete graphs should look like the ones in Figure 5.1.

Hint: Modify the plot_line function from 5.14 a) so that it only calls

plot() but not show(). The complete graph can then be drawn by looping

over the points and calling plot_line for each pair, and finally calling show()

after the loop.

Filename: graph1.py

Problem 5.15. Plotting graphs

Given a natural number n, make a function that plots the following graph:

? Two vertical rows with n points should be placed side by side.

? Each point on the left side should have a line to every point on the right

side and vice versa.

? No two points on the same side should be connected by a single line.

Filename: graph2.py

Problem 5.16. Inefficiency of primality checker

Consider the program from Problem 3.6. Use the timeit module and run

the program to find the time it takes to find a factorization of an n digit

number. Plot the time against the number of digits for the numbers in the file

prime_check.dat. You can use the following code to time the function for

different numbers:

str1 = "f(" + str(n) + ")"

str2= "g(" + str(n) + ")"

N = 100

time1 = timeit.timeit(str1, ’from __main__ import f’,number=N)

time2 = timeit.timeit(str2, ’from __main__ import g’,number=N)

21

Filename: prime_ineff.py

Problem 5.17. Animating a cycloid

One may create a curve by placing a circle on the x-axis, fixing a point on

the circle, and then drawing the trace of the point as the circle is rolling. The

resulting curve is called a cycloid. In mathematical language it is given as

r(θ) = [R(θ ? sin θ), R(1 ? cos θ)]

where R is the radius of the rotating circle and θ is the angle starting at 0 and

increasing.

a) Animate the cycloid as a function of θ starting at 0, ending at 15. Draw a

point at the end of the cycloid that varies with the animation.

Hint: A point can be added through a new plot using for example

point, =axes.plot([],[],’o’)

and updating during the animation.

b) Add the rolling circle defining the cycloid to the plot. You may use that at

a given time θ, the circle is given as s(θ) = (R · θ + cos θ, R + sin θ).

Filename: cycloid.py

22

Problem 5.18. Calibration curve

A tool in chemical analysis for measuring the concentration of a substance in

a sample (e.g. blood or urine), is making a calibration curve. Solutions with

known concentrations of a substance (standard solutions) are measured. The

X-axis is the concentrations of the standard solutions, and the Y-axis is for the

measured intensity of these solutions. The concentration of the substance in

the sample can then be determined using an equation that best describes the

calibration curve. This equation is determined using linear regression.

In Python, you can use the numpy module to obtain a linear regression. How

to do this will be shown.

Assume that you have made five standard solutions with concentrations of

10, 20, 30, 40 and 50 mg/L of the substance you wish to test for. You have used

equipment that detects this substance, and noted down the intensity for each of

your standard samples. The code below shows how the linear regression can be

performed.

# standard concentrations and height of their signals

I_stand = [9.19, 19.8, 27.0, 34.7, 44.9]

conc_stand = [10, 20, 30, 40, 50]

# linear regression

fit = np.poly1d(np.polyfit(conc_stand, I_stand, 1))

conc_curve = np.linspace(0, 60, 100)

signal_curve = fit(conc_curve)

a) Plot the calibration curve as a line, and plot the intensities of the calibration

solutions as points on top of the calibration curve.

The function created in the code above, fit(x), is now a linear function on

the form f(x) = ax + b = y, where x corresponds to the concentration and y

will correspond to the intensity of the signal at the given concentration.

b) Determine a and b by using fit(x) and use these to implement an inverse

function of f(x) such that g(y) = x. Use your function to calculate the

concentration of three different samples of the same unknown compound. The

intensities of the compound of the samples are given below.

I_unknown = np.array([19.9, 20.1, 19.8])

Print out the mean value, with the uncertainty (x ± sN ), of the samples of the

unknown compound. You may use np.mean and np.std to find the mean value

and the uncertainty.

Filename: calibration.py

23

Chapter 6

Dictionaries and Strings

Problem 6.1. A result on primes “dictionarized”

Consider the program from Problem 4.8. Since the entries correspond to each

other, working with two seperate lists is cumbersome. We may avoid that using

dictionaries. Modify the program such that the values are saved in a dictionary

instead of a list. Let the values of n be keys with values π(n).

Filename: primes_dict.py

Problem 6.2. Chemical elements in a dictionary

Consider the dictionary elements_10 consisting of the 10 first chemical elements

of the periodic table:

elements_10 = {1: ’-’, 2: ’Helium’, 3: ’Lithium’,

4: ’Beryllium’, 5: ’Boron’, 6: ’Carbon’,

7: ’Nitrogen’, 8: ’-’,

9: ’Fluorine’, 10: ’Neon’}

a) The chemical elements of number 1 (Hydrogen) and 8 (Oxygen) are missing.

Copy elements_10 into your file, and adjust the dictionary such that the keys

1 and 8 have their correct value. Use the technique as in this example:

dictionary[key] = ’value’

b) Copy the following code into your script, and run the file in your terminal.

Find the difference between the two dictionaries that are printed, and explain

why they are different from each other.

elements_10_copy = elements_10.copy()

elements_10_copy.update({11: ’Sodium’})

print(elements_10)

print(’\n’)

elements_11 = elements_10

elements_11.update({11: ’Sodium’})

print(elements_10)

Filename: chemical_elements_dict.py

24

Problem 6.3. Representation of polynomials

Let f(x) = Pn

i=0 aix

i and g(x) = Pm

j=0 bmx

m be two polynomials. Recall that

a polynomial can be expressed as a dictionary with keys equal to the degree of a

term, and the corresponding coefficient as value (so 3x

2 + 1/2 is represented by

the dictionary {2 : 3, 0 : 1/2}).

You will be asked to implement three functions in this exercise. Check that

each of your functions work as expected by creating two different polynomials

represented as dictionaries, call your functions and print the returned values.

a) Create a function that takes two dictionaries (corresponding to two polynomials

f and g) as arguments and returns a dictionary corresponding to the sum

of the two.

b) Create a function as above that returns the dictionary corresponding two

the product of two polynomials.

c) Add a function that evaluates a polynomial dictionary at a point.

Filename: poly_dict.py

Problem 6.4. Use string operations to create a pretty dictionary

The file atm_moon.txt contains information about the composition of elements

in the lunar atmosphere during nighttime. [Wil17] The values are given in

particles per cubic centimetre.

Write a function that reads the file, and returns a dictionary with the name

of the element as keys and particle density as value. Transform all characters to

their upper case equivalent. Strip off leading and trailing whitespaces in each of

the string keys. Remove all the commas that marks every three digits, and then

convert these values to float numbers.

For example, considering the information ’ Neon 20 -40,000 ’ extracted

from the file, the value of the key ’NEON 20’ should look like this:

’NEON 20’: 40000

Filename: atm_moon.py

Problem 6.5. Interpret output from a program

The program approx_derivative_sine.py calculates an approximation to the

derivative of sin,

for decreasing values of ?x. Direct the output of the program to a file (by

python approx_derivative_sine.py > filename). Write a function that

reads the file and returns three arrays consisting of numbers corresponding to

delta_x, abs_error and n. Plot delta_x and abs_error versus n. Use a

logarithmic scale on the y axis. Explain why the absolute error increases after

n = 8, i.e. after delta_x = 10?8

.

Hint: The function semilogy is an alternative to plot and gives logarithmic

scale on the y axis.

Filename: plot_round_off_error.py

25

Problem 6.6. Saving information in a nested dictionary

The file below contains information about various people. The first column is

the name, the second is the age, and the third is the gender.

John, 55, Male

Toney, 23, Male

Karin, 42, Female

Cathie, 29, Female

Rosalba, 12, Female

Nina, 50, Female

Burton, 16, Male

Joey, 90, Male

a) Create a function that reads the file and returns the information in a nested

dictionary. For example the key ’John’ has the dictionary {’Age’: 55, ’Gender’:

’Male’} as value. When reading the file, the name should read "John", not "John,

".

b) Create a function that takes as an argument a nested dictionary as above,

a person name, and optional arguments: a number (age), and a string (gender)

which returns the same dictionary with the new age and gender. Note that

neither should be changed if no age or gender is given.

c) Extend the function with the possibility to change the name of a person.

One should not be allowed to change a name to one that is taken (Can you think

of a way to allow this without overwriting another person?). Read the file above

and change the gender and name of John. Iterate the dictionary and print the

new information in a table format.

Remark. The nested dictionary here is a prototype of what is known as a class,

and the functions from b) is a prototype for what will be known as methods in

that class.

Filename: people_dict.py

Problem 6.7. Finding the frequency of words in a text

a) Write a function that reads the file RandomWords.dat and finds the frequency

of words of length n. Save the information in a dictionary with the length as

keys and the number of words of that length as values. You may assume that all

words are separated by spaces and that only punctuation marks appear in the

text.

Hint: For your program to be compatible with words of any length, it might

be helpful to use defaultdict imported from collections. See page 339 in the book.

Use the function dict() on such an object to convert it to an ordinary dictionary

b) Write a test function that generates a file of words and checks that the

function returns the correct values.

Filename: word_length.py

26

Problem 6.8. The Euler’s polyhedron formula

Let V, E, and F be the number of vertices, edges and faces in any polyhedron,

respectively. Then, Euler’s polyhedron formula tells us that

V ? E + F = 2.

In this exercise we shall check that the formula works for some given polyhedrons

in a file with this setup:

Polyhedron: cube

vertices: 8 edges: 12 faces: 6

Polyhedron: pyramid

vertices: 5 edges: 8 faces: 5

...

a) Write a function that can read such a file, and returns a dictionary with

the type of the polyhedron as a key, and a dictionary containing vertices, edges

and faces as a value. The function shall strip of leading and trailing whitespaces

in all strings. For example, the value of the key ’cube’ in the dictionary should

look like this:

’cube’: {’vertices’: 8, ’edges’: 12, ’faces’: 6}

Note that you must convert the string numbers to integers, as for example 8, not

’8’. Be aware of the fact that the value of each polyhedron in the dictionary is

again a dictionary. Print the (nested) dictionary that is returned when reading

the file polyhedrons.dat.

b) Write a test function test_polyhedron_formula() that checks that Euler’s

polyhedron formula works for the polyhedrons given in polyhedrons.dat.

Hint: Repeated indexing works for nested dictionaries as for nested lists. Below

is an example of how to access the value of the key ’vertices’ inside the value

of the key ’cube’.

cube_vertices = polyhedrons_dict[’cube’][’vertices’]

Filename: polyhedron_formula.py

27

Problem 6.9. Compute digital roots

Given a number, say 5282, we can compute the sum of the digits. In this case

5 + 2 + 8 + 2 = 17, and doing this again gives 1 + 7 = 8. The one digit number

we get by doing this is called the digital root of the number.

a) Make a function that calculates the digital root of a number.

Hint: Convert the number to a string in order to work with it.

b) Plot the digital root of numbers up to 500 with the digital root on the

x-axis and the frequency of digital roots on the y-axis. Use plt.scatter(x,y)

for the plot.

Filename: dig_root.py

Problem 6.10. Timezone converter

In the file timezones.dat you will find places and their timezone in GMT

format.

a) Make a function that reads the file and saves the information in a dictionary.

b) Create a function that takes local Norwegian time (GMT +1) in the string

format ’ddmmyy-hhmm’, a place, and returns the local time at that place. Your

program should display a message to the user if a place that is not saved in the

dictionary is used. Do the following conversions:

? March 21st 2018 05.34 in Vancouver

? December 31th 2017 20.03 in Sydney

? January 1st 2018 00.15 in London

Filename: timezones.py

28

Chapter 7

Introduction to Classes

Problem 7.1. Saving information in a class

In this problem you can use the program from 6.6.

a) Create a class Person where the constructor takes name, age, and gender

as arguments.

b) Add methods to the class for changing a persons name, age, and gender.

c) Add a method __str__ that returns a string with all the information of

that person. Create an instance of the class, using the information for ’John’ in

the table from Problem 6.6. Change the name and gender of John. Print the

information of the instance before and after changing.

Filename: class_people.py

Problem 7.2. Right triangle class

a) Make a class RightTriangle that represents a right triangle. The constructor

__init__ should take and store two numbers a and b. These are the two

catheti (shortest sides) that define the right triangle. The constructor should

also calculate and store the hypotenuse as a local variable in the class.

Remember that the hypotenuse c relates to the two short sides a and b by

the Pythagorean Theorem:

b) Make an object triangle1 of the class RightTriangle with both short

sides equal to 1. Make another object triangle2 with sides equal 3 and 4.

Check that your implementation is correct by printing the hypotenuse of both

objects. Use the usual object.variable convention to get the hypotenuse.

c) To make a robust program, we often want to code it such that it prevents

being used in ways that does not make sense. In our case, a natural thing to

prevent is making a triangle with sides having negative length, since length is

strictly positive.

29

Make changes to your class such that ValueError is raised if someone tries

to make a triangle with sides of negative length. To test that your class raises

the error correctly, you can test it with the following code:

def test_RightTriangle():

success = False

try:

triangle3 = RightTriangle(1,-1)

except ValueError:

success = True

assert success

d) Add a method plot_triangle() to your class which plots the triangle

when you call it. The corner where the two shortest sides meet should be in

origin. Side a should be along the x-axis, and side b along y. In order to plot,

you might want to find out what the coordinates of each corner must be. Also

make sure to make the axis equally long so that the triangle looks nice (you can

use plt.axis("equal"). Call the plotting method on the instance triangle2

that you created in b).

Filename: right_triangle.py

Problem 7.3. Make a function class

In this problem you will implement a class F which represents the function

f(x; n, m) = sin(nx) cos(mx).

Create the class and let n and m be parameters of the constructor, which

must be stored in the class.

Since the class represents a function, it should be a callable class. An

instance of a callable class can be called on like a function. The special method

for creating a callable class is __call__(self, args). Implement the special

method __call__(self, x) such that it returns the value of the function

evaluated at x.

Create two different instances u and v of class F. Choose values n and m

for both the instances. Plot the two instances u and v evaluated in x against

each other on the interval x ∈ [0, 2π]. In other words, plot u(x_values) on the

x-axis and v(x_values) on the y-axis. Make sure to have enough points on the

interval to ensure that the line is smooth.

Filename: F.py

Problem 7.4. Extending the AccountP class

Modify the class AccountP in the book to include a method transfer that

transfers an amount between two accounts. The method should take an amount

and the account you want to transfer to as arguments. Write a test function

that checks that the methods deposit, withdraw, transfer and get_balance

works properly.

Filename: AccountP.py

30

Problem 7.5. Approximating the square root of two

The square root of two can be represented by a so called continued fraction on

the following form:

In this exercise we will exhibit two possibilities for approximating the number

p

(2).

a) Make a class Square with a method approx_frac that takes an integer n,

an initial value and returns the first n fractions as above with initial value x0.

This can be done by iterating the function

f(x) = 1 + 1

1 + x

starting at x0. For n = 2 and x0 = 1 this gives.

b) Another way to approximate the square is by iterating the function f(x) =12x +2x. Add a method approx_iter that takes a number x0, an integer n,

and returns the value of the function at x0 iterated n times. For n = 2 we would

have f(f(x0)). From here on we assume for simplicity that 0.1 ≤ x0 ≤ 2.

c) Create a method that returns a nicely formatted table with the two approximations

and their difference  along with the exact value for n = 1, 2, 5, 10. Run

the program, which approximation is best?

d) To visualize the approximation plot the exact value as a line in the plane

and the two approximations as points (n, yn), where yn is the approximation.

Use n = 1, 2, 5, 10.

Filename: square_iteration.py

Problem 7.6. Tangent lines on a quadratic curve

Consider a quadratic polynomial on the form f(x) = x

2 + bx + c. At a point x0

the tangent line is given by l(x) = (2x0 +b)x+C where C = f(x0)?(2x0 +b)x0.

a) Make a class Quadratic with a function f(x) (a quadratic polynomial as

above) as initial argument. Make a method that computes the tangent at a

point and returns the function l(x).

Hint: You will need to extract the coefficients b = f(1)?f(0)?1 and c = f(0).

31

b) Create a method that plots the function along with its tangent at a point.

c) Make a method that animates the tangent line moving over the curve f(x).

Make the animation for uniformly distributed x values in the interval ?5 ≤ x ≤ 5.

Test the program with the function f(x) = x

2

.

Problem 7.7. Numerical approximations for the derivative

Let f(x) be a function and f

0

(x) its derivative. There are many ways to

approximate the derivative, some of which are:

a) Make a class Diff with a function f as initial argument and implement

three methods diff1, diff2, and diff3 approximating the derivative using the

above formulas.

b) Create an instance of the class Diff, using f(x) = sin (2πx). Visualise

the difference in accuracy of the three methods for computing the derivative

by comparing the results with the exact function for the derivative, f0(x) =2π cos(2πx). Use the four values h = 0.9, 0.6, 0.3, 0.1, and let x be on the interval

x ∈ [?1, 1].

Filename: class_diff.py

Problem 7.8. Visualizing functions

For a function f(x) we can plot the graph of the function as points (x, f(x)).

This results in a curve in the plane. Suppose we have a function

f(x, y) = (u(x, y), v(x, y)).

The graph of this function lives in four dimensions and is not easily visualized.

On way to visualize these functions is to instead of looking at the graph we look

at how f act on points. For example, how the grid lines in the plane look after

f is applied.

a) First we consider a specific function f(x, y) = (x2 ? y2, 2xy). Write a

program where you define the function f and make a figure with x and y axis

from -2 to 2 where you plot a number of the grid lines in x and y direction in

the same plot. You will need around 15 lines in each direction. Make a another

plot side-by-side in the same figure of all points (x2 ? y2, 2xy) where x and y

are the points in the first plot.

b) To make the construction more flexible, modify your program to be a class

Visualize taking a function f(x, y) = (u(x, y), v(x, y)) as initial argument. It

should contain a method grid(self, n) that generates two plots, one of grid

lines, and one of the image as in a).

32

c) We used grid lines of the plane to see how the function f behaved. We could

have used any curves in the plane. Extend the class with a method circ that

instead of using points corresponding to grid lines, uses circles with expanding

radii. Let the axes go from -5 to 5 and the radii be uniformly distributed

between 0 and 10 (15 circles should be sufficient). Test with the function

f(x, y) = (x2 ? y2 + x + 1, 2xy + y). The second plot should consist of circular

like objects with a self-crossing.

d) Add a method grid_anim that shows an animation of the image of the

functions

f(x, y) = [(1 ? )x + u(x, y),(1 ? )y ? v(x, y)]

where varies from 0 to 1.

e) Using the functions

f(x, y) = (x2 ? y2, 2xy) and g(x, y) = (x2 ? y2 + x + 1, 2xy + y),

test the grid and grid_anim methods on f, and the circ method on g. Use

15 gridlines and 15 circles.

Remark. For the student familiar with complex numbers, this is exactly how

one would visualize a complex function f(z). In our case we can use this for any

function f(x, y), but we usually restrict ourself to look at functions corresponding

to certain complex functions, namely the differentiable ones.

Filename: plot_functions.py

Problem 7.9. A class for coordinates

This exercise will focus on how to implement special methods. You will implement

the class Coords, which represents coordinates in three dimensions.

Hint: the class of complex numbers shown in the book is of similar nature to

the class you should implement in this problem.

a) Create the class Coords. Start by implementing the special methods

__init__(self, args) and __str__(self). The constructor should take

three parameters, x, y and z. The representation of the class should be (x, y, z).

The implementation should be such that the code below works.

sqrt3 = sqrt(3)

close = Coords(1/sqrt3, 1/sqrt3, 1/sqrt3)

far = Coords(3/sqrt3, 15/sqrt3, 21/sqrt3)

print("The coordinates close are at %s" %close)

print("The coordinates far are at %s" %far)

Output:

The coordinates close are at (0.58, 0.58, 0.58)

The coordinates far are at (1.73, 8.66, 12.12)

33

b) Implement the special methods __len__(self) and __abs__(self). The

length of coordinates should always be 3, as the coordinates will be in three

dimensions. The absolute value should yield the Euclidean norm (or the physical

length in space), which is given by

||(x, y, z)|| =px2 + y2 + z2.

The implementation should be such that the code below works.

print("The class Coords represents coordinates in \

%d dimensions" %len(close))

print("\nThe distance from the centre to the point \

close is %.2f" %abs(close))

print("The distance from the centre to the point \

far is %.2f" %abs(far))

Output:

The class Coords represents coordinates in 3 dimensions

The distance from the centre to the point close is 1.00

The distance from the centre to the point far is 15.00

c) Implement the special methods __add__(self, other) and

__sub__(self, other). When adding or subtracting two objects of class

Coords, a new object of class Coords should be returned.

The object returned when adding two coordinates should have the coordinates

xnew = xself + xother

ynew = yself + yother

znew = zself + zother,

and similarly should the method for subtracting should return an object of

Coords with coordinates at

xnew = xself ? xother

ynew = yself ? yother

znew = zself ? zother,

The implementation should be such that the code below works.

further = close + far

print("The coordinates further are at %s" %further)

distance = abs(far - close)

print("\nThe distance from far to close is %.2f" %distance)

centre = further - further

print("\nThe coordinates at the centre are %s" %centre)

34

Output:

The coordinates further are at (2.31, 9.24, 12.70)

The distance from far to close is 14.14

The coordinates at the centre are (0.00, 0.00, 0.00)

Filename: Coords.py

35

Chapter 8

Random Numbers and

Simple Games

Problem 8.1. Throw a die

Compute the probability of getting a 6 when throwing a die. Write a program

that throws a die N times and count how many times the die shows 6, let this

number be M. Then compute the probability of getting a 6 when throwing a

die as M/N.

Filename: die.py

Problem 8.2. Telephone number

A Norwegian telephone number consists of eight digits. We assume that all digits

from 0 to 9 are equally probable in every place of the telephone number. Make

a program that finds the probability of having a telephone number where the

digit 1 appears at least four times.

Filename: telephone.py

Problem 8.3. Coin-flip game

Two persons are playing a simple coin-flip game. They flip a coin in turn, and

whoever first gets a heads wins the game. Make a program to model 100 such

games. Estimate the probability for the first person to flip to win the game.

Filename: coin.py

Problem 8.4. Birthday probability

Make a function that generates a string of random integers between 0 and 9.

Estimate the probability that your birthday is contained in a string of random

numbers of length 1000. Let the format of the date be on the form ddmmyy.

Print the estimates in %.

Filename: birthday_prob.py

36

Problem 8.5. Approximate π by throwing darts

You are throwing darts at a square shaped target with an inscribed circle. Let the

length of the sides of the square be 2, which means that the circle has radius 1.

Assume that you throw the darts such that the darts gets uniformly distributed

on the target. Then, the number of darts which hits the target inside the circle

divided by the total number of darts that hits the target is approximately the

area of the circle divided by the total area of the target. This approximation

gets more accurate the more darts you throw.

number of darts inside circle

number of darts that hits target ≈

area of circle

area of target =π4.

Thus, π can be approximated byπ ≈ 4

number of darts inside circle

number of darts that hits target.

Write a program that throws M darts uniformly on the target. Then approximate

π. Read M from the command line.

Filename: approximate_pi.py

Problem 8.6. Wheel of fortune

At an amusement park they have a wheel of fortune where you can win 2kg of

chocolate. You get to choose one number between 1 and 20 for 20NOK. Assume

that you play on the same number until you win.

a) Write a program that finds the average number of times you have to play

before you win and check if you earn or loose money, compared to buying 2kg of

chocolate in the store.

b) Modify your program so that every time you lose you move one place to

the right, i.e. you increase n by one. If you are at n = 20 you go back to n = 1.

Does this make any difference to the result?

Filename: wheel_of_fortune.py

37

Chapter 9

Object-Oriented

Programming

Problem 9.1. Implement Newtons method

a) Make a subclass Function of the class Diff in problem 7.7 that takes a

function f as an initial variable. It should contain a method such that the

following code is compatible with your program and prints the value of f at 2.

def f(x):

return x**2

func = Function(f)

print(f(2))

b) We would like the class to give estimated values for roots of f. That is,

points such that f(x) = 0. To do this we implement Newton’s formula. It is

given recursively as

xn+1 = xn ?f(xn)f0(xn),

where we give a starting point x0. In some cases(not all) xn will approach a root

of f. Implement this in a method approx_root that takes a starting point and

a bound  < 1 as arguments and approximates xn such that f(xn) < .

Hint: Implement a simple convergence test. Check that f(xn) < 1 after

100 iterations. If not terminate the loop and inform the user that there is no

convergence for that starting point. It is still a possibility for convergence, but

unlikely.

c) Test the program with the function f(x) = x

2 ? 1 and starting value 5 with

bounds 10?i

, for i = 1, 2, 3, 4, 5, 6. Print the approximated value for x, f(x)

and the bound in a table format. Try to run the program with starting value 0.

What happens, can you see why?

Filename: newton.py

38

Problem 9.2. Implement Polynomials as a Class

a) Make a class Quadratic that implements second order polynomials. An object

of Quadratic should be initialised with a list containing the coefficients. Add

a __call__-method that evaluates the parabola at a point x, and a __repr__-

method such that you can print your polynomial.

Make an object of Quadratic with coefficients (1, 3, 2). Print it and evaluate it

at x = 1, x = 2 to make sure it works.

b) Make a subclass Cubic of the class Quadratic that implements third order

polynomials. You should make use of inheritance to extend the class Quadratic

you made in the previous exercise. Implement a method derivative for Cubic.

The method should return an object of type Quadratic that corresponds to the

derivative.

Make an object of Cubic with coefficients (1, 3, 2, 4). Print it and evaluate it at

x = 1, x = 2. Also call derivative and print the result.

Filename: polynomial.py

Problem 9.3. Vectors

a) Make a class Vector2D that implements 2D-vectors(two components). Add

__add__ and __sub__ methods so you can add and subtract vectors. Remember

that vectors add and subtract element-wise, e.i.:

(1, 2) + (4, 5) = (5, 7)

b) For two vectors ~a = (a1, a2) and ~b = (b1, b2), the dot product is defined as

~a ·

~b = a1b1 + a2b2

Implement a method dot that calculates the dot product of two vectors.

Make two vectors ~v = (1, 2) and ~w = (?2, 5) from Vector2D. Print ~v, ~w, ~v + ~w,

~v ? ~w and ~v · ~w

c) Make a subclass Vector3D of the class Vector2D, where Vector3D is extended

with an additional coordinate. Make sure all the methods are updated to work

with three coordinates. The dot product for 3D vectors extends as one would

expect:

~a ·

~b = a1b1 + a2b2 + a3b3

Use inheritance to reuse old code as much as possible(especially for the dot

product method).

Make two vectors ~v = (1, 2, 4) and ~w = (?2, 5, 1) from Vector3D. Print ~v, ~w,

~v + ~w, ~v ? ~w and ~v · ~w

39

d) There is a common vector operation that is defined for 3D vectors, but not

for 2D vectors. This is the cross product. It differs from the dot product in that

the result is not a number, but a new vector:

~c = ~a ×~b

where ~c = (c1, c2, c3). The cross product is define as

c1 = a2b3 ? a3b2

c2 = a3b1 ? b1a3

c3 = a1b2 ? b2a1

Implement a method cross for Vector3D only that calculates the cross product.

Make two vectors ~v = (2, 0, 0), ~w = (0, 2, 0) from Vector3D. Print v × w.

Filename: vector.py

Problem 9.4. Inheritance

In this exercise we will investigate how python handles inheritance by making

some intuitive classes.

a) Begin by making a class Mammal. Add a method info() that returns a

string stating something that is common among all mammals, e.i. "I have hair

on my body". Also add a method identity_mammal() that prints "I am a

mammal".

b) Make a subclass Primate of the class Mammal. Add a method info() that

returns the same as info() for Mammal, in addition to something new specific

for Primates. E.i. "I have a large brain". Do not include the string from Mammal

by copy-pasting, but use inheritance. Also add a method identity_primate()

equivalent to identity_mammal().

c) Now make two subclasses Human and Ape from Primate. Update info()

in the same manner for both Human and Ape, and also add their respective

identity method.

Make an object John of the class Human, and an object Julius of the

class Ape. Try calling info(), identity_mammal(), identity_primate(),

identity_human() and identity_ape() for both John and Julius. Does it

work as you expect? Some of these calls are meant to cause an error.

d) Python is able to check if an object is of a particular class with the

function isinstance. isinstance(object_name, class_name) returns True

if the object object_name is of class class_name. An example could be

isinstance(John, Primate), which returns True if John is of the class Primate.

Use isinstance to check if John is Mammal, Primate, Human and Ape. Do

the same for Julius.

Filename: inheritance.py

40

Appendix A

Sequences and Difference

Equations

Problem A.1. Computing Bell numbers

Let B0 = B1 = 1, the n’th Bell number is defined recursively as.

Make a function that returns the n first Bell numbers and print the first 10.

Filename: bell.py

Problem A.2. Solve a difference equation numerically

We study the difference equation

xn = xn?1 + xn?2.

Write a program that writes out the first 15 elements of the sequence for

x0 = x1 = 1.

Filename: fibonacci.py

Problem A.3. The spreading of a disease

We want to study the spreading of a disease. Assume that people recover at a

rate such that a ratio a of the people that are sick this week will still be sick

next week. It takes two weeks from when you get infected until you become sick,

and a person who is sick will on average infect b people each week, who then

become sick two weeks later.

Let xn be the number of sick people in week n. The number of sick people is

then given by the following difference equation

xn = axn?1 + bxn?2.

a) Write a function disease_weeks(a, b, x0, x1, N) that calculates the

number of people that are ill with a given disease (defined by values a and b)

and returns an array/list containing the number of sick people form the initial

week to week N. (In other words, the function should return an array or a list

with x0, x1, ..., xN ).

41

Let x0 = 100 , x1 = 150 and a = 0.25. Use the function calculate the number

of sick people in an array up to week N = 15 for both b = 0.5 and b = 0.75. Plot

both scenarios in the same plot. Remember to use labels on the curves.

b) We do not need to store all the N + 1 values. Since xn only depends on

xn?1 and xn?2, these are the only values we need to store. Write a function

disease_week_N(a, b, x0, x1, N) that does not use any lists or arrays, and

returns only the number of sick people in week N, xN . Use this function to

obtain the number of sick people in week N = 15 for both the cases you plotted

in a). Print out the results, and remember to include information about which

of the two cases you are printing. Verify that the result is the same as obtained

using arrays.

Filename: disease.py

Problem A.4. Finding π with Newton’s method

It is common knowledge that π ≈ 3.14, but in this exercise you will use Newton’s

method and knowledge of the sine function to improve an approximation of π.

Newton’s method can be written as a difference equation defined

xn = xn?1 ?f(xn?1)f0(xn?1),

and is used to find roots of f(x), based on the function f(x), the derivative of

the function f

0

(x) and an initial guess x0.

Consider the function f(x) = sin(x). Finding the functions derivative f

0

(x),

should be trivial. This function has infinitely many roots, which we know are

located at x = kπ, where k is an integer. This exercise will focus on the root for

k = 1, and use Newton’s method to approximate π. For Newton’s method to

converge towards the correct root, the initial condition, x0, needs to be set close

to the value of π.

Set x0 = 3.14 and calculate x1 and x2 following Newton’s method. Print out

x0, x1 and x2, as well as the value of numpy.pi with 13 decimals. You should

print them in a tidy way such that the values are easy to compare. If Newton’s

method was used correctly, your values for π should improve!

Filename: finding_pi.py

Problem A.5. Find difference equations for computing ln x

The Taylor expansion of ln x is

Introduce sj = S(x, j ? 1) and aj as the two sequences to compute. We have

the initial values s1 = 0 and a1 = (x ? 1).

a) Find the set of difference equations for sj and aj .

Hint: You can find an example on how this is done for e

x

in section A.1.8

in the book.

b) Implement the system of difference equations in a function ln_Taylor(n, x)

which returns sn+1 and |an+1|. The term |an+1| is the first neglected in the sum

and may act as a rough estimate of the size of the error in the Taylor polynomial

approximation.

c) Verify the implementation by computing the difference equations for n = 3

by hand and comparing with the output from the ln_Taylor function. Automate

this comparison in a test function.

d) Check that the accuracy of the Taylor polynomial improves as n increaces

and x is close to 1. What happens when x > 2?

Filename: ln_Taylor_series_diffeq.py

Problem A.6. Lotka-Volterra two species model

We have previously studied the logistic model for poulation growth. This is a

model showing the growth of a population in the abscence of preditors. The

Lotka-Volterra model describes interactions between two species in an ecosystem,

a predator and a prey. We will in the following take the preys to be rabbits and

the predators to be foxes. The number of rabbits and foxes in week n is denoted

by Rn and Fn respectively, and the population is modelled by the equations

Rn+1 = Rn + aRn ? cRnFn

Fn+1 = Fn + ecRnFn ? bFn,

where a is the natural growth rate of rabbits in the absence of predation, b is

the natural death rate of foxes in the absence of food (rabbits), c is the death

rate per encounter of rabbits due to predation and e is the efficiency of turning

predated rabbits into foxes.

Write a program that computes the number of rabbits and foxes up to

n = 500. Use a = 0.04, b = 0.1, c = 0.005 and e = 0.2. In the begining we have

R0 = 100 and F0 = 20. Plot how the number of individuals in the populations

vary with time.

Filename: Lotka_Volterra.py

Problem A.7. Difference equations for computing arctan(x)

The purpose of this exercise is to implement difference equations for computing

a Taylor polynomial approximation to arctan(x). For x ∈ (?1, 1), we have

arctan(x) ≈ S(x; n) = Xn

To compute S(x; n) efficiently, we can write the sum as S(x; n) = Pn

j=0 aj , and

derive a relation between two consecutive terms in the series:

aj = ?x2(2j ? 1)(2j + 1)aj?1.

We introduce sj = S(x; j ? 1) and aj as the two sequences to compute. We have

s0 = 0 and a0 = x.

a) Implement the system of difference equations in a function

arctan_Taylor(x,n). The function shall not use any lists or arrays. (Since

aj only depends on aj?1, and sj only depends on sj?1 and aj?1, these are the

only values that we need to store.) The function shall return sn+1 and |an+1|.

The latter is the first neglected term in the sum (since sn+1 =

Pn

j=0 aj ) and

may act as a rough measure of the size of the error in the Taylor polynomial

approximation.

The function arctan_Taylor(x, n) will give extremely inaccurate approximations

to arctan(x) for x 6∈ (?1, 1). To find a good approximation to arctan(x)

for all x, we can use the fact that.

b) Implement the following piecewise function as a python function

arctan_Taylor_improved(x, n). Your function shall return both the approximation

and the measure of the size of the error as in a).

d) Make a table or plot of the approximation from b) for various x and n to

illustrate that the accuracy improves as n increases.

In the table; include the x value and the order n, the approximation and the

exact value, and the measure of the error.

In the plot; calculate for x ∈ [?20, 20]. Include also the exact function for

comparison, and remember to label the graphs.

Filename: arctan_Taylor_series_diffeq.py

44

Appendix E

Programming of

Differential Equations

Problem E.1. Decrease the length of the time steps

We have the following differential equation

dx

dt =cos(6t)1 + t + x.

Use Forward Euler to solve this differential equation numerically. You should

solve it on the interval t ∈ [0, 10] for n = {20, 30, 35, 40, 50, 100, 1000, 10000}.

Plot all the solutions in the same plot.

Filename: decrease_dt.py

Problem E.2. Implement Euler’s midpoint method

Make a subclass Midpoint in the ODESolver hierarchy from Section E.3 for

solving ordinary differential equations with Euler’s midpoint method. This

method computes

uk+1/2 = uk + ?tf(uk, tk)/2,

uk+1 = uk + ?tf(uk+1/2, tk + dt/2).

Test your implementation by solving y

0 = f(y, x) = cos(x)?x sin(x) and plot the

numerical solutions obtained from both Euler’s midpoint method and Forward

Euler together with the analytical solution. Use 20 time steps on the interval

x ∈ [?5, 5], and y0 = ?5 cos(?5).

Filename: Midpoint.py

Problem E.3. Implement Heun’s method; function

a) Implement Heun’s method specified in formula (E.36-E.37) on page 779 in

the book. Use a plain function heuns_method of the type shown in Sect. E.1.3.

Construct a test problem where you know the analytical solution. Implement a

test function test_heuns_against_hand_calculations() where you compare

your own hand calculated results for u0, u1 and u2 with those in the program.

45

b) Plot the difference between the numerical and analytical solution of your

test problem. You should demonstrate that the numerical solution approaches

the exact solution as ?t decreases.

Filename: heuns_method_func.py

Problem E.4. Solve an ODE until constant solution

Newton’s law of cooling,

dT

dt = ?h(T ? Ts)

can be used to see how the temperature T of an object changes because

of heat exchange with the surroundings, which have a temperature Ts. The

parameter h, with unit s

?1

is an experimental constant (heat transfer coefficient)

telling how efficient the heat exchange with the surroundings is. In this exercise

we shall model the cooling of freshly brewed coffee. First we must find a measure

of h. Supposed we have measured T at t0 = 0 and t1. We can use a rough

Forward Euler approximation of dT

dt with one time step of length t1,

T(t1) ? T(0)

t1

= ?h(T(0) ? Ts)

to make the estimate

h =

T(t1) ? T(0)

t1(Ts ? T(0)).

a) Make a class Problem containing the parameters h and Ts as data attributes.

Let these parameters be set in the constructor. Implement the righthand

side of the ODE in a __call__(self, T, t) method. We will use a class

from the ODESolver hierarchy to solve the ODE, so you shall also include the

terminate function as a method in class Problem (reading the first part of Sect.

E.3.6 may be useful).

Create a stand-alone function estimate_h(t1, Ts, T0, T1) to estimate the h

parameter based on the initial temperature and one data point (t1, T(t1)). You

can use this function to estimate a value for h prior to calling the constructor in

the Problem class.

b) Implement a test function test_Problem() for testing that class Problem

works. It is up to you to define how to test the class.

c) The temperature of freshly brewed coffee is 95? C at t0 = 0 (when it is

poured into your cup), and 92? C after 15 seconds, in a room with temperature

Ts.

For each Ts = [20, 25], solve the ODE numerically by a method of your choice

from the ODESolver hierarchy. Remember to send both the array of t-values and

<instance>.terminate as parameters when calling the solve method (because

terminate=None by default). Plot your two T-arrays of different solutions in

the same plot. If the terminate method works all right, your graphs should be

cut when close enough to the asymptotic values Ts of room temperature.

Filename: coffee.py

46

Problem E.5. Compare numerical methods for solving ODEs

a) Make two subclasses Heun and RungeKutta2 in the ODESolver hierarchy

from Section E.3 in the book for solving ordinary differential equations with

Heun’s method and the 2nd-order Runge-Kutta method.

Heun’s method has the formula

uk+1 = uk +12?tf(uk, tk) + 12?tf(u?, tk+1),

where ?t = tk+1 ? tk is the time step and

u? = uk + ?tf(uk, tk).

The 2nd-order Runge-Kutta method has the formula

uk+1 = uk + K2,

b) Test your implementation in the main block of your program, using y0

(x) =x cos(x) ? sin(x). You shall evaluate y(x) using Heun’s method, 2nd-order

Runge-Kutta method, 4th-order Runge-Kutta method and the analytic solution.

Make one subplot for each method, and one subplot for the analytical solution.

You should solve and plot on the interval t ∈ [?17, 17] using n number of points

on the interval, for n = {20, 25, 50, 150}.

The analytical solution is y(x) = x sin(x) + 2 cos(x), from which you can find

the initial condition of the numerical methods, y0 = y(?17).

Remember to label the subplots with both the names of the methods and

the different values of n.

Filename: compare_methods.py

Problem E.6. Solving a system of ODEs; motions of a spring

The ODESolver hierarchy is adapted to cope with both scalar ODEs and systems

of ODEs. In this exercise we will solve a system of ODEs using ODESolver.

Any ordinary differential equation of n

th order can be written as a system of

1

st order equations.

We will see how a 2

nd order ODE can be used to study the motions of a

spring. We have a block of mass m hanging from a spring. The block is pulled

downwards before it is released, causing the block to oscillate vertically. We will

study the oscillation of the block in this exercise.

If you want to read the mathematical reasoning concerning this physical

phenomena in detail, you may look up section 10.7 Svingninger og resonans in

Kalkulus (Lindstr?m, 2016). This is however not necessary to solve the exercise.

47

a) Consider the following 2

which describes the motion of the spring. The initial condition is

The parameter k is a constant factor that describes the stiffness of the spring.

The parameter q describes the friction, for example air resistance. When q = 0

there is no friction. In this entire exercise we will consider the case where m = 1

and k = 2.

Rewrite the equation to a system of ODEs by hand. Then create a class

ProblemSpring that takes the parameters m, k, and q as instance attributes in

the constructor. The vector u0 should also be defined in the constructor, as it

depends on m, k and q. Create a special method __call__(self, u, t) that

returns the vector

that you calculated by hand.

b) Create a new class SolverSpring that takes problem, T (time stop) and

n (time steps) as parameters in the constructor. The problem attribute is

supposed to be an instance of the ProblemSpring class. Write a method

solve(self, method) that solves our system of ODEs. This is an example of

how it can be done:

def solve(self, method=RungeKutta4):

self.solver = method(self.problem)

self.solver.set_initial_condition(self.problem.U0)

time_points = linspace(0, self.T, self.n+1)

U, self.t = self.solver.solve(time_points)

self.u, self.u_der= U[:,0], U[:,1]

You can choose whether you want to use RungeKutta4 or another method from

ODESolver. Notice how we extract the inital condition U0 from the problem

instance of class ProblemSpring. Also note that the array U contains the

values of the approximations to both u(t) and u

0

(t). Make sure you understand

how to extract the two columns vectors self.u and self.u_der.

c) Write a method exact(self) that returns the analytical solution to our

differential equation. The analytic solution is given by:

d) Write a method plot(self) which plots the exact solution of u(t) together

with your approximations to u(t) and u0

(t).

48

e) In the main block, create two problem instances of ProblemSpring that

represents the two cases where q = 0 and q = 0.3

km, respectively.

Also, create a solver instance of SolverSpring for each of the two problem

instances you just created. Let T = 30 and n = 1000.

Call the solve and the plot methods for both of your solver instances.

f) You will now create two test functions to ensure your implementation of the

problem class and solver class is correct.

Make one test function where you compare the computed solution with the

analytical solution of u(t) for some given parameters, and let the test pass if the

maximum error is less than some given tolerance.

Make another test function where you compare the computed derivative u

0

(t)

with the exact u

0

(t) for the case where m = 1, k = 2 and q = 0. In this case the

analytical solution is

Let the test pass if the maximum error is less than some given tolerance.

Filename: spring_diffeq.py

Problem E.7. Modelling war between nations

We consider the interaction between two nations C1 and C2 and a system of

equations for modelling a conflict between these [Bra13, p.396]. Assuming that

each nation is determined to defend itself against a possible attack, let x(t)

and y(t) denote the armaments of the first and second nation respectively. The

change x

0

(t) depends on the armaments of y(t). We assume that it’s proportional

to it, say ky(t) for some positive constant k. It also depends on the relationship

of the two. Assuming anger leads to increased armaments, let g measure the

relationship between them, positive numbers meaning anger towards the other

nation and 0 means neutral, and negative numbers meaning disarmament. The

cost of having an army will restrain x(t), represented by a term ?αx for some

positive constant α. Similar setup for y(t) yields a system of differential equations:

(t) = 0 we have reached a stable point where neither

nation is increasing armies. We interpret such a fixed point as peace. In the case

were x(t) and y(t) diverges we have an arms race, and we interpret this as war.

a) Make a function that solves the system (E.1) with a numerical method of

your own choice (you may use ODESolver to do this) and a function that plots

the solution curves of x(t) and y(t) for given initial values. Your program should

not solve beyond the point where either x or y is zero. We want to allow the

value zero, so have your program check whether x and y are larger than a very

small negative number. If you use ODESolver this can be done by defining a

terminate function to send with the solve method. Until otherwise specified, we

let t be the time measured in years.

Filename: C_model.py

49

b) Modify your program to instead consist of two classes. The first class

ProblemConflict should contain the following:

? An init method saving all information relevant to the problem (parameters

etc)

? A call method such that the class can be called as a function. It should

take an array [x, y] of specific values of x and y at time t, the time t, and

return the right hand side of the ODE system.

The second class Solver should consist of the following:

? An init method that takes a problem instance on the form above, and a

step length dt.

? A method that solves the problem, with the same restrictions as in Problem

a). It should solve any problem on the same form as ProblemConflict

that is given by two differential equations.

? A method that plots the solutions as in Problem a).

? A method that saves an image of the plot in .png format. When this is

called, no plot should be visible to the user.

Use the parameter values α = β = 0.2, g = h = 0, x0 = 10000, y0 = 20000. Run

the program once with k = l = 0.2, and once for k = l = 0.3 plotting the first 10

years. What is the interpreted difference between these two?

Hint: You may need to convert step length to the number of time points to

use. This can be done as

n =int(round("Last time step"/ "Step length"))+1

Filename: C_model_class.py

c) Let us consider the parameter values k = l = 0.9, α = β = 0.2, and

g = h = 0. One can argue that these give rough estimates for the arms race

between 1909 to 1914 between the alliances of France and Russia, and Germany

and Austria-Hungary [Bra13]. Assuming stability prior to this and negligible

armies, we assume x0 ≈ 0 and y0 ≈ 0 (This does not mean that neither nation

had armies, but that they were much smaller prior to the arms race). Solve the

problem when x0 and y0 are zero versus when they are small positive numbers.

Plot the next 5 years of both in the same figure. What happens?

Filename: C_model_c.py

d) So far we have seen a model intended to describe a conflict situation prior

to war. The preceding model doesn’t describe what happens during a war, but

similar equations can.

First of all, we will work with two types of warfare. The conventional

one, what we know as regular warfare, and guerrilla warfare, where groups

of combatants use military tactics such as ambush, raids, hit-and-run, among

others.

We first consider two conventional armies engaging. Let x(t) and y(t) denote

the respective forces (the number of soldiers) and t denote the time measured

in days. The rate of change of x(t) is affected by combat loss, operational loss

50

(non-combat related. e.g. disease, accidents), and reinforcements. Combat loss

should be proportional to the size of the opponent, represented by a term ?αy(t),

α > 0. The operational losses should depend only on x(t), represented by a

term ?kx(t), k > 0. The reinforcements are represented by a function f(t). In

short-term warfare, the operational losses are negligible, and we will assume it

to be zero. We get equations

dx

dt= ?αy(t) + f(t),

dy

dt= ?βx + g(t).

For a conventional-guerrilla combat, y(t) representing the guerrilla army, we

assume that the combat losses also depend on the size of its own army. As

guerrilla armies often use strategies of surprise and hidden attacks, it is safe

to assume that bigger losses are experienced when the army is larger. Let

?βx(t)y(t) denote the combat loss of a guerrilla army. By the same arguments

as above, we get equations.

Make two classes ProblemCCWar and ProblemGCWar on the same form as

ProblemConflict representing the new problems. Note that f and g are now

functions. To handle the case of the provided f and/or g not being functions,

you may need the commands callable(f) which checks if f a callable, and

isinstance(f, (float, int)) that checks if f is a float or int, in order to

convert a constant to a constant function.

Filename: CW_model.py

e) The battle of Iwo Jima is a famous battle during World War II. It was fought

on an island just outside of Japan. America invaded the island on February 19,

1945, and the fight lasted for 36 days. The Japanese army consisted of around

21500 soldiers, while the Americans had a number above 50000 by the 36th day.

During the war, the Japanese had no reinforcements. The Americans started

with no soldiers, but landed 54000 soldiers the first day, 6000 the third, 13000

the sixth, and none for the remaining. The reinforcements is therefore given as

It can be shown that good estimates for the parameter values are α = 0.0544

and β = 0.0106 [Bra13]. The exact values on a day to day basis is given in the

file Casualties.dat. Plot the modeled American army vs the exact numbers,

and y(t) vs x(t). Both plots should have the x-axis corresponding to the first

T = 36 days.

Filename: iwo_jima.py

51

f) Find the least number of soldiers Japan would need (according to the model)

to have won the fight. You may round to the nearest hundred. Hint: Check

which army decreases to zero first. You might want to extend the variable T for

this.

Filename: least_number.py

g) Suppose the Japanese army was interpreted as a guerrilla army. Find a

value for β such that the fight is close. Is it likely that the outcome would be

different if America met a large guerrilla army?

Filename: guerrilla.py

Problem E.8. Simulate the spreading of a disease

In this exercise we will model epidemiological diseases by implementing the

SIRD model. Suppose we have four categories of people: susceptible (S) who are

healthy but may contract the disease, infected (I) who have developed the disease

and can infect the susceptible population, recovered (R) who have recovered

from the disease and become immune, and the deceased (D) who did not survive

the disease. Let S(t), I(t), R(t) and D(t) be the number of people in category S,

I, R and D, respectively at time t. We have that S(t) + I(t) + R(t) + D(t) = N,

where N is the size of the initial population. For simplicity, we will assume that

the populations otherwise remains constant.

Normal interaction between infected and susceptible members of the population

causes a fraction of the susceptible to contract the disease. The fraction of

the susceptible population that becomes infected will depend on the likelihood

of encountering an infected individual as well as how contagious the disease is.

This will be proportional to the number of infected members of the population.

During a time interval ?t starting at time t, the fraction of the susceptible

population which contracts the disease is αI(t)?t. The number of people that

moves from the S to the I category is given by

S(t + ?t) = S(t) ? αS(t)I(t)?t.

Divide by ?t and let ? → 0 to get the differential equation:

(t) = ?αS(t)I(t). (E.2)

Per time unit a fraction β of the infected will recover from the disease, and

a fraction γ of the infected will die as a result of the disease. In a time ?t,

βI(t)?t people of the infected population will recover and move from the I to

the R category, and γI(t)?t dies and move from the I to the D category. In the

same time interval, αS(t)I(t)?t from the S category will be infected and moved

to the I category. The accounting for the I category therefore becomes

I(t + ?t) = I(t) + αS(t)I(t)?t ? βI(t)?t ? γI(t)?t,

which in the limit ?t → 0 becomes the differential equation

(t) = αS(t)I(t) ? βI(t) ? γI(t). (E.3)

The R category gets contributions from the I category:

R(t + ?t) = R(t) + βI(t)?t,

and the corresponding ODE for R reads

(t) = βI(t). (E.4)

Finally, the D category gets contributions from the I category as well:

D(t + ?t) = D(t) + γI(t)?t,

and the corresponding ODE for D reads

D0

(t) = γI(t). (E.5)

The system (E.2)-(E.5) is what we will call a SIRD model.

a) Make a function for solving the system of of differential equations in the

SIRD model by a numerical method of your choice from the ODESolver. Make

a separate function for visualising S(t), I(t), R(t) and D(t) in the same plot.

Make sure to make labels.

b) Adding the equations shows that S, which means

that S(t) +I(t) +R(t) + D(t) must be constant. Perform a test at each time step

by checking that S(t)+I(t)+R(t)+D(t) equals S(0)+I(0)+R(0)+D(0) within

some small tolerance. Since ODESolver is used to solve the ODE system, the

test should be implemented as a user-specified terminate(u, t, k) function

that is called by the solve function at every time step. The terminate function

should simply print an error message and return True for if S + I + R + D is

not constant.

We will now consider the spreading of one of the most devastating pandemics

in human history, the Black Death. The Black Death, also called the Plague,

was evidently spread to Norway in 1349 after a ship from England arrived in

Bj?rgvin (today Bergen), carrying the disease.

There lived approximately 7000 people in Bj?rgvin in 1349. [Hor14] Let’s say

the ship crew consisted of 30 men, which were all infected with the Plague. For

simplicity we only consider human-to-human transmission of the disease (we do

not consider the rats and fleas). We are interested in how the disease developed

in Bj?rgvin the first 8 weeks after the ship arrived.

We assume that the Plague was 90% deadly and that death usually occurred

four days after being infected.

c) Visualise first how the disease develops when α = 6.5 · 10?5 and print out

the number of deceased after 63 days.

Certain precautions, like staying inside, will reduce α. Try α = 5.5 · 10?5

and compare the plot with the plot where α = 6.5 · 10?5

. Comment how the α

influences S(t).

Use the constants β = 0.1/4 and γ = 0.9/4 to describe the Plague in Bj?rgvin.

The initial conditions would be S(0) = 7000, I(0) = 30, R(0) = 0, D(0) = 0,

?t = 1, and t ∈ [0, 63]. Time t here is given in days.

As there do not exist any exact data from the condition in Norway during

the Black Death, the parameters above are all fictional. However, there

53

is a broad consensus that the disease killed more than half of Norway’s population.

https://sml.snl.no/svartedauden.

Filename: bjorgvin.py

Problem E.9. Introduce classes in the SIRD model

Implement the SIRD model from exercise E.8 in a module called SIRD.py. First

we will create a class Region which can represent any region given the specific

initial conditions. Then we will create a problem class ProblemSIRD and a solver

class SolverSIRD to solve the SIRD system of differential equations for a given

region.

a) Create a class Region which has three methods, a constructor, a method

set_SIRD_values(self, u, t) and a method for plotting the SIRD values

plot(self, x_label).

The constructor should take in the name of the region and the initial conditions

S(0), I(0), R(0) and D(0). All five parameters should be stores as

attributes in the class. You will also need an attribute self.population which

is the total population of the region at t0.

The method set_SIRD_values(self, u, t) should take out the SIRD

values from u and store S, I, R, D and t as attributes of the class.

The method plot(self, x_label) should plot S, I, R and D in the same

plot. The string for the plt.xlabel should be given as a parameter (as the

units of time may vary), while the plt.ylabel should always be set to e.g

"Population". Set the title of the plot to be the name of the region. Specify color

and label for all the different SIRD categories (an example could be

plt.plot(self.t, self.S, label=’Susceptible’, color=’Blue’)). Do

not include plt.legend() or plt.show() in the code. This is because we may

later want to add labels to the graphs, and because will use this method to plot

several subplots.

In a main block in the bottom of the file, create an instance of class Region

called bjorgvin using the parameters found in Problem E.8.

b) Create the class ProblemSIRD. The constructor should take in the parameters

α, β, γ and a region, which must be an instance of the class Region.

The parameter α in the SIRD model can be constant or function of time. The

implementation of ProblemSIRD should be such that α can be given as either a

constant or a Python function. The constructor should therefore look like this:

def __init__(self, region, alpha, beta, gamma):

if isinstance(alpha, (float, int)): # number?

self.alpha = lambda t: alpha # wrap as function

elif callable(alpha):

self.alpha = alpha

# Store the other parameters

self.set_initial_condition() # method call

Write the method set_initial_condition(self) which stores a list

self.initial_condition containing the initial values of S(0), I(0), R(0), D(0)

(in this particular order). The initial values should be extracted from the class

attribute region.

54

Write a method get_population(self) which simply returns the value of

the population of the region, which is stored in the class attribute region.

Write a method solution(self, u, t) which simply calls the method

set_SIRD_values(u, t) of the class attribute region.

Finally, write a special method __call__(self, u, t) which represents

the right-hand side function of our SIRD system of ODEs. This method will

return a list of the calculated values of S.

In the main block, create an instance of class ProblemSIRD called problem

using the parameters found in Problem E.8 and the Region bjorgvin.

c) Now we will create a class SolverSIRD. The class constructor should take

the parameters problem (which must be an instance of class ProblemSIRD), T

(final time) and ?t, and store them as attributes. The constructor should also

store an attribute called total_population, which is obtained by calling the

get_population method of problem.

Implement a method terminate which at each time step t checks that

S(t) +I(t) +R(t) + D(t) equals total_population within some small tolerance.

Simply print an error message and return True for termination if the total

population is not constant. As the ODESolver will be used to solve the ODE

system, this method will be called by the solve method in ODESolver at every

time step.

Write a method solve(self, method) that solves the SIRD system of

ODEs by a method of your choice from the ODESolver. Use the following sketch

for this method:

def solve(self, method=RungeKutta4):

solver = method(self.problem)

solver.set_initial_condition(...)

t = np.linspace(...)

# Remember to "activate" terminate in the solve call:

u, t = solver.solve(t, self.terminate)

# set the values of S, I, R, D, and t via

# Problem class to the Region class:

self.problem.solution(u, t)

In the main block, create an instance of class SolverSIRD called solver using

the parameters found in Problem E.8 and the ProblemSIRD problem. Plot the

results by calling the plot(x_label) method of the Region bjorgvin. Label

the graphs by calling plt.legend() before you call plt.show().

Filename: SIRD.py

Problem E.10. The SIRD model across regions

The problem class from exercise E.9 will only be able to model the spreading of

a disease within one region. In this exercise we will improve our program with

subclasses of ProblemSIRD and Region that permits people in one region to get

infected by people from another region. The likelihood of transmission of disease

across regions will depend on the distance between the regions.

We now denote the categories to specify which region they belong, such that

e.g. Si(t) would be the number of susceptible in the i-th region at time t.

(t) belong to the i-th region. In this new,

interacting SIRD model the expressions of R0

(t) must consider the possibility of members of the

population in the i-th region contracting the disease due to contact with another

region. The expression of the i-th region is given by,

where M is the number of regions and dij is the distance between the i-th and

the j-th region. Note that the distance from a region to itself, dii, is always zero,

which leaves the expression unchanged from the previous SIRD model.

The derivative for the infected category is then.

a) Create a subclass of Region called RegionInteraction. In addition to the

parameters in Region, the RegionInteraction needs two parameters latitude

(φ) and longitude (λ). The constructor should store the two values and convert

them from degrees to radii, which is done by multiplying by π

180? . Use the super

class to store the rest of the parameters as attributes.

Create a method distance(self, other) which calculates the distance

between the self region (i) and another region (j). The distance between two

regions is given by the arc length d between the regions:

dij = REarth?σij ,

where the radius of the Earth is REarth = 64 given in units of 104 m and ?σ is

given by

?σij = arccos

sin φi sin φj + cos φi cos φj cos (|λi ? λj |).

b) Create a subclass ProblemInteraction of class ProblemSIRD. This class

should take in a list of regions that must be instances of RegionInteraction.

In adition to all the same parameters as its super class, the constructor of

ProblemInteraction should take in and store region_name. Send all other

parameters to the constructor in the super class.

The method get_population(self) should store the total population of

all the regions combined as self.total_population.

The method set_initial_condition(self) must create a (not nested)

list self.initial_condition with the initial values from all the regions. Loop

over all the regions in the list region to get the list on the form

[S1(0), I1(0), R1(0), D1(0), S2(0), I2(0), R2(0), D2(0), ..., DM(0)].

The special method __call__(self, u, t) should return a list with the

derivatives at time t, in the same order as the list self.initial_condition.

Below is a sketch of the implementation could look like:

56

def __call__(self, u, t):

n = len(self.region)

# create a nested: SIRD_list[i] = [S_i, I_i, R_i, D_i]:

SIRD_list = [u[i:i+4] for i in range(0, len(u), 4)]

# crate a list containing all the I(t)-values:

I_list = ...

derivative = []

for i in range(n):

S, I, R, D = SIRD_list[i]

dS = 0

for j in range(n):

I_other = I_list[j]

dS += ...

# calculate dI, dR and dD

# put the values in the end of derivative

return derivative

The method solution must provide the SIRD lists to all the regions. The

example below shows how a nested list where the each element in the list contains

the SIRD lists for a certain region. You do not have to use the code, but the

result must be the same.

def solution(self, u, t):

n = len(t)

SIRD_list = [u[:, i:i+4] for i in range(0, n, 4)]

for i in range(n):

SIRD = SIRD_list[i]

part_in_region = ...

part_in_region.set_SIRD_values(SIRD, t)

# store t as an attribute

# store S, I, R and D as attributes

The attributes self.S, self.I, self.R and self.D must be the total values

for all the regions combined.

Create a new method plot(self, x_label). the method should create

the same kind of plot as class Region’s method plot(self, x_label), as

explained in Problem E.9. The method in ProblemInteraction should plot

the SIRD values for all the regions combined, and the title of the plot should be

self.region_name.

Filename: SIRD_interaction.py

Problem E.11. Simulate the spreading of the Plague in Norway

In this exercise we will use the classes ProblemInteraction, SolverSIRD and

RegionInteraction from Problem E.10 to simulate the spread of the Black

Death in Norway.

We assume that the disease first broke out in Bj?rgvin in 1349, and that this

was the only case of the Plague in country at that time. We assume that the

disease ravaged for two years (104 weeks).

We divide Norway into five regions: Vestlandet, S?rlandet, Tr?ndelag,

?stafjells and Nord-Norge. We assume that there lived about 370000 peo-

57

ple in Norway; 90 000 in Vestlandet, 65 000 in S?rlandet, 80 000 in ?stlandet,

70 000 in Tr?ndelag and 65 000 in Nord-Norge.

The positions of the regions are given by the latitude ang longitude of certain

city in each region, which are given in the table below:

Region city latitude longitude

Vestlandet Bj?rgvin: 60° 5.3°

S?rlandet ?ysleb?: 58° 7.6°

?stlandet Brandbu: 60° 11°

Tr?ndelag Steinkjer: 64° 11°

Nord-Norge Bardufoss: 69° 19°

a) Create a function for simulating the Plague in Norway. The function must

contain one instance of class RegionInteraction for each region. Let S(0) be

the population in each region. Except for I(0) = 30 in Vestlandet, let all the other

initial conditions be set to zero. Create a list of the five RegionInteraction

instances.

The function will plot one subplot of the disease progress of each region, and

also one subplot for the total progress for all regions combined. The function

could look like this:

def plague_Norway(alpha, beta, gamma, num_Weeks, dt):

# create all five instances of regionInteraction

# create a list containing all five regions

# create problem, instance of ProblemInteraction

# create solver, instance of SolverSIRD & call method solve

plt.figure(figsize=(.., ..)) # set figsize

index = 1

# for each part in problem’s attribute region:

plt.subplot(2, 3, index)

# Call plot method from current part

index += 1

plt.subplot(2, 3, index)

# Call plot method from problem

plt.legend()

plt.show()

# print the total percentage of deceased after two years

Hint: Write the percent sign twice (%%) to get % in a string.

Call the function using the parameters α = 7 · 10?6

, β = 0.1/4 and γ = 0.9/4.

The unit of time is weeks. Solve for 104 weeks and set ?t = 1/7 such that one

time step represent one day.

b) Until now we have assumed that α is constant. α is the parameter which

describes the possibility that one susceptible meets an infected during the time

interval ?t with the result that the infected "successfully" infect the susceptible.

In reality, this α may probably not be constant. In times of bad weather, more

people would stay at home and α would be lower. Other factors could also

decrease alpha, like better hygiene and wearing masks over a certain period of

time. On the other hand; factors like nice weather, festivities and other reasons

58

for people to gather would increase α. As the Plague ravaged in Norway so long

ago it is hard to reproduce an accurate approximation of α. Let us therefore

assume that the weather and other factors made α look like this piecewise

function:

Implement α(t) as a piecewise function alpha(t). Call plague_Norway using

the piecewise function for α(t). Keep the rest of the values from part a).

You may notice that by this model, the Plague barely reached Nord-Norge

at all. In fact, we aren’t even sure today whether the Plague ravaged in the

Northern part of Norway or not.

c) In the Norwegian legends, Pesta was the literal personification of the Plague.

Pesta was depicted as an old woman who wandered the countryside, carrying

either a rake or a broom. Legend has it that if she passed by your home while

carrying a rake, someone in the household would be infected by the plague while

the rest would be spared. However, if she was carrying a broom then there were

not much hope for anyone.

Implement a function pesta(t) where you generate a random number between

0 and 20. You can use the randint function imported from the Random

module like this: number =randint(0, 20). If the generated number is 13,

Pesta is at your door carrying a broom. The function should return ten times the

value of alpha(t) from part b). If the generated number is 4, then Pesta has

brought her rake. The function should return five times the value of alpha(t).

Otherwise, Pesta is not around today, and your function should return 0.4 times

the value of alpha(t).

Call plague_Norway using the pesta(t) function for the parameter α. Keep

the rest of the values from part a). In which region was Pesta most present?

Filename: plague.py

59

Bibliography

[Bra13] M Braun. Differential Equations and Their Applications: An Introduction

to Applied Mathematics. Vol. 15. Springer Science & Business

Media, 2013.

[Hor14] Kulturhistorisk vegbok Hordaland. Middelalderbyen Bergen. https:

//digitaltmuseum.no/011085442898/middelalderbyen-bergen.

2014.

[Ins19] Norwegian Meteorological Institute. Homogenized daily values. https:

//eklima.met.no. 2019.

[Wil17] Dr. David R. Williams. Moon Fact Sheet. https://nssdc.gsfc.nasa.

gov/planetary/factsheet/moonfact.html. 2017.

60