FIGlet – as you can read here – is a program for making large letters out of ordinary text.

A very quick note on how to add FontArt to Python’s output view.I will give a few tips on how to add FontArt to the python output display.

The method is quite easy, the first thing we have to do is install the PyFiglet library. PyFiglet converts ASCII text to ASCII art font.

pip install pyfiglet

After the installation is complete, then open the code editor and type the following code:

# You can get the type of font itself at the following link.

import pyfiglet
result = pyfiglet.figlet_format('wow', font='banner3')

The result of the code is like this:

In addition, we can also change the displayed font, by adding a font parameter to the .figlet_format function:

# You can get the type of font itself at the following link.

import pyfiglet
result = pyfiglet.figlet_format('wow', font='colossal')

The result of the code is like this:

Below, a complete code for the project:

import pyfiglet
import argparse

# python
# python --text "Ciao"
# python --text "Ciao" --font banner3-d

# You can see the font list here:

parser = argparse.ArgumentParser()
parser.add_argument('-t', '--text', help='text')
parser.add_argument('-f', '--font', help='font')

text = 'This is sample text'

args = parser.parse_args()
if args.text:
    text = args.text

if args.font:
    result = pyfiglet.figlet_format(text, font=args.font)
    result = pyfiglet.figlet_format(text)


The Fibonacci sequence…


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.

processing time of the recoursiveFibo routine

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
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 cache with recursion

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}}\)


\(\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)
Using Binet’s formula

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.

Basic Plotting with matplotlib

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$')

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)
A matplotlib.pyplot module example

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.

Some of the most important line styles/markers in matplotlib

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 q0 located at r0 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 q0 to charge Q. The electric field is then the ratio of the force F0 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 q0 at r0.

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 q0, q1,…, qn −1 located at r0, r1,…, rn−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')


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 ri at which the point charges are located. We have opted to store these in a dictionary, which maps from charge qi to position ri. For each position r we need to evaluate E(r): in two dimensions, this is made up of Ex (r) and Ey (r), namely the two Cartesian components of the total electric field. Focusing on only one of these for the moment, say Ex (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 qi ) 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 qi and their locations ri; 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 Ex (r) and Ey (r) for many different choices of r.

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

Visualizing the electric fields resulting from two point charges.

1-D Linear Convection

The 1-D Linear Convection equation is the simplest, most basic model that can be used to learn something about CFD.

\(\frac{{\partial u}}{{\partial t}} + c\frac{{\partial u}}{{\partial x}} = 0\)

With given initial conditions (understood as a wave), the equation represents the propagation of that initial wave with speed \(c\), without change of shape.

Let the initial condition be \(u(x,0)=u_0(x)\). Then the exact solution of the equation is \(u(x,t)=u_0(x-ct)\).

We discretize this equation, in both space and time, using the Forward Difference scheme for the time derivative and the Backward Difference scheme for the space derivative.

Our discrete equation is:

\(\frac{{u_i^{n + 1} – u_i^n}}{{\Delta t}} + c\frac{{u_i^n – u_{i – 1}^n}}{{\Delta x}} = 0\)

Now let’s try implementing this in Python.

We’ll start by importing a few libraries to help us out.

  • numpy is a library that provides a bunch of useful matrix operations akin to MATLAB;
  • matplotlib is a 2D plotting library that we will use to plot our results;
  • time and sys provide basic timing functions that we’ll use to slow down animations for viewing.
# Remember: comments in python are denoted by the pound sign
import numpy as np                      #here we load numpy
from matplotlib import pyplot as plt      #here we load matplotlib
import time, sys                   #and load some utilities

nx = 41         # spatial discretization: number of grid points
dx = 2 / (nx-1) # the distance between any pair of adjacent grid points
nt = 25         # nt is the number of timesteps we want to calculate
dt = .025       # dt is the amount of time each timestep covers (delta t)
c = 1           # assume wavespeed of c = 1

We have define an evenly spaced grid of points within a spatial domain that is 2 units of length wide, i.e., \(x_i \in [0,2]\). We have define a variable nx, which represent the number of grid points we want and dx is the distance between any pair of adjacent grid points.

We also need to set up our initial conditions. The initial velocity \(u_0\) is given as \(u=2\) in the interval \(0.5 \leq x \leq 1\) and \(u=1\) everywhere else in \((0,2)\) (i.e., a hat function). Here, we use the function ones() defining a numpy array which is nx elements long with every value equal to 1.

u = np.ones(nx)                      #numpy function ones()
u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s

x = np.linspace(0,2,nx)
plt.plot(x,u, 'red', label='Initial shape')

Now it’s time to implement the discretization of the convection equation using a finite-difference scheme.

We’ll store the result in a new (temporary) array un, which will be the solution \(u\) for the next time-step. We will repeat this operation for as many time-steps as we specify and then we can see how far the wave has convected.

We first initialize our placeholder array un to hold the values we calculate for the \(n+1\) timestep, using once again the NumPy function ones().

Then, we may think we have two iterative operations: one in space and one in time (we’ll learn differently later), so we’ll start by nesting one loop inside the other. Note the use of the nifty range() function. When we write: for i in range(1,nx) we will iterate through the u array, but we’ll be skipping the first element (the zero-th element).

un = np.ones(nx) #initialize a temporary array

for n in range(nt):        #loop for values of n from 0 to nt, so it will run nt times
    un = u.copy()          ##copy the existing values of u into un
    for i in range(1, nx): ## you can try commenting this line and...
        u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])

plt.plot(x,u, 'green', label='Final shape')

pplt.title('1-D Linear Convection : nx=41, nt=25, dt=.025 and c=1')


OK! So our hat function has definitely moved to the right, but it’s no longer a hat. What’s going on?

Basic data manipulation in Python…

In this post, we will deal with data from ECDC and we will explain basic data manipulation in Python with the Pandas package.

In our day, data is everywhere in enormous size and depth. Data science is an emerging field that penetrates every aspect of our life and, lately, it has proved to be an extraordinary weapon for predicting infections from Covid-19 and organizing strategies to limit the damage.

To import Pandas and Matplotlib packages we code:

import pandas as pd
import matplotlib.pyplot as plt

We download the excel file locally from ECDC site and open it using the read_excel function of Pandas library. We have named the file as data.xls in our case.

df=pd.read_excel("data.xlsx", engine="openpyxl")

We can first explore the data and the columns of the dataframe df:

We observe the columns of the dataframe — in our case, we will use the columns: dateRep, cases and deaths. Additionally, the name of the country is stored in column countriesAndTerritories.

We next select ‘Italy’ as the country under study. A new column is created named DateTime of type datetime where we store the day. In the following, we create a new dataframe with the name df_italia which is the same as the dataframe df_italia_sorted but it is sorted according to the column DateTime


#We sort according to DateTime


We are interested in data after the month of April (i.e., May, June, July, August, … etc) so we choose to filter using the column month and create a new dataframe df_italia_selected.

Since the data in columns cases and deaths may have great variation, it is practical in order to understand the trend to use a moving average. We choose a moving average of seven days and we create two new columns (Moving Average Cases and Moving Average Deaths) where we store the average values of cases and deaths.

#Calculate moving average

df_italia_selected['Moving Average Cases']=df_italia_selected.cases.rolling(7,min_periods=1).mean()
df_italia_selected['Moving Average Deaths']=df_italia_selected.deaths.rolling(7,min_periods=1).mean()

We now plot the cases and deaths as functions of time. We choose the red color for cases and blue for deaths. It is useful to plot cases and deaths in the same figure with common x-axis in order to understand possible connection and relation. So, we use the subplots function and first create figure fig and axis ax1 (this will be the axis for the cases and it will be the left axis). We then create ax2 using twinx function. The values for deaths will be our right axis. A dashed line is used for the average values.


fig, ax1=plt.subplots()

ax1.plot(df_italia_selected['DateTime'], df_italia_selected['cases'], color=color1)

ax1.plot(df_italia_selected['DateTime'], df_italia_selected['Moving Average Cases'], color=color1,linestyle='dashed')



locs, labels=plt.xticks()

ax2=ax1.twinx() #instantiate a second axes that shares the same x-axis

ax2.plot(df_italia_selected['DateTime'], df_italia_selected['deaths'], color=color2)
ax2.plot(df_italia_selected['DateTime'], df_italia_selected['Moving Average Deaths'], color=color2,linestyle='dashed')


fig.tight_layout() #otherwise the right y-label is slightly clipped

The figure below is the program output.

Cases and deaths as a function of data for Italy

A spiral of first 10k prime numbers…

A prime number is a natural number greater than 1 that is not a product of two smaller natural numbers.

Euclid proved the infinity of primes by contradiction.

  • Assume there are a finite number, n , of primes , the largest being \(p_n\);
  • Consider the number that is the product of these, plus one: \(N = \prod\limits_{i = 1}^n {{p_i}} + 1\)
  • By construction, \(N\) is not divisible by any of the \(p_i\).
  • Hence it is either prime itself, or divisible by another prime greater than \(p_n\), contradicting the assumption.

Euclid’s proof is the easiest to understand, especially if you’re not well-versed in mathematics, but now I am going to introduce another proof of the infinitude of primes, which is shorter and I also find it to be the most appealing, at least in my eyes.

In mathematics, the prime-counting function is the function counting the number of prime numbers less than or equal to some real number \(x\). E.g. for \(x=10.124\), we have that 4 number of prime number: 2, 3, 5, 7.

Let \(\pi(x)\) be a prime-counting function.

The Prime Number Theorem suggests that:

\(\mathop {\lim }\limits_{x \to \infty } \frac{{\pi \left( x \right)}}{{x/\ln (x)}} = 1\)

In other words, we can say that \(\pi(x)\) “behaves” like \(x/\ln (x)\) when \(x\) is approaching infinity — there’s an asymptotic ‘equality’.

So, there are infinite primes.

For Python code, we use SymPy – a Python library for symbolic mathematics. In particular, we’ll use primerange(a, b) for generate all prime numbers in the range [a, b).

The code:

import sympy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as ani'dark_background')'default')

def convert_to_polar_coordinates(num):
    return num*np.cos(num), num*np.sin(num)

lim = 10000

fig, ax = plt.subplots(figsize = (8, 8))
primes_numbers = sympy.primerange(0, lim)
primes_numbers = np.array(list(primes_numbers))
x, y = convert_to_polar_coordinates(primes_numbers)

def init():
    ax.plot([], [])

def plot_in_polar(i):
    ax.plot(x[:i], y[:i], linestyle='', marker='o', markersize=0.75, c='#FC0000')
    ax.set_xlim(-lim, lim)
    ax.set_ylim(-lim, lim)


animator = ani.FuncAnimation(fig=fig, func=plot_in_polar, interval=1, frames=len(primes_numbers))"prime.gif")
A spiral of first 10k prime numbers in a polar coordinate system.

Spiral of Theodorus…

The spiral of Theodorus is a spiral composed of right triangles. Hundreds of years ago, Theodorus of Cyrene constructed continuous right triangles and got a beautiful spiral. He used that spiral to prove that all non-square integers from 3–17 are irrational.

How would you plot this spiral? At each step, you need to draw a segment of length 1, perpendicular to the hypotenuse of the previous triangle. There are two perpendicular directions, and you want to choose the one that moves counterclockwise.

the spiral of Theodorus

If we step outside the xy plane, we can compute the cross product of the unit vector in the z direction with the vector (x, y). The cross product will be perpendicular to both, and by the right-hand rule, it will point in the counterclockwise direction.

The cross product of (0, 0, 1) and (x, y, 0) is (-y, x, 0), so the direction we want to go in the xy plane is (-y, x). We divide this vector by its length to get a vector of length 1, then add it to our previous point.

Here is a code written in Python to plot the spiral

import matplotlib.pyplot as plt
def vertex(x, y):
    h = (x**2 + y**2)**0.5
    return (x - y/h, y + x/h)
# base of the first triangle
plt.plot([0, 1], [0, 0])
N = 17
x_old, y_old = 1, 0

for n in range(1, N):
    x_new, y_new = vertex(x_old, y_old)
    # draw short side
    plt.plot([x_old, x_new], [y_old, y_new])
    # draw hypotenuse
    plt.plot([0, x_new], [0, y_new])
    x_old, y_old = x_new, y_new

How to generate and read QR Code in Python

QR code is a type of matrix barcode that is machine readable optical label which contains information about the item to which it is attached. In practice, QR codes often contain data for a locator, identifier, or tracker that points to a website or application, etc.

Problem Statement :

Generate and read QR codes in Python using qrcode and OpenCV libraries

Installing required dependencies:

