Diary

The Fibonacci sequence is one of the most well known mathematical sequences and is the most basic example of recurrence relations. Each number in the sequence consists of the sum of the previous two numbers and the starting two numbers are \(0\) and \(1\). It goes something like this:

\(1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144 \text{ and so on forever…}\)

We will go through a few different approaches for the optimal solution:

1. Using simple recursion
2. Using cache with recursion
3. Using Binet’s formula

Using simple recursion

This is a very simple and easy way to return the nth Fibonacci number in Python:

def recursiveFib(n):
    if n == 1 or n == 2:
        return 1

    return recursiveFib(n - 1) + recursiveFib(n - 2)

It uses recursion to repeatedly call itself multiple times calculating the previous number and using it to work out the next. But that is also its downside, as the function extremely inefficient and resource intensive, this is because at every stage it calculates the previous two numbers, and the previous two numbers of those number etc. Soon you reach a point where it takes too long to calculate the next number – for example on my computer it took me 1.61 seconds to calculate the 35th number. This will obviously be extremely slow, and basically impossible, to calculate higher values.

Using cache with recursion

Since we constantly calculate the previous two numbers, we can harness the power of caching to store the number, so we do not need to calculate them anymore. The built-in functools module allows us to use a least recently used cache; a type of cache which organises the items in order of use. This can speed up the process by a huge amount.

# procedura ricorsiva con cache
@lru_cache()
def recursiveFibCached(n):
    if n == 1 or n == 2:
        return 1

    return recursiveFibCached(n - 1) + recursiveFibCached (n - 2)

Firstly, we need to import the lru_cache decorator from the functools module and place it before our function. We can supply a maxsize value to tell the cache how many items to store, but by default it is 128 and that works perfectly fine.

Using Binet’s formula

Binet’s formula is a formula which can be used to calculate the nth term of the Fibonacci sequence, which is exactly what we want to do, and is named after it was derived by the french mathematician Jacques Philippe Marie Binet. The formula is shown below:

\(F_n = \frac{\phi^n-(-\phi)^{-n}}{\sqrt{5}}\)

where:

\(\phi = \frac{1+\sqrt{5}}{2}\)

We can convert thing formula into Python and start using it straight away:

def formulaFib(n):
    root_5 = 5 ** 0.5
    phi = ((1 + root_5) / 2)

    a = ((phi ** n) - ((-phi) ** (-n))) / root_5

    return round(a)

Note for the Python implementation we need to return the rounded version of the number we calculate, this is because when we calculate a large number Python will return to us a number where there can be 20+ trailing 9’s after the decimal place.

This is all well and good since now we do not have any loops and can instantly compute the answer, right? Well, there is a small catch to this method. If we try to calculate anything above the 1475th number, we will run into an error; OverflowError: (34, 'Numerical result out of range'). This is due to the way floats are implemented in Python and they can only have a certain maximum value, which we exceed using this method.

However, a fix to this is very easy to implement. We can use a built-in module called decimal to create a decimal object with a much higher precision and size to use in our equation:

import decimal

def formulaFibWithDecimal(n):
    decimal.getcontext().prec = 10000

    root_5 = decimal.Decimal(5).sqrt()
    phi = ((1 + root_5) / 2)

    a = ((phi ** n) - ((-phi) ** (-n))) / root_5

    return round(a)

In this new function we set the precision value to 10’000 digits long and we convert our root 5 value into a decimal object value and use that in our equation. This allows us to easily calculate the 10’000th number in the sequence in an astonishing 41.9 seconds, a huge improvement over all our previous methods.

The Python ecosystem gives us the option of visualizing relationships between numbers via matplotlib, a plotting library (i.e., not part of core Python) which can produce quality figures. Inside the matplotlib package is the matplotlib.pyplot module, which is used to produce figures in a MATLAB-like environment.

Here a simple example code for the matplotlib.pyplot module.

import matplotlib.pyplot as plt

def plotex(cxs,cys,dxs,dys):
    plt.xlabel('$x$', fontsize=10)
    plt.ylabel('$f(x)$', fontsize=10)
    plt.plot(cxs, cys, 'r-', label='quadratic: $y=x^2$')
    plt.plot(dxs, dys, 'b--^', label='other function: $y=x^{1.8}-1/2$')
    plt.legend()
    plt.show()

x1 = [0.1*i for i in range(60)]
y1 = [x**2 for x in x1]

x2 = [i for i in range(7)]
y2 = [x**1.8 - 0.5 for x in x2]

plotex(x1, y1, x2, y2)

We start by importing matplotlib.pyplot in the (standard) way which allows us to use it below without repeated typing of unnecessary characters.

We then define a function, plotex(), that takes care of the plotting, whereas the main program simply introduces four list comprehensions and then calls our function. If you’re still a beginner, you may be wondering why we defined a Python function in this code. An important design principle in computer science goes by the name of separation of concerns (or sometimes information hiding or encapsulation): each aspect of the program should be handled separately. In our case, this means that each component of our task should be handled in a separate function.

Let’s discuss this function in more detail. Its parameters are (meant to be) four lists, namely two pairs of \(x_i\) and \(y_i\) values. The function body starts by using xlabel() and ylabel() to provide labels for the x and y axes. It then creates individual curves/sets of points by using matplotlib’s function plot(), passing in the x-axis values as the first argument and the y-axis values as the second argument. The third positional argument to plot() is the format string: this corresponds to the color and point/line type. In the first case, we used r for red and - for a solid line.

The fourth argument to plot() is a keyword argument containing the label corresponding to the curve. In the second call to plot() we pass in a different format string and label (and, obviously, different lists); observe that we used two style options in the format string: -- to denote a dashed line and ˆ to denote the points with a triangle marker. The function concludes by calling legend(), which is responsible for making the legend appear, and show(), which makes the plot actually appear on our screen.

We could fine-tune almost all aspects of our plots, including basic things like line width, font size, and so on. For example, we could get LaTeX-like equations by putting dollar signs inside our string, e.g., ‘$x i$’ appears as \(x_i\).

Electric Field of a Distribution of Point Charges

Very briefly, let us recall Coulomb’s law: the force on a test charge Q located at point P (at the position r), coming from a single point charge q₀ located at r₀ is given by:

\(\mathbf{F}_{0}=k \frac{q_{0} Q}{\left(\mathbf{r}-\mathbf{r}_{0}\right)^{2}} \frac{\mathbf{r}-\mathbf{r}_{0}}{\left|\mathbf{r}-\mathbf{r}_{0}\right|}\)

where Coulomb’s constant is \(k=1 /\left(4 \pi \epsilon_{0}\right)\) in SI units (and \(\epsilon_0\) is the permittivity of free space).

The force is proportional to the product of the two charges, inversely proportional to the square of the distance between the two charges, and points along the line from charge q₀ to charge Q. The electric field is then the ratio of the force F₀ with the test charge Q in the limit where the magnitude of the test charge goes to zero. In practice, this given us:

\(\mathbf{E}_{0}(\mathbf{r})=k q_{0} \frac{\mathbf{r}-\mathbf{r}_{0}}{\left|\mathbf{r}-\mathbf{r}_{0}\right|^{3}}\)

where we cancelled out the Q and also took the opportunity to combine the two denominators. This is the electric field at the location r due to the point charge q₀ at r₀.

If we were faced with more than one point charge, we could apply the principle of superposition: the total force on Q is made up of the vector sum of the individual forces acting on Q. As a result, if we were dealing with the n point charges q₀, q₁,…, q_{n −1} located at r₀, r₁,…, r_n−1 (respectively) then the electric field at the location r is:

\(\mathbf{E}(\mathbf{r})=\sum_{i=0}^{n-1} \mathbf{E}_{i}(\mathbf{r})=\sum_{i=0}^{n-1} k q_{i} \frac{\mathbf{r}-\mathbf{r}_{i}}{\left|\mathbf{r}-\mathbf{r}_{i}\right|^{3}}\)

Note that you can consider this total electric field at any point in space, r. Note, also, that the electric field is a vector quantity: at any point in space this E has a magnitude and a direction. One way of visualizing vector fields consists of drawing field lines, namely imaginary curves that help us keep track of the direction of the field. More specifically, the tangent of a field line at a given point gives us the direction of the electric field at that point. Field lines do not cross; they start at positive charges (“sources”) and end at negative charges (“sinks”).

We will plot the electric field lines in Python; while more sophisticated ways of visualizing a vector field exist (e.g., line integral convolution), what we describe below should be enough to give you a qualitative feel for things.

We are faced with two tasks: first, we need to produce the electric field (vector) at several points near the charges and, second, we need to plot the field lines in such a way that we can physically interpret what is happening.

Code below is a Python implementation, where Coulomb’s constant is divided out for simplicity.

import numpy as np
import matplotlib.pyplot as plt
from math import sqrt
from copy import deepcopy

