Exercises Laboratorio di Calcolo: Practicing NumPy + Matplotlib¶

By Hendrik Speleers, AY 2024-2025

Import all modules used in the exercises¶

In [1]:
import numpy as np
from matplotlib import pyplot as plt
In [2]:
%matplotlib inline

Exercise 1¶

Represent the following matrices in NumPy. Write for each of the matrices a Python function that generates a NumPy array in the right format. Avoid the use of explicit for-loops.

A matrix of size $(m, n)$ with checkerboard pattern¶

In [3]:
def matrix_checkerboard(m, n):
    matrix = np.zeros((m, n), dtype=int)
    matrix[::2, ::2] = 1
    matrix[1::2, 1::2] = 1
    return matrix
In [4]:
matrix_checkerboard(5, 6)
Out[4]:
array([[1, 0, 1, 0, 1, 0],
       [0, 1, 0, 1, 0, 1],
       [1, 0, 1, 0, 1, 0],
       [0, 1, 0, 1, 0, 1],
       [1, 0, 1, 0, 1, 0]])

A matrix of size $(m, n)$ containing the integer numbers $1, 2, ..., m·n$¶

In [5]:
def matrix_enumeration(m, n):
    matrix = np.arange(1, m*n+1, dtype=int).reshape((m, n))
    return matrix
In [6]:
matrix_enumeration(5, 5)
Out[6]:
array([[ 1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10],
       [11, 12, 13, 14, 15],
       [16, 17, 18, 19, 20],
       [21, 22, 23, 24, 25]])

A square matrix of size $(n, n)$ with the numbers $1, 2, ..., n$ on both the diagonal and the anti-diagonal¶

In [7]:
def matrix_cross(n):
    vals = np.arange(1, n+1, dtype=int)
    matrix = np.diag(vals)
    matrix[vals-1, n-vals] = vals
    return matrix
In [8]:
matrix_cross(5)
Out[8]:
array([[1, 0, 0, 0, 1],
       [0, 2, 0, 2, 0],
       [0, 0, 3, 0, 0],
       [0, 4, 0, 4, 0],
       [5, 0, 0, 0, 5]])

Exercise 2¶

Write Python functions with the described utilities. Avoid the use of explicit for-loops.

A function that returns the closest value (to a given scalar) in a 1D array¶

In [9]:
def find_closest(vals, val):
    index = np.abs(vals - val).argmin()
    return vals[index]
In [10]:
vals = np.array([1, 2, 5, 9, 7])
val = 4.7
find_closest(vals, val)
Out[10]:
5

A function that inserts 3 consecutive zeros between each element in a 1D array¶

In [11]:
def insert_zeros(vals, nz=3):
    len_zeros = (len(vals)-1) * (nz+1) + 1
    vals_zeros = np.zeros(len_zeros, dtype=vals.dtype)
    vals_zeros[::nz+1] = vals
    return vals_zeros
In [12]:
vals = np.array([1, 2, 5, 9, 7])
insert_zeros(vals, 3)
Out[12]:
array([1, 0, 0, 0, 2, 0, 0, 0, 5, 0, 0, 0, 9, 0, 0, 0, 7])

A function that swaps two rows of a 2D array¶

In [13]:
def swap_rows(matrix, row1, row2):
    matrix[[row1, row2]] = matrix[[row2, row1]]
    return matrix
In [14]:
matrix = np.array([[1, 2, 5], [4, 9, 7], [6, 3, 8]])
swap_rows(matrix, 0, 2)
Out[14]:
array([[6, 3, 8],
       [4, 9, 7],
       [1, 2, 5]])

Exercise 3¶

Write a Python program that visualizes a spiral (in red color) with bullets (in blue color) on top of it.

In [15]:
def plot_spiral(theta):
    t_line = np.linspace(0, theta, int(10*theta))
    x_line = t_line * np.cos(t_line)
    y_line = t_line * np.sin(t_line)
    plt.plot(x_line, y_line, 'r')
    t_mark = np.linspace(0, theta, int(theta))
    x_mark = t_mark * np.cos(t_mark)
    y_mark = t_mark * np.sin(t_mark)
    plt.plot(x_mark, y_mark, 'bo')
    plt.axis('equal')
    plt.axis('off')
In [16]:
theta = 20*np.pi
plot_spiral(theta)

Exercise 4¶

Write a Python program to visualize the first $n$ Legendre polynomials in the interval $[-1,1]$.

In [17]:
def plot_legendre(n):
    x = np.linspace(-1, 1, 100)
    for i in range(n):
        Li = np.polynomial.Legendre.basis(i)
        plt.plot(x, Li(x))
    plt.legend(['$L_{%d}$' % i for i in range(n)])
    plt.grid()
In [18]:
n = 6
plot_legendre(n)
In [19]:
def plot_legendre_alternative(n):
    x = np.linspace(-1, 1, 100)
    for i in range(n):
        Li = np.polynomial.Legendre.basis(i)
        plt.plot(x, Li(x), label='$L_{%d}$' % i)
    plt.legend()
    plt.grid()
In [20]:
n = 6
plot_legendre_alternative(n)
In [21]:
def plot_legendre_grid(n):
    m1 = int(np.sqrt(n))
    m2 = (n-1) // m1 + 1
    x = np.linspace(-1, 1, 100)
    for i in range(n):
        Li = np.polynomial.Legendre.basis(i)
        plt.subplot(m1, m2, i+1)
        plt.plot(x, Li(x))
        plt.axis([-1.2, 1.2, -1.2, 1.2])
        plt.text(0, 0.6, '$L_{%d}$' % i)
In [22]:
n = 6
plt.figure(figsize=(10, 6))
plot_legendre_grid(n)

Exercise 5: the Mandelbrot set¶

Write a Python function that returns a matrix containing for each complex point the minimal iteration step $i$ where the corresponding value $|z_{i+1}| > 2$ in the Mandelbrot recurrence relation.

In [23]:
def mandelbrot_iter(nx, ny, max_it=20, x1=-2.0, x2=1.0, y1=-1.5, y2=1.5):
    x = np.linspace(x1, x2, nx)
    y = np.linspace(y1, y2, ny)
    C = x[:,np.newaxis] + 1.0j*y[np.newaxis,:]
    Z = np.zeros((nx, ny), dtype=complex)
    M = np.ones((nx, ny), dtype=bool)
    M_it = np.full((nx, ny), max_it, dtype=int)
    for i in range(max_it):
        Z[M] = Z[M]**2 + C[M]
        Z_div = np.abs(Z) > 2
        M_it[M & Z_div] = i
        M[Z_div] = False
    return M_it

Then, visualize this matrix using plt.imshow and choose a high resolution (dpi) in order to fully appreciate the fractal structure.

In [24]:
M_it = mandelbrot_iter(2501, 2501, 200)
plt.figure(dpi=1200)
plt.imshow(M_it.T, cmap='flag')
plt.axis('off')
Out[24]:
(-0.5, 2500.5, 2500.5, -0.5)

Enhance your function so that you can also specify the considered rectangle in the complex plane. This will allow for making nice zooms of the mandelbrot fractal.

In [25]:
M_it = mandelbrot_iter(2501, 2501, 1200, -0.7450, -0.7425, 0.1300, 0.1325)
plt.figure(dpi=1200)
plt.imshow(M_it.T, cmap='flag')
plt.axis('off')
Out[25]:
(-0.5, 2500.5, 2500.5, -0.5)
In [ ]: