5. Le module numpy II

Dans ce TP, nous allons continuer à apprendre à manipuler les tableaux numpy.

import numpy as np

5.1. Exercice 1

Dans cet exercice, il est demandé de tester les fonctions linspace et arange du module numpy. Ne passez pas par les listes… Et faites attention au type des données stockées dans les tableaux (dtype) !

Question

1. Définir de la façon la plus rapide possible les vecteurs

\[\begin{split} \begin{aligned} v_1 &= \begin{pmatrix} 1&2&3&\ldots&16 \end{pmatrix}, \\ v_2 &= \begin{pmatrix} 0&0.2&0.4&\ldots&2 \end{pmatrix}, \\ v_3 &= \begin{pmatrix} 1&2&4&8&\ldots&64 \end{pmatrix}. \end{aligned} \end{split}\]

Show code cell source Hide code cell source

# solution avec linspace
v1 = np.linspace(1, 16, 16)
print(v1)
v1 = np.linspace(1, 16, 16, dtype=int)
print(v1)

# solution avec arange
v1 = np.arange(1, 17, 1)
print(v1)
v1 = np.arange(1, 17, 1, dtype=float)
print(v1)
v1 = np.arange(1, 17, 1.)
print(v1)

[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16.]
[ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16]
[ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16]
[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16.]
[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16.]

Show code cell source Hide code cell source

# solution avec linspace
v2 = np.linspace(0, 2, 11)
print(v2)
v2 = np.linspace(0, 2, 11, dtype=int)  # ATTENTION ICI
print(v2)

# solution avec arange
v2 = np.arange(0, 2.1, 0.2)
print(v2)
v2 = np.arange(0, 2.1, 0.2, dtype=int)  # ATTENTION ICI
print(v2)

[0.  0.2 0.4 0.6 0.8 1.  1.2 1.4 1.6 1.8 2. ]
[0 0 0 0 0 1 1 1 1 1 2]
[0.  0.2 0.4 0.6 0.8 1.  1.2 1.4 1.6 1.8 2. ]
[0 0 0 0 0 0 0 0 0 0 0]

Show code cell source Hide code cell source

# solution avec linspace
v3 = 2**np.linspace(0, 6, 7)
print(v3)
v3 = 2**np.linspace(0, 6, 7, dtype=int)
print(v3)

# solution avec arange
v3 = 2**np.arange(7)
print(v3)
v3 = 2**np.arange(7, dtype=float)
print(v3)
v3 = 2.**np.arange(7)
print(v3)

# solution avec logspace
v3 = np.logspace(0, 6, 7, base=2)
print(v3)

[ 1.  2.  4.  8. 16. 32. 64.]
[ 1  2  4  8 16 32 64]
[ 1  2  4  8 16 32 64]
[ 1.  2.  4.  8. 16. 32. 64.]
[ 1.  2.  4.  8. 16. 32. 64.]
[ 1.  2.  4.  8. 16. 32. 64.]

Question

2. Définir de la façon la plus simple possible les tableaux :

\[\begin{split} \begin{aligned} X &= \begin{pmatrix} 0.0 & 0.2 & 0.4 & \cdots & 1. \\ 0.0 & 0.2 & 0.4 & \cdots & 1. \\ 0.0 & 0.2 & 0.4 & \cdots & 1. \end{pmatrix}, \\ Y &= \begin{pmatrix} 0.0 & 0.0 & 0.0 & \cdots & 0.0 \\ 0.5 & 0.5 & 0.5 & \cdots & 0.5 \\ 1.0 & 1.0 & 1.0 & \cdots & 1.0 \end{pmatrix}, \\ A &= \begin{pmatrix} 1 & 2 & 4 & \cdots & 2^8 \\ 1 & 3 & 9 & \cdots & 3^8 \\ 1 & 5 & 25 & \cdots & 5^8 \end{pmatrix}. \end{aligned} \end{split}\]

Show code cell source Hide code cell source

X = np.tile(np.arange(0, 1.1, 0.2), (3, 1))
print(X)

[[0.  0.2 0.4 0.6 0.8 1. ]
 [0.  0.2 0.4 0.6 0.8 1. ]
 [0.  0.2 0.4 0.6 0.8 1. ]]