def makefield(xs, ys):
    qtopos = {1: (-1,0), -1: (1,0)}
    n = len(xs)
    Exs = [[0. for k in range(n)] for j in range(n)]
    Eys = deepcopy(Exs)
    for j,x in enumerate(xs):
        for k,y in enumerate(ys):
            for q,pos in qtopos.items():
                posx, posy = pos
                R = sqrt((x - posx)**2 + (y - posy)**2)
                Exs[k][j] += q*(x - posx)/R**3
                Eys[k][j] += q*(y - posy)/R**3
    return Exs, Eys

def plotfield(boxl,n):
    xs = [-boxl + i*2*boxl/(n-1) for i in range(n)]
    ys = xs[:]
    Exs, Eys = makefield(xs, ys)
    xs=np.array(xs); ys=np.array(ys)
    Exs=np.array(Exs); Eys=np.array(Eys)
    plt.streamplot(xs, ys, Exs, Eys, density=1.5, color='m')
    plt.xlabel('$x$')
    plt.ylabel('$y$')
    plt.show()

plotfield(2.,20)

We start by importing numpy and matplotlib, since the heavy lifting will be done by the function streamplot(), which expects NumPy arrays as input. We also import the square root and the deepcopy() function, which can create a distinct list-of-lists.

The function makefield() takes in two lists, xs and ys, corresponding to the coordinates at which we wish to evaluate the electric field (x and y together make up r). We also need some way of storing the r_i at which the point charges are located. We have opted to store these in a dictionary, which maps from charge q_i to position r_i. For each position r we need to evaluate E(r): in two dimensions, this is made up of E_x (r) and E_y (r), namely the two Cartesian components of the total electric field. Focusing on only one of these for the moment, say E_x (r), we realize that we need to store its value for any possible r, i.e., for any possible x and y values. We decide to use a list-of-lists, produced by a nested list comprehension. We then create another list-of-lists, for Ey(r). We need to map out (i.e., store) the value of the x and y components of the total electric field, at all the desired values of the vector r, namely, on a two-dimensional grid made up of xs and ys. This entails computing the electric field (contribution from a given point charge q_i ) at all possible y’s for a given x, and then iterating over all possible x’s. We also need to iterate over our point charges q_i and their locations r_i; we do this by saying for q, pos in qtopos.items(): at which point we unpack pos into posx and posy. We thus end up with three nested loops: one over possible x values, one over possible y values, and one over i. All three of these are written idiomatically, employing items() and enumerate().

Our second function, plotfield(), is where we build our two-dimensional grid for the xs and ys. We take in as parameters the length L and the number of points n we wish to use in each dimension and create our xs using a list comprehension; all we’re doing is picking x’s from −L to L. We then create a copy of xs and name it ys. After this,
we call our very own makefield() to produce the two lists-of-lists containing E_x (r) and E_y (r) for many different choices of r.

The result of running this code is shown in the following figure.

• We have a giant vase, which is infinitely big, that can hold an infinite amount of balls.
• We also have an infinite amount of balls that we can insert into the vase.
• The balls are numbered in order 1, 2, 3, …
• We are adding balls at a rate of 10 balls per each task.
• Tasks are performed at decreasing timing, with first task starting at 11:59:00 AM:
  • The first task would is completed in half a min (1/2 min).
  • The second task, would be completed in half the amount, so a quarter of a min (1/4 min).
  • The third one, in half the last one’s amount, so one eigth of a min (1/8 min).
  • and so on and so forth…until, the last task to be completed by noon.

So basically, it drills down to performing an infinite amount of tasks, in an finite amount of time, which is the definition of what we refer as a supertask.

We start with task 1, where we add balls 1 to 10, and remove the first, which is ball 1. At task 2, we add balls 11-20, and remove the remaining first, which is ball 2. At task 3, we add balls 21-30, and remove the remaining first, ball 3… and so on and so forth.

Now the question is: how many balls would be in the vase at noon?

Let us look into what happened on the actual the steps. At step 1, we removed ball 1, and then step 2, removed ball 2, step 3, removed ball 3,… so at step n, we would remove ball n, and hence all the balls will be removed from the vase by the end of the tasks. So, the answer to our problem is none.

Apparently, the order in which we add balls, and remove them, affects how many balls are in the vase at noon. Let us look at an alternative scenario.

• At task 1, we add balls 1 to 10, and remove the last, which is ball 10
• At task 2, we add balls 11-20, remove the last, which is ball 20
• At task 3, we add balls 21-30, remove the last, ball 30,
• And so on and so forth…

So again, how many balls would be in the vase at noon?
This time, the answer is different, and it is Infinity. This is actually due to the fact that only balls which are multiples of ten are the ones being removed, so the vase will not contain balls 10, 20, 30,… yet all the other balls will still be there, hence the infinite amount.