Suppose we wish to evaluate exp x to a very high accuracy. As discussed in the lecture, floating-point arithmetic provides only limited precision, whi
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 <. 1 2 3 4 5 6 7 finally, test the methods floor, frac and to_float. check you understand what they do, how do it. 1.2.3 partial sums adapt function exp_taylor to compute first n terms of exp x in rational numbers. term is then given by rational(1, 1). probably need a number for n! each sum. use this give ten approximants e. 1.2.4 arbitrary accuracy download file expand.py from vle which contains code three methods: expand, expand_integer_part expand_fractional_part. place them class are used convert any base (like decimal), digits. example, [1]: 4).expand(10, 4) out[1]: '0.25 (base 10)' [2]: 7).expand(10, 30) out[2]: '0.142857142857142857142857142857 experiment with these conversion functions. 100th decimal e? many sum necessary establish this? 1.2.5 approximations π leonhard euler gave following expression π: write computes series, it confirm that 50th digit expansion 0. 1.2.6 advanced material srinivasa ramanujan was self-taught indian mathematician who had remarkable affinity he discovered extraordinary series π. where k!! double factorial, − 2)(k 4)· · × odd k. nth using class, zero. find digit? current record computation places 62.8 trillion (6.28 1013) digits (established august 2021). computed variation infinite 1910: class. will obtain an exact approximation √2, perhaps continued fractions. required π? chapter root-finding 2.1 tutorial one topic endlessly useful computational mathematics root-finding. is, f, value(s) f(x)="0." simple cases root can be found explicitly algebra, but general numerical technique needed. there several different ways approach problem, workshop we investigate well-known ones — bisection, secant, newton. since always finding roots functions, take opportunity review extend our knowledge python 2.1.1 functions now familiar defining python, formulation def function_name(argument_1, argument_2, argument_3, ...): return value point often block repeatedly, varying input parameters. help organise large program coherently, aid debugging. returns something; if no statement, none. 2.1.2 default argument values possible define or more arguments. test(a, b="1," c="True" ): print(a, b, c) produce test(4, 5, 6) test(3) true note called fewer arguments than defined allow. here thought as mandatory (in case a), optional c). not given, are, viewed order: [3]: 'false') false important evaluated only once. f(a, l="[]):" l.append(a) [4]: print(f(1)); print(f(2)); print(f(3)) [1] [1, 2] 2, 3] list emptied time called. 2.1.3 variable numbers occasionally might want accepts very example (and assuming better way accomplish this), suppose mean numbers, may change. has syntax achieve this. mean_of_n_numbers(*args): i args: +="i" len(args) call as: mean_of_n_numbers(1, 3, 4), inside function, args interpreted tuple. similarly, 4, 5] mean_of_n_numbers(*numbers) * unpacks lists variables into positional variables. course, could have been coded lists, *args standard method cope 2.1.4 lambdas another mechanism designed relatively simple, one-line procedures. alternative lambda. equivalent cube_plus_a(x, a): ** lambda x, a: a. does name so anonymous intended declarations typically once, forgotten. their simplicity manifested fact must single only. straightforward exactly meant “single expression”, statement something. assignment statements (such anything, cannot offer anything offered apart arguable simplicity, largely down personal preference. however, natural when passing such root-finder. fixed-point iteration until two successive iterates differ less user-supplied tolerance: 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.