# Exercises Laboratorio di Calcolo: Practicing NumPy + Matplotlib¶

By Hendrik Speleers, AY 2025-2026

## 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 [ ]: