Suggestions for lab exercises.
In [1]:
fifteen_factorial = 15*14*13*12*11*10*9*8*7*6*5*4*3*2*1
print(fifteen_factorial)
In [2]:
import math
print(math.factorial(15))
print("Result correct?", math.factorial(15) == fifteen_factorial)
Stirling's approximation gives that, for large enough $n$,
\begin{equation} n! \simeq \sqrt{2 \pi} n^{n + 1/2} e^{-n}. \end{equation}Using functions and constants from the math library, compare the results of $n!$ and Stirling's approximation for $n = 5, 10, 15, 20$. In what sense does the approximation improve?
In [3]:
print(math.factorial(5), math.sqrt(2*math.pi)*5**(5+0.5)*math.exp(-5))
print(math.factorial(10), math.sqrt(2*math.pi)*10**(10+0.5)*math.exp(-10))
print(math.factorial(15), math.sqrt(2*math.pi)*15**(15+0.5)*math.exp(-15))
print(math.factorial(20), math.sqrt(2*math.pi)*20**(20+0.5)*math.exp(-20))
print("Absolute differences:")
print(math.factorial(5) - math.sqrt(2*math.pi)*5**(5+0.5)*math.exp(-5))
print(math.factorial(10) - math.sqrt(2*math.pi)*10**(10+0.5)*math.exp(-10))
print(math.factorial(15) - math.sqrt(2*math.pi)*15**(15+0.5)*math.exp(-15))
print(math.factorial(20) - math.sqrt(2*math.pi)*20**(20+0.5)*math.exp(-20))
print("Relative differences:")
print((math.factorial(5) - math.sqrt(2*math.pi)*5**(5+0.5)*math.exp(-5)) / math.factorial(5))
print((math.factorial(10) - math.sqrt(2*math.pi)*10**(10+0.5)*math.exp(-10)) / math.factorial(10))
print((math.factorial(15) - math.sqrt(2*math.pi)*15**(15+0.5)*math.exp(-15)) / math.factorial(15))
print((math.factorial(20) - math.sqrt(2*math.pi)*20**(20+0.5)*math.exp(-20)) / math.factorial(20))
We see that the relative error decreases, whilst the absolute error grows (significantly).
Write a function to calculate the volume of a cuboid with edge lengths $a, b, c$. Test your code on sample values such as
In [4]:
def cuboid_volume(a, b, c):
"""
Compute the volume of a cuboid with edge lengths a, b, c.
Volume is abc. Only makes sense if all are non-negative.
Parameters
----------
a : float
Edge length 1
b : float
Edge length 2
c : float
Edge length 3
Returns
-------
volume : float
The volume a*b*c
"""
if (a < 0.0) or (b < 0.0) or (c < 0.0):
print("Negative edge length makes no sense!")
return 0
return a*b*c
In [5]:
print(cuboid_volume(1,1,1))
print(cuboid_volume(1,2,3.5))
print(cuboid_volume(0,1,1))
print(cuboid_volume(2,-1,1))
In later cases, after having covered exceptions, I would suggest raising a NotImplementedError for negative edge lengths.
Write a function to compute the time (in seconds) taken for an object to fall from a height $H$ (in metres) to the ground, using the formula
\begin{equation}
h(t) = \frac{1}{2} g t^2.
\end{equation}
Use the value of the acceleration due to gravity $g$ from scipy.constants.g. Test your code on sample values such as
In [6]:
def fall_time(H):
"""
Give the time in seconds for an object to fall to the ground
from H metres.
Parameters
----------
H : float
Starting height (metres)
Returns
-------
T : float
Fall time (seconds)
"""
from math import sqrt
from scipy.constants import g
if (H < 0):
print("Negative height makes no sense!")
return 0
return sqrt(2.0*H/g)
In [7]:
print(fall_time(1))
print(fall_time(10))
print(fall_time(0))
print(fall_time(-1))
In [8]:
def triangle_area(a, b, c):
"""
Compute the area of a triangle with edge lengths a, b, c.
Area is sqrt(s (s-a) (s-b) (s-c)).
s is (a+b+c)/2.
Only makes sense if all are non-negative.
Parameters
----------
a : float
Edge length 1
b : float
Edge length 2
c : float
Edge length 3
Returns
-------
area : float
The triangle area.
"""
from math import sqrt
if (a < 0.0) or (b < 0.0) or (c < 0.0):
print("Negative edge length makes no sense!")
return 0
s = 0.5 * (a + b + c)
return sqrt(s * (s-a) * (s-b) * (s-c))
In [9]:
print(triangle_area(1,1,1)) # Equilateral; answer sqrt(3)/4 ~ 0.433
print(triangle_area(3,4,5)) # Right triangle; answer 6
print(triangle_area(1,1,0)) # Not a triangle; answer 0
print(triangle_area(-1,1,1)) # Not a triangle; exception or 0.
Computers cannot, in principle, represent real numbers perfectly. This can lead to problems of accuracy. For example, if
\begin{equation} x = 1, \qquad y = 1 + 10^{-14} \sqrt{3} \end{equation}then it should be true that
\begin{equation} 10^{14} (y - x) = \sqrt{3}. \end{equation}Check how accurately this equation holds in Python and see what this implies about the accuracy of subtracting two numbers that are close together.
In [10]:
from math import sqrt
x = 1.0
y = 1.0 + 1e-14 * sqrt(3.0)
print("The calculation gives {}".format(1e14*(y-x)))
print("The result should be {}".format(sqrt(3.0)))
We see that the first three digits are correct. This isn't too surprising: we expect 16 digits of accuracy for a floating point number, but $x$ and $y$ are identical for the first 14 digits.
The standard quadratic formula gives the solutions to
\begin{equation} a x^2 + b x + c = 0 \end{equation}as
\begin{equation} x = \frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}. \end{equation}Show that, if $a = 10^{-n} = c$ and $b = 10^n$ then
\begin{equation} x = \frac{10^{2 n}}{2} \left( -1 \pm \sqrt{1 - 10^{-4n}} \right). \end{equation}Using the expansion (from Taylor's theorem)
\begin{equation} \sqrt{1 - 10^{-4 n}} \simeq 1 - \frac{10^{-4 n}}{2} + \dots, \qquad n \gg 1, \end{equation}show that
\begin{equation} x \simeq -10^{2 n} + \frac{10^{-2 n}}{4} \quad \text{and} \quad -\frac{10^{-2n}}{4}, \qquad n \gg 1. \end{equation}This is pen-and-paper work; each step should be re-arranging.
In [11]:
a = 1e-3
b = 1e3
c = a
formula1_n3_plus = (-b + sqrt(b**2 - 4.0*a*c))/(2.0*a)
formula1_n3_minus = (-b - sqrt(b**2 - 4.0*a*c))/(2.0*a)
formula2_n3_plus = (2.0*c)/(-b + sqrt(b**2 - 4.0*a*c))
formula2_n3_minus = (2.0*c)/(-b - sqrt(b**2 - 4.0*a*c))
print("For n=3, first formula, solutions are {} and {}.".format(formula1_n3_plus,
formula1_n3_minus))
print("For n=3, second formula, solutions are {} and {}.".format(formula2_n3_plus,
formula2_n3_minus))
a = 1e-4
b = 1e4
c = a
formula1_n4_plus = (-b + sqrt(b**2 - 4.0*a*c))/(2.0*a)
formula1_n4_minus = (-b - sqrt(b**2 - 4.0*a*c))/(2.0*a)
formula2_n4_plus = (2.0*c)/(-b + sqrt(b**2 - 4.0*a*c))
formula2_n4_minus = (2.0*c)/(-b - sqrt(b**2 - 4.0*a*c))
print("For n=4, first formula, solutions are {} and {}.".format(formula1_n4_plus,
formula1_n4_minus))
print("For n=4, second formula, solutions are {} and {}.".format(formula2_n4_plus,
formula2_n4_minus))
There is a difference in the fifth significant figure in both solutions in the first case, which gets to the third (arguably the second) significant figure in the second case. Comparing to the limiting solutions above, we see that the larger root is definitely more accurately captured with the first formula than the second (as the result should be bigger than $10^{-2n}$).
In the second case we have divided by a very small number to get the big number, which loses accuracy.
The standard definition of the derivative of a function is
\begin{equation} \left. \frac{\text{d} f}{\text{d} x} \right|_{x=X} = \lim_{\delta \to 0} \frac{f(X + \delta) - f(X)}{\delta}. \end{equation}We can approximate this by computing the result for a finite value of $\delta$:
\begin{equation} g(x, \delta) = \frac{f(x + \delta) - f(x)}{\delta}. \end{equation}Write a function that takes as inputs a function of one variable, $f(x)$, a location $X$, and a step length $\delta$, and returns the approximation to the derivative given by $g$.
In [12]:
def g(f, X, delta):
"""
Approximate the derivative of a given function at a point.
Parameters
----------
f : function
Function to be differentiated
X : real
Point at which the derivative is evaluated
delta : real
Step length
Returns
-------
g : real
Approximation to the derivative
"""
return (f(X+delta) - f(X)) / delta
In [13]:
from math import exp
for n in range(1, 8):
print("For n={}, the approx derivative is {}.".format(n, g(exp, 0.0, 10**(-2.0*n))))
We have a combination of floating point inaccuracies: in the numerator we have two terms that are nearly equal, leading to a very small number. We then divide two very small numbers. This is inherently inaccurate.
This does not mean that you can't calculate derivatives to high accuracy, but alternative approaches are definitely recommended.
This is a "simple" solution, but not efficient.
In [14]:
def isprime(n):
"""
Checks to see if an integer is prime.
Parameters
----------
n : integer
Number to check
Returns
-------
isprime : Boolean
If n is prime
"""
# No number less than 2 can be prime
if n < 2:
return False
# We only need to check for divisors up to sqrt(n)
for m in range(2, int(n**0.5)+1):
if n%m == 0:
return False
# If we've got this far, there are no divisors.
return True
In [15]:
for n in range(50):
if isprime(n):
print("Function says that {} is prime.".format(n))
We could do this many ways. This "elegant" solution says:
In [16]:
n = 2
while (not isprime(n)) or (isprime(2**n-1)):
n += 1
print("The first n such that 2^n-1 is not prime is {}.".format(n))
In [17]:
for n in range(2, 41):
if isprime(n) and isprime(2**n-1):
print("n={} is such that 2^n-1 is prime.".format(n))
Write a function to compute all prime factors of an integer $n$, including their multiplicities. Test it by printing the prime factors (without multiplicities) of $n = 17, \dots, 20$ and the multiplicities (without factors) of $n = 48$.
One effective solution is to return a dictionary, where the keys are the factors and the values are the multiplicities.
This solution uses the trick of immediately dividing $n$ by any divisor: this means we never have to check the divisor for being prime.
In [18]:
def prime_factors(n):
"""
Generate all the prime factors of n.
Parameters
----------
n : integer
Number to be checked
Returns
-------
factors : dict
Prime factors (keys) and multiplicities (values)
"""
factors = {}
m = 2
while m <= n:
if n%m == 0:
factors[m] = 1
n //= m
while n%m == 0:
factors[m] += 1
n //= m
m += 1
return factors
In [19]:
for n in range(17, 21):
print("Prime factors of {} are {}.".format(n, prime_factors(n).keys()))
print("Multiplicities of prime factors of 48 are {}.".format(prime_factors(48).values()))
Here we will do it directly.
In [20]:
def divisors(n):
"""
Generate all integer divisors of n.
Parameters
----------
n : integer
Number to be checked
Returns
-------
divs : list
All integer divisors, including 1.
"""
divs = [1]
m = 2
while m <= n/2:
if n%m == 0:
divs.append(m)
m += 1
return divs
In [21]:
for n in range(16, 21):
print("The divisors of {} are {}.".format(n, divisors(n)))
We can do this much more efficiently than the code below using packages such as numpy, but this is a "bare python" solution.
In [22]:
def isperfect(n):
"""
Check if a number is perfect.
Parameters
----------
n : integer
Number to check
Returns
-------
isperfect : Boolean
Whether it is perfect or not.
"""
divs = divisors(n)
sum_divs = 0
for d in divs:
sum_divs += d
return n == sum_divs
In [23]:
for n in range(2,10000):
if (isperfect(n)):
factors = prime_factors(n)
print("{} is perfect.\n"
"Divisors are {}.\n"
"Prime factors {} (multiplicities {}).".format(
n, divisors(n), factors.keys(), factors.values()))
In fact we did this above already:
Investigate the timeit function in python or IPython. Use this to measure how long your function takes to check that, if $k$ on the Mersenne list then $n = 2^{k-1} \times (2^k - 1)$ is a perfect number, using your functions. Stop increasing $k$ when the time takes too long!
You could waste considerable time on this, and on optimizing the functions above to work efficiently. It is not worth it, other than to show how rapidly the computation time can grow!
In [24]:
%timeit isperfect(2**(3-1)*(2**3-1))
In [25]:
%timeit isperfect(2**(5-1)*(2**5-1))
In [26]:
%timeit isperfect(2**(7-1)*(2**7-1))
In [27]:
%timeit isperfect(2**(13-1)*(2**13-1))
It's worth thinking about the operation counts of the various functions implemented here. The implementations are inefficient, but even in the best case you see how the number of operations (and hence computing time required) rapidly increases.
Partly taken from Newman's book, p 120.
The logistic map builds a sequence of numbers $\{ x_n \}$ using the relation
\begin{equation} x_{n+1} = r x_n \left( 1 - x_n \right), \end{equation}where $0 \le x_0 \le 1$.
In [28]:
def logistic(x0, r, N = 1000):
sequence = [x0]
xn = x0
for n in range(N):
xnew = r*xn*(1.0-xn)
sequence.append(xnew)
xn = xnew
return sequence
In [29]:
import numpy
from matplotlib import pyplot
%matplotlib inline
x0 = 0.5
N = 2000
sequence1 = logistic(x0, 1.5, N)
sequence2 = logistic(x0, 3.5, N)
pyplot.plot(sequence1[-100:], 'b-', label = r'$r=1.5$')
pyplot.plot(sequence2[-100:], 'k-', label = r'$r=3.5$')
pyplot.xlabel(r'$n$')
pyplot.ylabel(r'$x$')
pyplot.show()
This suggests that, for $r=1.5$, the sequence has settled down to a fixed point. In the $r=3.5$ case it seems to be moving between four points repeatedly.
Fix $x_0 = 0.5$. For each value of $r$ between $1$ and $4$, in steps of $0.01$, calculate the first 2,000 members of the sequence. Plot the last 1,000 members of the sequence on a plot where the $x$-axis is the value of $r$ and the $y$-axis is the values in the sequence. Do not plot lines - just plot markers (e.g., use the 'k.' plotting style).
In [30]:
import numpy
from matplotlib import pyplot
%matplotlib inline
r_values = numpy.linspace(1.0, 4.0, 401)
x0 = 0.5
N = 2000
for r in r_values:
sequence = logistic(x0, r, N)
pyplot.plot(r*numpy.ones_like(sequence[1000:]), sequence[1000:], 'k.')
pyplot.xlabel(r'$r$')
pyplot.ylabel(r'$x$')
pyplot.show()
For iterative maps such as the logistic map, one of three things can occur:
Using just your plot, or new plots from this data, work out approximate values of $r$ for which there is a transition from fixed points to limit cycles, from limit cycles of a given number of values to more values, and the transition to chaos.
The first transition is at $r \approx 3$, the next at $r \approx 3.45$, the next at $r \approx 3.55$. The transition to chaos appears to happen before $r=4$, but it's not obvious exactly where.
The Mandelbrot set is also generated from a sequence, $\{ z_n \}$, using the relation
\begin{equation} z_{n+1} = z_n^2 + c, \qquad z_0 = 0. \end{equation}The members of the sequence, and the constant $c$, are all complex. The point in the complex plane at $c$ is in the Mandelbrot set only if the $|z_n| < 2$ for all members of the sequence. In reality, checking the first 100 iterations is sufficient.
Note: the python notation for a complex number $x + \text{i} y$ is x + yj: that is, j is used to indicate $\sqrt{-1}$. If you know the values of x and y then x + yj constructs a complex number; if they are stored in variables you can use complex(x, y).
In [31]:
def in_Mandelbrot(c, n_iterations = 100):
z0 = 0.0 + 0j
in_set = True
n = 0
zn = z0
while in_set and (n < n_iterations):
n += 1
znew = zn**2 + c
in_set = abs(znew) < 2.0
zn = znew
return in_set
In [32]:
c_values = [0.0, 2+2j, 2-2j, -2+2j, -2-2j]
for c in c_values:
print("Is {} in the Mandelbrot set? {}.".format(c, in_Mandelbrot(c)))
In [33]:
import numpy
def grid_Mandelbrot(N):
x = numpy.linspace(-2.0, 2.0, N)
X, Y = numpy.meshgrid(x, x)
C = X + 1j*Y
grid = numpy.zeros((N, N), int)
for nx in range(N):
for ny in range(N):
grid[nx, ny] = int(in_Mandelbrot(C[nx, ny]))
return grid
In [34]:
from matplotlib import pyplot
%matplotlib inline
pyplot.imshow(grid_Mandelbrot(100))
Out[34]:
In [35]:
from math import log
def log_Mandelbrot(c, n_iterations = 100):
z0 = 0.0 + 0j
in_set = True
n = 0
zn = z0
while in_set and (n < n_iterations):
n += 1
znew = zn**2 + c
in_set = abs(znew) < 2.0
zn = znew
return log(n)
def log_grid_Mandelbrot(N):
x = numpy.linspace(-2.0, 2.0, N)
X, Y = numpy.meshgrid(x, x)
C = X + 1j*Y
grid = numpy.zeros((N, N), int)
for nx in range(N):
for ny in range(N):
grid[nx, ny] = log_Mandelbrot(C[nx, ny])
return grid
In [36]:
from matplotlib import pyplot
%matplotlib inline
pyplot.imshow(log_grid_Mandelbrot(100))
Out[36]:
This is a simple example:
In [37]:
pyplot.imshow(log_grid_Mandelbrot(1000)[600:800,400:600])
Out[37]:
An equivalence class is a relation that groups objects in a set into related subsets. For example, if we think of the integers modulo $7$, then $1$ is in the same equivalence class as $8$ (and $15$, and $22$, and so on), and $3$ is in the same equivalence class as $10$. We use the tilde $3 \sim 10$ to denote two objects within the same equivalence class.
Here, we are going to define the positive integers programmatically from equivalent sequences.
Define a python class Eqint. This should be
__repr__ function) to be the integer length of the sequence;__eq__ function) so that two eqints are equal if their sequences have same length.
In [38]:
class Eqint(object):
def __init__(self, sequence):
self.sequence = sequence
def __repr__(self):
return str(len(self.sequence))
def __eq__(self, other):
return len(self.sequence)==len(other.sequence)
Define a zero object from the empty list, and three one objects, from a single object list, tuple, and string. For example
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')
Check that none of the one objects equal the zero object, but all equal the other one objects. Print each object to check that the representation gives the integer length.
In [39]:
zero = Eqint([])
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')
print("Is zero equivalent to one? {}, {}, {}".format(zero == one_list,
zero == one_tuple,
zero == one_string))
print("Is one equivalent to one? {}, {}, {}.".format(one_list == one_tuple,
one_list == one_string,
one_tuple == one_string))
print(zero)
print(one_list)
print(one_tuple)
print(one_string)
Redefine the class by including an __add__ method that combines the two sequences. That is, if a and b are Eqints then a+b should return an Eqint defined from combining a and bs sequences.
Adding two different types of sequences (eg, a list to a tuple) does not work, so it is better to either iterate over the sequences, or to convert to a uniform type before adding.
In [40]:
class Eqint(object):
def __init__(self, sequence):
self.sequence = sequence
def __repr__(self):
return str(len(self.sequence))
def __eq__(self, other):
return len(self.sequence)==len(other.sequence)
def __add__(a, b):
return Eqint(tuple(a.sequence) + tuple(b.sequence))
In [41]:
zero = Eqint([])
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')
sum_eqint = zero + one_list + one_tuple + one_string
print("The sum is {}.".format(sum_eqint))
print("The internal sequence is {}.".format(sum_eqint.sequence))
We will sketch a construction of the positive integers from nothing.
positive_integers.Eqint called zero from the empty list. Append it to positive_integers.Eqint called next_integer from the Eqint defined by a copy of positive_integers (ie, use Eqint(list(positive_integers)). Append it to positive_integers.Use this procedure to define the Eqint equivalent to $10$. Print it, and its internal sequence, to check.
In [42]:
positive_integers = []
zero = Eqint([])
positive_integers.append(zero)
N = 10
for n in range(1,N+1):
positive_integers.append(Eqint(list(positive_integers)))
print("The 'final' Eqint is {}".format(positive_integers[-1]))
print("Its sequence is {}".format(positive_integers[-1].sequence))
print("That is, it contains all Eqints with length less than 10.")
Instead of working with floating point numbers, which are not "exact", we could work with the rational numbers $\mathbb{Q}$. A rational number $q \in \mathbb{Q}$ is defined by the numerator $n$ and denominator $d$ as $q = \frac{n}{d}$, where $n$ and $d$ are coprime (ie, have no common divisor other than $1$).
In [43]:
def normal_form(numerator, denominator):
from fractions import gcd
factor = gcd(numerator, denominator)
return numerator//factor, denominator//factor
In [44]:
print(normal_form(3, 2))
print(normal_form(15, 3))
print(normal_form(20, 42))
Define a class Rational that uses the normal_form function to store the rational number in the appropriate form. Define a __repr__ function that prints a string that looks like $\frac{n}{d}$ (hint: use len(str(number)) to find the number of digits of an integer). Test it on the cases above.
In [45]:
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
In [46]:
q1 = Rational(3, 2)
print(q1)
q2 = Rational(15, 3)
print(q2)
q3 = Rational(20, 42)
print(q3)
In [47]:
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
In [48]:
print(Rational(1,2) + Rational(1,3) + Rational(1,6))
In [49]:
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __mul__(a, b):
numerator = a.numerator * b.numerator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
In [50]:
print(Rational(1,3)*Rational(15,2)*Rational(2,5))
Overload the __rmul__ function so that you can multiply a rational by an integer. Check that $\frac{1}{2} \times 2 = 1$ and $\frac{1}{2} + (-1) \times \frac{1}{2} = 0$. Also overload the __sub__ function (using previous functions!) so that you can subtract rational numbers and check that $\frac{1}{2} - \frac{1}{2} = 0$.
In [51]:
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __mul__(a, b):
numerator = a.numerator * b.numerator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __rmul__(self, other):
numerator = self.numerator * other
return Rational(numerator, self.denominator)
def __sub__(a, b):
return a + (-1)*b
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
In [52]:
half = Rational(1,2)
print(2*half)
print(half+(-1)*half)
print(half-half)
In [53]:
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __mul__(a, b):
numerator = a.numerator * b.numerator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __rmul__(self, other):
numerator = self.numerator * other
return Rational(numerator, self.denominator)
def __sub__(a, b):
return a + (-1)*b
def __float__(a):
return float(a.numerator) / float(a.denominator)
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
In [54]:
print(float(Rational(1,2)))
print(float(Rational(1,3)))
print(float(Rational(1,11)))
In [55]:
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __mul__(a, b):
numerator = a.numerator * b.numerator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __rmul__(self, other):
numerator = self.numerator * other
return Rational(numerator, self.denominator)
def __sub__(a, b):
return a + (-1)*b
def __float__(a):
return float(a.numerator) / float(a.denominator)
def __lt__(a, b):
return a.numerator * b.denominator < a.denominator * b.numerator
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = '\n'+str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
In [56]:
q_list = [Rational(n//2, n) for n in range(2, 12)]
print(sorted(q_list))
The Wallis formula for $\pi$ is
\begin{equation} \pi = 2 \prod_{n=1}^{\infty} \frac{ (2 n)^2 }{(2 n - 1) (2 n + 1)}. \end{equation}We can define a partial product $\pi_N$ as
\begin{equation} \pi_N = 2 \prod_{n=1}^{N} \frac{ (2 n)^2 }{(2 n - 1) (2 n + 1)}, \end{equation}each of which are rational numbers.
Construct a list of the first 20 rational number approximations to $\pi$ and print them out. Print the sorted list to show that the approximations are always increasing. Then convert them to floating point numbers, construct a numpy array, and subtract this array from $\pi$ to see how accurate they are.
In [57]:
def wallis_rational(N):
"""
The partial product approximation to pi using the first N terms of Wallis' formula.
Parameters
----------
N : int
Number of terms in product
Returns
-------
partial : Rational
A rational number approximation to pi
"""
partial = Rational(2,1)
for n in range(1, N+1):
partial = partial * Rational((2*n)**2, (2*n-1)*(2*n+1))
return partial
In [58]:
pi_list = [wallis_rational(n) for n in range(1, 21)]
print(pi_list)
print(sorted(pi_list))
In [59]:
import numpy
print(numpy.pi-numpy.array(list(map(float, pi_list))))
A candidate for the shortest mathematical paper ever shows the following result:
\begin{equation} 27^5 + 84^5 + 110^5 + 133^5 = 144^5. \end{equation}This is interesting as
This is a counterexample to a conjecture by Euler ... that at least $n$ $n$th powers are required to sum to an $n$th power, $n > 2$.
In [60]:
lhs = 27**5 + 84**5 + 110**5 + 133**5
rhs = 144**5
print("Does the LHS {} equal the RHS {}? {}".format(lhs, rhs, lhs==rhs))
The more interesting statement in the paper is that
\begin{equation} 27^5 + 84^5 + 110^5 + 133^5 = 144^5. \end{equation}[is] the smallest instance in which four fifth powers sum to a fifth power.
Interpreting "the smallest instance" to mean the solution where the right hand side term (the largest integer) is the smallest, we want to use python to check this statement.
You may find the combinations function from the itertools package useful.
In [61]:
import numpy
import itertools
The combinations function returns all the combinations (ignoring order) of r elements from a given list. For example, take a list of length 6, [1, 2, 3, 4, 5, 6] and compute all the combinations of length 4:
In [62]:
input_list = numpy.arange(1, 7)
combinations = list(itertools.combinations(input_list, 4))
print(combinations)
We can already see that the number of terms to consider is large.
Note that we have used the list function to explicitly get a list of the combinations. The combinations function returns a generator, which can be used in a loop as if it were a list, without storing all elements of the list.
How fast does the number of combinations grow? The standard formula says that for a list of length $n$ there are
\begin{equation} \begin{pmatrix} n \\ k \end{pmatrix} = \frac{n!}{k! (n-k)!} \end{equation}combinations of length $k$. For $k=4$ as needed here we will have $n (n-1) (n-2) (n-3) / 24$ combinations. For $n=144$ we therefore have
In [63]:
n_combinations = 144*143*142*141/24
print("Number of combinations of 4 objects from 144 is {}".format(n_combinations))
In [64]:
from matplotlib import pyplot
%matplotlib inline
In [65]:
n = numpy.arange(5, 51)
N = numpy.zeros_like(n)
for i, n_c in enumerate(n):
combinations = list(itertools.combinations(numpy.arange(1,n_c+1), 4))
N[i] = len(combinations)
In [66]:
pyplot.figure(figsize=(12,6))
pyplot.loglog(n, N, linestyle='None', marker='x', color='k', label='Combinations')
pyplot.loglog(n, n**4, color='b', label=r'$n^4$')
pyplot.xlabel(r'$n$')
pyplot.ylabel(r'$N$')
pyplot.legend(loc='upper left')
pyplot.show()
With 17 million combinations to work with, we'll need to be a little careful how we compute.
One thing we could try is to loop through each possible "smallest instance" (the term on the right hand side) in increasing order. We then check all possible combinations of left hand sides.
This is computationally very expensive as we repeat a lot of calculations. We repeatedly recalculate combinations (a bad idea). We repeatedly recalculate the powers of the same number.
Instead, let us try creating the list of all combinations of powers once.
numpy array containing all integers in $1, \dots, 144$ to the fifth power. in keyword).
In [67]:
nmax=145
range_to_power = numpy.arange(1, nmax)**5
lhs_combinations = list(itertools.combinations(range_to_power, 4))
Then calculate the sums:
In [68]:
lhs_sums = []
for lhs_terms in lhs_combinations:
lhs_sums.append(numpy.sum(numpy.array(lhs_terms)))
Finally, loop through the sums and check to see if it matches any possible term on the RHS:
In [69]:
for i, lhs in enumerate(lhs_sums):
if lhs in range_to_power:
rhs_primitive = int(lhs**(0.2))
lhs_primitive = (numpy.array(lhs_combinations[i])**(0.2)).astype(int)
print("The LHS terms are {}.".format(lhs_primitive))
print("The RHS term is {}.".format(rhs_primitive))
The Lorenz system is a set of ordinary differential equations which can be written
\begin{equation} \frac{\text{d} \vec{v}}{\text{d} \vec{t}} = \vec{f}(\vec{v}) \end{equation}where the variables in the state vector $\vec{v}$ are
\begin{equation} \vec{v} = \begin{pmatrix} x(t) \\ y(t) \\ z(t) \end{pmatrix} \end{equation}and the function defining the ODE is
\begin{equation} \vec{f} = \begin{pmatrix} \sigma \left( y(t) - x(t) \right) \\ x(t) \left( \rho - z(t) \right) - y(t) \\ x(t) y(t) - \beta z(t) \end{pmatrix}. \end{equation}The parameters $\sigma, \rho, \beta$ are all real numbers.
In [70]:
def dvdt(v, t, sigma, rho, beta):
"""
Define the Lorenz system.
Parameters
----------
v : list
State vector
t : float
Time
sigma : float
Parameter
rho : float
Parameter
beta : float
Parameter
Returns
-------
dvdt : list
RHS defining the Lorenz system
"""
x, y, z = v
return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]
In [71]:
import numpy
from scipy.integrate import odeint
In [72]:
v0 = [1.0, 1.0, 1.0]
sigma = 10.0
beta = 8.0/3.0
t_values = numpy.linspace(0.0, 100.0, 5000)
rho_values = [13.0, 14.0, 15.0, 28.0]
v_values = []
for rho in rho_values:
params = (sigma, rho, beta)
v = odeint(dvdt, v0, t_values, args=params)
v_values.append(v)
In [73]:
%matplotlib inline
from matplotlib import pyplot
from mpl_toolkits.mplot3d.axes3d import Axes3D
In [74]:
fig = pyplot.figure(figsize=(12,6))
for i, v in enumerate(v_values):
ax = fig.add_subplot(2,2,i+1,projection='3d')
ax.plot(v[:,0], v[:,1], v[:,2])
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
ax.set_title(r"$\rho={}$".format(rho_values[i]))
pyplot.show()
Fix $\rho = 28$. Solve the Lorenz system twice, up to $t=40$, using the two different initial conditions $\vec{v}(0) = \vec{1}$ and $\vec{v}(0) = \vec{1} + \vec{10^{-5}}$.
Show four plots. Each plot should show the two solutions on the same axes, plotting $x, y$ and $z$. Each plot should show $10$ units of time, ie the first shows $t \in [0, 10]$, the second shows $t \in [10, 20]$, and so on.
In [75]:
t_values = numpy.linspace(0.0, 40.0, 4000)
rho = 28.0
params = (sigma, rho, beta)
v_values = []
v0_values = [[1.0,1.0,1.0],
[1.0+1e-5,1.0+1e-5,1.0+1e-5]]
for v0 in v0_values:
v = odeint(dvdt, v0, t_values, args=params)
v_values.append(v)
In [76]:
fig = pyplot.figure(figsize=(12,6))
line_colours = 'by'
for tstart in range(4):
ax = fig.add_subplot(2,2,tstart+1,projection='3d')
for i, v in enumerate(v_values):
ax.plot(v[tstart*1000:(tstart+1)*1000,0],
v[tstart*1000:(tstart+1)*1000,1],
v[tstart*1000:(tstart+1)*1000,2],
color=line_colours[i])
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
ax.set_title(r"$t \in [{},{}]$".format(tstart*10, (tstart+1)*10))
pyplot.show()
This shows the sensitive dependence on initial conditions that is characteristic of chaotic behaviour.
We are interested in the solution of
\begin{equation} \frac{\text{d} y}{\text{d} t} = e^{-t} - y^n, \qquad y(0) = 1, \end{equation}where $n > 1$ is an integer. The "minor" change from the above examples mean that sympy can only give the solution as a power series.
In [77]:
import sympy
sympy.init_printing()
In [78]:
y, t = sympy.symbols('y, t')
In [79]:
sympy.dsolve(sympy.diff(y(t), t) + y(t)**2 - sympy.exp(-t), y(t))
Out[79]:
In [80]:
for n in range(2, 11):
ode_solution = sympy.dsolve(sympy.diff(y(t), t) + y(t)**n - sympy.exp(-t), y(t),
ics = {y(0) : 1})
print(ode_solution)
In [81]:
%matplotlib inline
for n in range(2, 11):
ode_solution = sympy.dsolve(sympy.diff(y(t), t) + y(t)**n - sympy.exp(-t), y(t),
ics = {y(0) : 1})
sympy.plot(ode_solution.rhs.removeO(), (t, 0, 1));
A twin prime is a pair $(p_1, p_2)$ such that both $p_1$ and $p_2$ are prime and $p_2 = p_1 + 2$.
Write a generator that returns twin primes. You can use the generators above, and may want to look at the itertools module together with its recipes, particularly the pairwise recipe.
In [82]:
def all_primes(N):
"""
Return all primes less than or equal to N.
Parameters
----------
N : int
Maximum number
Returns
-------
prime : generator
Prime numbers
"""
primes = []
for n in range(2, N+1):
is_n_prime = True
for p in primes:
if n%p == 0:
is_n_prime = False
break
if is_n_prime:
primes.append(n)
yield n
Now we can generate pairs using the pairwise recipe:
In [83]:
from itertools import tee
def pair_primes(N):
"Generate consecutive prime pairs, using the itertools recipe"
a, b = tee(all_primes(N))
next(b, None)
return zip(a, b)
We could examine the results of the two primes directly. But an efficient solution is to use python's filter function. To do this, first define a function checking if the pair are twin primes:
In [84]:
def check_twin(pair):
"""
Take in a pair of integers, check if they differ by 2.
"""
p1, p2 = pair
return p2-p1 == 2
Then use the filter function to define another generator:
In [85]:
def twin_primes(N):
"""
Return all twin primes
"""
return filter(check_twin, pair_primes(N))
Now check by finding the twin primes with $N<20$:
In [86]:
for tp in twin_primes(20):
print(tp)
Again there are many solutions, but the itertools recipes has the quantify pattern. Looking ahead to exercise 3 we'll define:
In [87]:
def pi_N(N):
"""
Use the quantify pattern from itertools to count the number of twin primes.
"""
return sum(map(check_twin, pair_primes(N)))
In [88]:
pi_N(1000)
Out[88]:
We've now done all the hard work and can use the solutions above.
In [89]:
import numpy
from matplotlib import pyplot
%matplotlib inline
In [90]:
N = numpy.array([2**k for k in range(4, 17)])
twin_prime_fraction = numpy.array(list(map(pi_N, N))) / N
In [91]:
pyplot.semilogx(N, twin_prime_fraction)
pyplot.xlabel(r"$N$")
pyplot.ylabel(r"$\pi_N / N$")
pyplot.show()
For those that have checked Wikipedia, you'll see Brun's theorem which suggests a specific scaling, that $\pi_N$ is bounded by $C N / \log(N)^2$. Checking this numerically on this data:
In [92]:
pyplot.semilogx(N, twin_prime_fraction * numpy.log(N)**2)
pyplot.xlabel(r"$N$")
pyplot.ylabel(r"$\pi_N \times \log(N)^2 / N$")
pyplot.show()
In the section on classes we defined a Monomial class to represent a polynomial with leading coefficient $1$. As the $N+1$ monomials $1, x, x^2, \dots, x^N$ form a basis for the vector space of polynomials of order $N$, $\mathbb{P}^N$, we can use the Monomial class to return this basis.
Define a generator that will iterate through this basis of $\mathbb{P}^N$ and test it on $\mathbb{P}^3$.
Again we first take the definition of the crucial class from the notes.
In [93]:
class Polynomial(object):
"""Representing a polynomial."""
explanation = "I am a polynomial"
def __init__(self, roots, leading_term):
self.roots = roots
self.leading_term = leading_term
self.order = len(roots)
def __repr__(self):
string = str(self.leading_term)
for root in self.roots:
if root == 0:
string = string + "x"
elif root > 0:
string = string + "(x - {})".format(root)
else:
string = string + "(x + {})".format(-root)
return string
def __mul__(self, other):
roots = self.roots + other.roots
leading_term = self.leading_term * other.leading_term
return Polynomial(roots, leading_term)
def explain_to(self, caller):
print("Hello, {}. {}.".format(caller,self.explanation))
print("My roots are {}.".format(self.roots))
return None
In [94]:
class Monomial(Polynomial):
"""Representing a monomial, which is a polynomial with leading term 1."""
explanation = "I am a monomial"
def __init__(self, roots):
Polynomial.__init__(self, roots, 1)
def __repr__(self):
string = ""
for root in self.roots:
if root == 0:
string = string + "x"
elif root > 0:
string = string + "(x - {})".format(root)
else:
string = string + "(x + {})".format(-root)
return string
Now we can define the first basis:
In [95]:
def basis_pN(N):
"""
A generator for the simplest basis of P^N.
"""
for n in range(N+1):
yield Monomial(n*[0])
Then test it on $\mathbb{P}^N$:
In [96]:
for poly in basis_pN(3):
print(poly)
This looks horrible, but is correct. To really make this look good, we need to improve the output. If we use
In [97]:
class Monomial(Polynomial):
"""Representing a monomial, which is a polynomial with leading term 1."""
explanation = "I am a monomial"
def __init__(self, roots):
Polynomial.__init__(self, roots, 1)
def __repr__(self):
if len(self.roots):
string = ""
n_zero_roots = len(self.roots) - numpy.count_nonzero(self.roots)
if n_zero_roots == 1:
string = "x"
elif n_zero_roots > 1:
string = "x^{}".format(n_zero_roots)
else: # Monomial degree 0.
string = "1"
for root in self.roots:
if root > 0:
string = string + "(x - {})".format(root)
elif root < 0:
string = string + "(x + {})".format(-root)
return string
then we can deal with the uglier cases, and re-running the test we get
In [98]:
for poly in basis_pN(3):
print(poly)
An even better solution would be to use the numpy.unique function as in this stackoverflow answer (the second one!) to get the frequency of all the roots.
An alternative basis is given by the monomials
\begin{align} p_0(x) &= 1, \\ p_1(x) &= 1-x, \\ p_2(x) &= (1-x)(2-x), \\ \dots & \quad \dots, \\ p_N(x) &= \prod_{n=1}^N (n-x). \end{align}Define a generator that will iterate through this basis of $\mathbb{P}^N$ and test it on $\mathbb{P}^4$.
In [99]:
def basis_pN_variant(N):
"""
A generator for the 'sum' basis of P^N.
"""
for n in range(N+1):
yield Monomial(range(n+1))
In [100]:
for poly in basis_pN_variant(4):
print(poly)
I am too lazy to work back through the definitions and flip all the signs; it should be clear how to do this!
Hopefully by now you'll be aware of how useful itertools is!
In [101]:
from itertools import product
In [102]:
def basis_product():
"""
Basis of the product space
"""
yield from product(basis_pN(3), basis_pN_variant(4))
In [103]:
for p1, p2 in basis_product():
print("Basis element is ({}) X ({}).".format(p1, p2))
I've cheated here as I haven't introduced the yield from syntax (which returns an iterator from a generator). We could write this out instead as
In [104]:
def basis_product_long_form():
"""
Basis of the product space (without using yield_from)
"""
prod = product(basis_pN(3), basis_pN_variant(4))
yield next(prod)
In [105]:
for p1, p2 in basis_product():
print("Basis element is ({}) X ({}).".format(p1, p2))