Show code cell source Hide code cell source

Y = np.tile(np.linspace(0, 1, 3).reshape((3, 1)), (1, 10))
print(Y)

[[0.  0.  0.  0.  0.  0.  0.  0.  0.  0. ]
 [0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5]
 [1.  1.  1.  1.  1.  1.  1.  1.  1.  1. ]]

Show code cell source Hide code cell source

ligne = np.array([2, 3, 5]).reshape(3, 1)
colonne = np.arange(9).reshape(1, 9)

# première solution avec tile
l = np.tile(ligne, (1, 9))
c = np.tile(colonne, (3, 1))
A = l ** c
print(A)

# deuxième solution avec repeat
l = np.repeat(ligne, 9, axis=1)
c = np.repeat(colonne, 3, axis=0)
A = l ** c
print(A)

# troisième solution en utilisant le broadcast
A = ligne ** colonne
print(A)

# quatrième solution avec une fonction
A = np.fromfunction(lambda i, j: (2+i+i//2)**j, shape=(3, 9), dtype=int)
print(A)

[[     1      2      4      8     16     32     64    128    256]
 [     1      3      9     27     81    243    729   2187   6561]
 [     1      5     25    125    625   3125  15625  78125 390625]]
[[     1      2      4      8     16     32     64    128    256]
 [     1      3      9     27     81    243    729   2187   6561]
 [     1      5     25    125    625   3125  15625  78125 390625]]
[[     1      2      4      8     16     32     64    128    256]
 [     1      3      9     27     81    243    729   2187   6561]
 [     1      5     25    125    625   3125  15625  78125 390625]]
[[     1      2      4      8     16     32     64    128    256]
 [     1      3      9     27     81    243    729   2187   6561]
 [     1      5     25    125    625   3125  15625  78125 390625]]

Question

3. Créer un tableau à 2 dimensions et \(10 \times 10\) éléments dont tous les éléments sont égaux à 0 sauf ceux aux bords et sur la diagonale qui seront pris égaux à 1.

Show code cell source Hide code cell source

T = np.eye(10)
T[[0, -1], :] = 1
T[:, [0, -1]] = 1
print(T)

[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 0. 0. 0. 0. 0. 0. 0. 1.]
 [1. 0. 1. 0. 0. 0. 0. 0. 0. 1.]
 [1. 0. 0. 1. 0. 0. 0. 0. 0. 1.]
 [1. 0. 0. 0. 1. 0. 0. 0. 0. 1.]
 [1. 0. 0. 0. 0. 1. 0. 0. 0. 1.]
 [1. 0. 0. 0. 0. 0. 1. 0. 0. 1.]
 [1. 0. 0. 0. 0. 0. 0. 1. 0. 1.]
 [1. 0. 0. 0. 0. 0. 0. 0. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]

Question

4. Créer un tableau à 2 dimensions et \(10 \times 10\) éléments dont tous les éléments sont égaux à 0 sauf les éléments diagonaux qui seront égaux à 2 et les éléments immédiatement au-dessus ou en-dessous de la diagonale qui seront égaux à \(-1\):

\[\begin{split} L = \begin{pmatrix} 2 & -1 & 0 & 0 & \dots & 0 & 0 \\ -1 & 2 & -1 & 0 & \ddots & & \vdots\\ 0 & \ddots & \ddots & \ddots & \ddots & & \\ \vdots & \dots & & \ddots & \ddots & \ddots & -1 \\ \\ 0 & \dots & & & 0 & -1 & 2 \end{pmatrix}. \end{split}\]

Indication

Pensez à utiliser la commande np.eye

Show code cell source Hide code cell source

n = 10
L = 2*np.eye(n, n, k=0) - np.eye(n, n, k=-1) - np.eye(n, n, k=1)
print(L)

[[ 2. -1.  0.  0.  0.  0.  0.  0.  0.  0.]
 [-1.  2. -1.  0.  0.  0.  0.  0.  0.  0.]
 [ 0. -1.  2. -1.  0.  0.  0.  0.  0.  0.]
 [ 0.  0. -1.  2. -1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0. -1.  2. -1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0. -1.  2. -1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0. -1.  2. -1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0. -1.  2. -1.  0.]
 [ 0.  0.  0.  0.  0.  0.  0. -1.  2. -1.]
 [ 0.  0.  0.  0.  0.  0.  0.  0. -1.  2.]]

Question

5. Créer un tableau T à 2 dimensions tel que \(T[i,j] = 1\) si \(i+j\) est paire, et \(T[i,j]=0\) sinon.

Show code cell source Hide code cell source

# première solution
T = np.fromfunction(lambda i, j: 1. - (i+j) % 2, shape=(10, 10))
print(T)

# deuxième solution
n = 10
l = np.arange(n).reshape((1, n))
c = np.arange(n).reshape((n, 1))
T = 1. - (l+c) % 2
print(T)

[[1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
 [1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
 [1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
 [1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
 [1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]]
[[1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
 [1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
 [1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
 [1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]
 [1. 0. 1. 0. 1. 0. 1. 0. 1. 0.]
 [0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]]

5.2. Exercice 2

Question

1. Déclarez des tableaux numpy permettant de stocker les matrices et vecteurs suivants

\[\begin{split} A = \begin{pmatrix} 7&0 \\ -1&5 \\ -1&2 \end{pmatrix}, \quad B = \begin{pmatrix} 1&4 \\ -4&0 \end{pmatrix}, \quad C = \begin{pmatrix} 7&3 \end{pmatrix}, \quad D = \begin{pmatrix} 8&2 \end{pmatrix}. \end{split}\]

Show code cell source Hide code cell source

A = np.array([[7, 0], [-1, 5], [-1, 2]])
B = np.array([[1, 4], [-4, 0]])
C = np.array([7, 3])
D = np.array([8, 2])

Question

2. Calculez avec python les produits \(AB\), \(AC^T\), \(CD^T\), \(DC^T\), \(DBC^T\), \(A^TA\) et \(AA^T\).

Indication

Vous pourrez utiliser le symbole @ ou la fonction dot pour les produits et la méthode transpose() pour calculer la transposée.

Show code cell source Hide code cell source

print("AB = ")
print(A @ B)
print("AC.T = ")
print(A @ C)
print("CD.T et DC.T valent")
print(C @ D, D @ C)
print("DBC.T = ")
print(D @ B @ C)
print("A.T A = ")
print(A.T @ A)
print("A A.T = ")
print(A @ A.T)

AB = 
[[  7  28]
 [-21  -4]
 [ -9  -4]]
AC.T = 
[49  8 -1]
CD.T et DC.T valent
62 62
DBC.T = 
96
A.T A = 
[[51 -7]
 [-7 29]]
A A.T = 
[[49 -7 -7]
 [-7 26 11]
 [-7 11  5]]

5.3. Exercice 3

Soit \(A\) une matrice carrée. Nous définissons l’exponentielle de la matrice \(A\), notée \(e^A\) la matrice définie par la somme infinie

\[ e^A = \sum_{n=0}^\infty \frac{1}{n!} A^N.\]

Question

• Proposez une fonction exp_matrice qui prend en argument une matrice A, qui vérifie que cette matrice est carrée et qui retourne une valeur approchée de l’exponentielle de A en tronquant la série.

• Vérifiez que votre fonction est correcte en calculant les exponentielles de ces trois matrices :

\[\begin{split} A_1 = \begin{pmatrix} 1&0 \\ 0&2 \end{pmatrix}, \quad A_2 = \begin{pmatrix} 1&1 \\ 0&1 \end{pmatrix}, \quad A_3 = \begin{pmatrix} 1&-1 \\ 1&2 \end{pmatrix}. \end{split}\]

Indication

Vous calculerez tous les termes \(A^n/n!\) jusqu’à ce que ce terme soit plus petit qu’une précision donnée (vous pourrez ajouter cette précision en argument optionnel à votre fonction). La matrice \(A^n/n!\) sera dite petite si np.linalg.norm(A^n/n!) < epsilon.

Show code cell source Hide code cell source

def exp_matrice(A, epsilon=1.e-7):
    """calcule l'exponentielle de la matrice A"""
    if A.shape[0] != A.shape[1]:
        print("Problem in exp_matrice: the matrix is not square")
        return None
    eA = np.eye(A.shape[0])
    An = A.copy()
    n = 1
    while np.linalg.norm(An) > epsilon:
        eA += An
        n += 1
        An[:] = An @ (A/n)
    return eA

A = np.array([[1, 0], [0, 2]], dtype=float)
print(exp_matrice(A))

[[2.71828183 0.        ]
 [0.         7.38905607]]

A = np.array([[1, 1], [0, 1]], dtype=float)
print(exp_matrice(A))

[[2.71828183 2.7182818 ]
 [0.         2.71828183]]

A = np.array([[1, -1], [1, 2]], dtype=float)
print(exp_matrice(A))

[[ 0.93244685 -3.94211457]
 [ 3.94211457  4.87456142]]

5.4. Exercice 4

Etant donnée deux matrices \(A\) de taille \((n_A,m_A)\) et \(B\) de taille \((n_B,m_B)\), le produit de Kronecker (ou produit tensoriel) des matrices \(A\) et \(B\) est la matrice de taille \((n_An_B,m_Am_B)\) définie par

\[ (A\otimes B)_{p,q} = A_{\alpha,\beta}B_{\gamma,\delta}, \]

où

\[\begin{split} \begin{aligned} p &=\alpha n_B + \gamma, && 0\leq\gamma<n_B, && 0\leq\alpha<n_A, && 0\leq p<n_An_B,\\ q &=\beta m_B + \delta, && 0\leq\delta<m_B, && 0\leq\beta<m_A && 0\leq q<m_Am_B. \end{aligned} \end{split}\]

Questions

1. Proposez une fonction qui calcule le produit de Kronecker de deux matrices entrées en argument.

2. Testez votre fonction en calculant le produit \(A\otimes B\) avec

\[\begin{split} A = \begin{pmatrix} 1&2\\0&-1 \\ 1&0\end{pmatrix} \qquad\text{et}\qquad B = \begin{pmatrix} 1&2&3\\2&0&-1 \end{pmatrix}.\end{split}\]

Show code cell source Hide code cell source

def kronecker(A, B):
    nA, mA = A.shape
    nB, mB = B.shape
    AB = np.zeros((nA*nB, mA*mB))
    ###########################################################################
    # solution 1 (tranquille avec des boucles)
    #--------------------------------------------------------------------------
    #for alpha in range(nA):
    #    for beta in range(mA):
    #        for gamma in range(nB):
    #            for delta in range(mB):
    #                p = alpha*nB + gamma
    #                q = beta*mB + delta
    #                AB[p, q] = A[alpha, beta]*B[gamma, delta]
    ###########################################################################
    # solution 2 (construction par blocs)
    #--------------------------------------------------------------------------
    #for alpha in range(nA):
    #    for beta in range(mA):
    #        AB[alpha*nB:(alpha+1)*nB, beta*mB:(beta+1)*mB] = A[alpha, beta]*B
    ###########################################################################
    # solution 3 (un peu plus sophistiquée)
    #--------------------------------------------------------------------------
    #AB = A.copy()
    #for k in range(B.ndim):
    #    AB = np.repeat(AB, B.shape[k], axis=k)
    #AB *= np.tile(B, A.shape)
    ###########################################################################
    # solution 4 (en utilisant le broadcast)
    #--------------------------------------------------------------------------
    AB = (A[:, None, :, None] * B[None, :, None, :]).reshape((nA*nB, mA*mB))
    ###########################################################################
    return AB

A = np.array([[1, 2], [0, -1], [1, 0]])
B = np.array([[1, 2, 3], [2, 0, -1]])
print(np.kron(A, B))
print(kronecker(A, B))

[[ 1  2  3  2  4  6]
 [ 2  0 -1  4  0 -2]
 [ 0  0  0 -1 -2 -3]
 [ 0  0  0 -2  0  1]
 [ 1  2  3  0  0  0]
 [ 2  0 -1  0  0  0]]
[[ 1  2  3  2  4  6]
 [ 2  0 -1  4  0 -2]
 [ 0  0  0 -1 -2 -3]
 [ 0  0  0 -2  0  1]
 [ 1  2  3  0  0  0]
 [ 2  0 -1  0  0  0]]