pyqrcode module is a QR code generator. The module automates most of the building process for creating QR codes. This module attempts to follow the QR code standard as closely as possible. The terminology and the encoding used in pyqrcode come directly from the standard.

pip install pyqrcode

Install an additional module pypng to save image in png format:

pip install pypng

Import Libraries

import pyqrcode
import png
from pyqrcode import QRCode
import cv2
import numpy as np

Create QR Code:


# String which represents the QR code
s = ""

# output file name

filename = "qrcode.png"

# Generate QR Code

img = pyqrcode.create (s)

# Create and save the svg file naming "brqr.svg"

img.svg("brqr.svg", scale=8)

# Create and save the svg file naming "brqr.png"

img.png("brqr.png", scale=6)
qr code file named brqr.png

Read QR Code

Here we will be using OpenCV for that, as it is popular and easy to integrate with the webcam or any video.


# read the QRCODE image

img = cv2.imread('brqr.png')

# initialize the cv2 QRCode detector
detector = cv2.QRCodeDetector()

val, pts, st_code=detector.detectAndDecode(img)


Here the output

The output of code

Estimating Pi…

We will be exploring a very interesting, and simple for that matter, application of statistics to help us estimate the value of Pi.

For this method, we will be imagining a simple scenario. Imagine for a moment, that we have the unit circle, inside a square.

Unit circle inside a square

By having the unit circle, we immediately figure out that the area of the square will be four since the radius of the circle is defined at one, which means that our square will have sides with a value of two. Now here’s where things get interesting, if we take the ratio of both of areas, we end up getting the following:

\(\frac{{{A_c}}}{{{A_s}}} = \frac{{\pi {r^2}}}{{{{\left( {2r} \right)}^2}}} = \frac{\pi }{4}\)

Both of these geometric figures end up having a ratio of pi over four between them, which is an important value for our next step in which we use a bit of imagination.

from “Ordine e disordine” by Luciano De Crescenzo

For a moment, imagine that you have a circle inside a square on the ground; suppose it starts raining. Some drops will most likely fall inside the circle and others will likely fall inside the square but outside the circle. Using this concept is how we will code our estimator pi, throwing some random numbers into the unit circle equation, as shown below:


Furthermore, taking a ratio of throws that landed inside our circle and the total number of throws, we can then formulate the following:

\(\frac{{N{\text{of throws inside circle}}}}{{N {\text{of total throws}}}} = \frac{{{N_C}}}{{{N_T}}}\)

And by combining our ratio between the unit circle with the square with this new equation, we can assume the next equation

\(\frac{{{N_C}}}{{{N_T}}} \simeq \frac{\pi }{4} \Leftrightarrow \frac{{4{N_C}}}{{{N_T}}} \simeq \pi \)

With this equation, we can finally start coding our estimator and see how close we can get to pi’s actual value.

The code

import matplotlib.pyplot as plt
import numpy as np
import time as t
from import Bar

tin = t.time()

# numeri di punti della simulazione

# Vettore coordinate  x e y dei punti casuali
x = np.random.rand(n)
y = np.random.rand(n)

Pi_Greco = np.zeros(n)
#Vettore distanza

d= (x**2 + y**2)**(1/2)

Ps = 0
Pq = 0

bar = Bar('Processing', max=n)

for i in range(n):
	Pq = Pq + 1
	if d[i] < 1:
		Ps = Ps+1
	Pi_Greco [i] = 4*(Ps/Pq)


Pi_Greco_Reale = np.ones(n)*np.pi

tfin = t.time()

print('Valore di Pi Greco: ', 4*Ps/Pq)
print('elapsed time: ', round(tfin-tin, 3))

plt.plot(Pi_Greco, 'red', label='estimate of Pi')
plt.plot(Pi_Greco_Reale, 'green', label='Exact value')

plt.title('Monte Carlo simulation')

Results for our pi estimation. Pi is represented by the green line, while red represents our estimations.

It’s quite interesting to see our estimation start with low accuracy but as we increase our attempts, we start to get a convergence on the the value of pi, as shown by our green line. Statistical methods like this, and other more complex versions, are nice tools to understand and experiment within the world of physics. I highly recommend taking a look into some Statistical Mechanic concepts to see the beauty behind the application of statistics and probability in physics, and maybe take some time play with these concepts in Python!