Suggestions for lab exercises.
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?
Write a function to calculate the volume of a cuboid with edge lengths $a, b, c$. Test your code on sample values such as
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
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.
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}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$.
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.
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!
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$.
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).
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 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).
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.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.
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.
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.
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$).
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.
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$.
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.
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$.
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 [35]:
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 [36]:
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 [37]:
n_combinations = 144*143*142*141/24
print("Number of combinations of 4 objects from 144 is {}".format(n_combinations))
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).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.
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.
This shows the sensitive dependence on initial conditions that is characteristic of chaotic behaviour.