MATH2601: Computer Tutorials and Assignments
Academic Year 2024/25
Organisation
Tutorials
In the tutorials, we will introduce new tools and concepts in Python, in order to use methods from numerical analysis to solve a range of mathematical problems.
Assignments
Computer assignments will be set fortnightly. During the Computer Practical classes on Monday, you will be able to work independently on these assignments, under the supervision of academic staff.
You should submit your work by 4.pm on the Monday in the week following the week of the Practical class to receive feedback. Write your work in a Python Jupyter Notebook called xyz_N .ipynb, where xyz is your ISS username (i.e., your login) and N the number of the assignment (1, 2, 3, 4, 5), then submit it to the Assessment area for MATH2601 on the VLE. Jupyter Notebooks will naturally allow you to write a short report, i.e., to produce an analysis of your results and to comment on them.
Assessments
Two of the computer assignments, 3 and 5, will be assessed. They will be graded and you will get feedback on your work. The overall Computer Assessment grade counts for 30% of the MATH2601 grade and students need to achieve at least 40% on the Computer Assessments to pass the module.
Chapter 1
Computational accuracy
1.1 Tutorial
1.1.1 Motivation
Suppose we wish to evaluate exp x to a very high accuracy. As discussed in the lecture, floating-point arithmetic provides only limited precision, which will limit the accuracy of calculation, so we need an alternative approach. One is based on the Taylor series
(1.1)
If we truncate the series and x is a rational number, then the right-hand side is also a rational number. In Python, integers can be arbitrarily large and there is no round-off error, so we have no problem, in principle, representing any accuracy in this way.
To make this work, we need to introduce a new type to represent rational numbers, complementing the int and float types in standard Python. We will achieve this be defining a new class. In fact, the real purpose of this lecture is to learn the Python class structure, which is a powerful and flexible programming technique.
1.1.2 Classes
Python has types like int and list, and instances of those types like 42 and [1, 4, 9]. A class is a mechanism for defining a new type. A class defines attributes (data associated to every object of that type) and methods (functions that act on the type). In this tutorial, we will define a class called Rational, the object are the rational numbers, the attributes are the numerator and the denominator (which define the rational number) and the methods are operations like addition.
class Rational:
def __init__(self, num, den): # constructor: x = Rational(num, den)
self.num = num # num = numerator
self.den = den # den = denominator
def display(self):
print(str(self.num) + ‘ / ‘ + str(self.den))
This is a first attempt to create the class Rational. It defines two methods, a special method called __init__ and a normal method called display. Here is how to use the new class:
In [1]: two_thirds = Rational(2, 3) # construct variable of type Rational
In [2]: two_thirds.num # retrieve attribute called `num`
Out[2]: 2
In [3]: two_thirds.display() # call `display` method
2 / 3
The special __init__ method must appear in any class definition and is used to construct instances of the class. When we write Rational(2, 3), the __init__ method of the class Rational is called. The first argument, customarily called self, is the instance being constructed. The other arguments are given by the user (in the example, they are 2 and 3). In our class, the __init__ method takes the arguments passed by the user and assigns them to the num and den attributes of the instance which is being constructed.
Once the new instance is constructed and assigned to the variable two_thirds, we can refer to the two attributes of the new object created as two_thirds.num and two_thirds.den. We use a similar syntax to call the display method, which displays the rational number on the screen. Note that the definition of the display method has one argument, self. This argument is set automatically by Python to the instance that the method is called on. The user should call the display method without any arguments.
Before we go any further, let us improve our class definition in two ways. Firstly, we ensure that rational numbers numbers are expressed in their lowest possible form, by dividing through by the highest common factor (hcf). Secondly, we want to display rational numbers using the print function, just like any other variable in Python. We can achieve this by defining another special method, called __repr__. Whenever Python prints a variable, it looks for this method to see how to display the variable.
class Rational:
def __init__(self, num, den):
hcf = gcd(num, den)
self.num = num // hcf
self.den = den // hcf
def __repr__(self):
return str(self.num) + ‘ / ‘ + str(self.den)
We now wish to endow our class with other methods, in order to manipulate rational numbers. For example, let us define a method for multiplying two rational numbers. Again, there is a special method which Python calls if we use the * operator. This method is called __mul__. We can define this method in our class as follows:
def __mul__(self, rhs):
new_num = self.num * rhs.num
new_den = self.den * rhs.den
return Rational(new_num, new_den)
The first argument is self and represents the instance on the left-hand side of the * operator, while the second argument is the right-hand side. To use this method we would need to type:
In [1]: x = Rational(3, 5)
In [2]: y = Rational(4, 7)
In [3]: x * y
12 / 35
In the Practical session, we will further extend the Rational class.
1.2 Assignment
1.2.1 Accuracy
The math module includes a function exp which computes the exponential of a floating point number. Import this function, and note how many decimal places are given in the evaluation of exp(1). Using Taylor series to approximate the exponential function, we could define the following function:
def exp_taylor(x, n):
exp_approx = 0
for i in range(n):
exp_approx += x ** i / factorial(i)
return exp_approx
which computes the value of n terms of the sum in equation (1.1) from the tutorial notes. How many terms are needed before this function offers no further increase in accuracy in the evaluation of exp(1)?
1.2.2 Defining the Rational class
Download the file rational_start.py from the VLE, which is the Rational class discussed during the tutorial. Test that it works by trying commands like in the tutorial.
Your task is to extend this class. Add a method __add__ which defines the + operation (note that it must be called __add__), which produces the sum of two rationals, self and rhs.
Next, add a method __neg__ which defines the negation operation -. (Remember, you only need to negate the numerator or the denominator!). Now define the subtraction operation (also -) with a method __sub__. You can use the fact that subtraction is the same as adding the negation.
After you defined the __sub__ method, the __gt__ method should work. Test that you can compare rational numbers using > and <. Finally, test the methods floor, frac and to_float. Check you understand what they do, and how they do it.
1.2.3 Partial sums
Adapt the function exp_taylor to compute the first n terms of exp x in rational numbers. The first term is then given by Rational(1, 1). You probably need to compute a rational number for 1/n! for each term of the sum.
Use this function to give the first ten rational approximants to e.
1.2.4 Arbitrary accuracy
Download the file expand.py from the VLE which contains the code for the three methods: expand, expand_integer_part and expand_fractional_part. Place them in the Rational class and check you understand how they are used to convert a rational number to any base (like decimal), and to any number of digits. For example,
In [1]: Rational(1, 4).expand(10, 4)
Out[1]: '0.25 (base 10)'
In [2]: Rational(1, 7).expand(10, 30)
Out[2]: '0.142857142857142857142857142857 (base 10)'
Experiment with these conversion functions. What is the 100th decimal place of e? How many terms in the partial sum are necessary to establish this?
1.2.5 Approximations of π
Leonhard Euler gave the following expression for π:
Write a function which computes the first n terms of this series, and use it to confirm that the 50th digit in the decimal expansion for π is 0.
1.2.6 Advanced material
Srinivasa Ramanujan was a self-taught Indian mathematician who had a remarkable affinity for numbers. He discovered many extraordinary series for π. For example,
where k!! is the double factorial, given by
k!! = k(k − 2)(k − 4)· · · 7 × 5 × 3 × 1 for odd k.
Write a function to compute the nth partial sum of this series, using the Rational class, and confirm that the 50th digit of π is zero. How many terms in the series are necessary to find this digit?
The current record for computation of decimal places of π is 62.8 trillion (6.28 · 1013) digits (established in August 2021). This was computed using a variation of the following infinite series given by Srinivasa Ramanujan in 1910:
Write a function to compute the nth partial sum of this series using the Rational class. You will need to obtain an exact rational approximation to √2, perhaps using continued fractions. How many terms are required to obtain the 50th digit of π?
Chapter 2
Root-finding
2.1 Tutorial
One topic that is endlessly useful in computational mathematics is root-finding. That is, given a function f, to find the value(s) of x for which f(x) = 0. In simple cases the root can be found explicitly by algebra, but in general a numerical technique is needed. There are several different ways to approach this problem, and in this workshop we investigate three well-known ones — bisection, secant, and Newton. Since we will always be finding roots of functions, we will take this opportunity to review and extend our knowledge of Python functions.
2.1.1 Functions
By now we are familiar with defining functions in Python, using the formulation
def function_name(argument_1, argument_2, argument_3, ...):
return value
The point of using a function is that often a block of code is to be used repeatedly, with varying input parameters. Functions can help organise a large program coherently, and can aid in debugging. A function always returns something; if there is no return statement, the function returns None.
2.1.2 Default argument values
It is possible to define a default value for one or more function arguments. For example,
def test(a, b=1, c='True'):
print(a, b, c)
which will produce
In [1]: test(4, 5, 6)
4 5 6
In [2]: test(3)
3 1 True
Note that the function test can now be called with fewer arguments than it is defined to allow. Here the arguments can be thought of as mandatory (in the case of a), or optional (in the case of b and c). The optional arguments need not be given, but if they are, they are viewed in order:
In [3]: test(4, 'False')
4 False True
One important point is that the default value is evaluated only once. For example, the function
def f(a, L=[]):
L.append(a)
return L
will produce
In [4]: print(f(1)); print(f(2)); print(f(3))
[1]
[1, 2]
[1, 2, 3]
since the list L is not emptied each time the function is called.
2.1.3 Variable numbers of arguments
Occasionally you might want to define a function which accepts different numbers of arguments. As a very simple example (and assuming that there is no better way to accomplish this), suppose we want a function to compute the mean of n numbers, where n may change. Python has a syntax to achieve this.
def mean_of_n_numbers(*args):
mean = 0
for i in args:
mean += i
return mean / len(args)
We can call this function as: mean_of_n_numbers(1, 2, 3, 4), for example, or with any number of arguments. Inside the function, args is interpreted as a Python tuple. Similarly, the function can be called using
numbers = [1, 2, 3, 4, 5]
mean_of_n_numbers(*numbers)
The * unpacks lists of variables into positional variables. Of course, this simple example could have been coded using only Python lists, but the *args syntax is a standard method to cope with arbitrary numbers of arguments.
2.1.4 Lambdas
There is another mechanism for defining functions in Python, designed for relatively simple, one-line procedures. This alternative to def is called a lambda. As an example, the equivalent of the function
def cube_plus_a(x, a):
return x ** 3 + a
is given by lambda x, a: x ** 3 + a.
The lambda formulation does not give a name to the function, so they are often called anonymous functions. Lambdas are intended for function declarations which will typically be used once, and then forgotten. Their simplicity is manifested in the fact that they must be a single expression only. It is not straightforward to define exactly what is meant by a “single expression”, but it must be a statement which returns something. Assignment statements (such as x = 3) do not return anything, so cannot be used in a lambda.
Lambdas do not offer anything that is not offered by def functions, apart from arguable simplicity, and so their use is largely down to personal preference. However, one natural place to use a lambda is when passing an arbitrary function to another function, such as a root-finder. Here is a function which does fixed-point iteration until two successive iterates differ by less than the user-supplied tolerance:
def fixed_point_iteration(function, x_0, tol):
x_old = x_0
x_new = function(x_0)
while abs(x_new - x_old) > tol:
x_old = x_new
x_new = function(x_old)
return x_new
To call this function we might first define the function we wish to finds roots of:
def my_function(x):
return 1 / x + x / 2
and then call the function with fixed_point_iteration(my_function, 1, 1e-10), to find a fixed point to within 10−10 with an initial guess of 1. However, if we want to repeatedly find roots of different functions, it will get tedious having to define a new function every time we want to call the rootfinder. Instead we can call it directly with a lambda as an argument: fixed_point_iteration(lambda x: 1 / x + x / 2, 1, 1e-10).
2.2 Assignment
2.2.1 Bisection method
The bisection method was discussed in the lecture. Write a function bisect(f, left, right, tol) which takes four arguments: the function you want to find a root of, the left and right endpoints of the interval at the start, and the tolerance. The function should apply the bisection method until endpoints of the interval are within the given tolerance.
Use lambda functions to call bisect with different input functions. Use it to solve the following to some reasonable tolerance:
1. Find three values of x for which f(x) = x3 + 2×2 − x − 1 is zero.
2. Find the cube root of 7.
3. Find both solutions to tan−1 (x) = 3 − x2.
4. Find all solutions to log(x4 ) = x3 − 1.
For each of these you may want to first plot the function in order to select sensible initial estimates for the roots. Remember, for the bisection method you must supply initial guesses which bracket the root.
Explain what happens when you try to apply bisection method to the function f(x) = 1/x, with initial guesses −2 and 3.
Try to rewrite your function so that it uses a default tolerance of 10−6 if the user does not specify the tolerance, and check that this works.
2.2.2 Secant method
Write a secant method root-finding function,secant(f, x0, x1, tol), based on the material from the lecture. Check that it works by verifying the solutions to the problems above.
The function g(x) = cos(πx) has roots at x = n + 1/2, n ∈ N. Why do the three function calls: secant(g, 0.1, 1,0.0001), secant(g, 0.1, 1.2, 0.0001) and secant(g, 0.1, 1.4, 0.0001) find three different roots?
Practise the syntax for variable numbers of arguments by adapting the secant method program to accept an arbitrary number of coefficients of a polynomial.
2.2.3 Newton’s method
Write a Newton method root-finding function. Remember that for the Newton method you only need to give a single initial estimate, but you also need to supply the derivative. You can give both input function and its derivative as lambdas.
1. Use Newton’s method to verify that the function h(x) = 4 tan−1 x has a root at x = 0.
2. Try to find this root with a starting guess of x0 = 1, and then with an starting guess of x0 = 1.5. What goes wrong in the second case?
3. Now consider the function h(x) = x2 − 2 (which clearly has a root at x = √2) with an initial guess of x0 = 0. What is the problem?
2.2.4 Convergence
Here we will investigate, briefly, the convergence properties of each root-finding method, by finding the cube root of 1/2 with increasing accuracy.
1. For each method, alter your root-finding function to return a list (or better a “numpy array”) of all iterations, rather than the final one only.
2. We shall solve f(x) = x3 − 1/2 = 0 to find sequences of approximations of the cube root of 1/2. Check that all three methods work, given sensible starting parameters.
3. For each method, run the algorithm for a small value of the tolerance (10−15, say).
4. Plot the exact error against the number of steps for the three methods on the same graph. You should find that the secant and Newton methods are markedly quicker than the bisection method. In particular, the bisection method should appear roughly linear on a semilogy plot, while the others are faster than linear. Explain these results and find a method to infer the rate of convergence of the secant and Newton methods from the data.
