4. Le module numpy I

Essayez pour chaque exercice de trouver plusieurs façons de répondre et comparez les temps d’exécution en utilisant la commande magique timeit. Cette commande magique se place en début d’une cellule et mesure le temps d’exécution de la cellule complète lorsqu’elle est précédée de %%. Si elle est précédée de %, elle n’agit que sur la commande qui suit.

Voici un exemple :

%%timeit
commande1
commande2
...

ou bien

%timeit commande

import numpy as np

4.1. Exercice 1

Question

Créez et affichez le tableau numpy permettant de stocker les valeurs du vecteur \(x\) avec \(x_i = 2^i\), \(1\leq i \leq 100\).

Show code cell source Hide code cell source

def sol1(n):
    """Solution 1 (compréhension de liste)"""
    return np.array([2**(i+1) for i in range(n)], dtype=float)

def sol2(n):
    """Solution 2 (boucle à la main)"""
    x = np.zeros(n)
    for i in range(n):
        x[i] = 2**(i+1)
    return x

def sol3(n):
    """Solution 3 (à l'aide d'une fonction)"""
    return np.fromfunction(lambda i: 2**(i+1), (n,))

def sol4(n):
    """Solution 4 (calcul vectoriel)"""
    return 2.**np.arange(1,n+1)

##############################################################################
############### désactiver pour ne pas alourdir la compilation ###############
# n = 1000
# for k, solk in enumerate([sol1, sol2, sol3, sol4]):
#     print(f"Solution {k+1} : ", end='')
#     %timeit x = solk(n)
##############################################################################

Réponse

Sur mon ordinateur, voici à titre informatif les temps relevés pour n=1000:

• Solution 1 : 636 µs ± 218 ns per loop

• Solution 2 : 664 µs ± 2.46 µs per loop

• Solution 3 : 10.5 µs ± 8.09 ns per loop

• Solution 4 : 1.66 µs ± 2.33 ns per loop

4.2. Exercice 2

Rappel : une matrice symétrique est une matrice telle que \(A_{i,j} = A_{j, i}\), pour tout \(i, j\).

Questions

1. Proposez une fonction mat_alea(N, a, b) qui prend en argument un entier N et deux réels a et b et qui retourne une matrice aléatoire de taille \(N\times N\) constituée de réels compris entre \(a\) et \(b\).

2. Proposez une fonction mat_alea_sym(N, a, b) qui prend en argument un entier N et deux réels a et b et qui retourne une matrice aléatoire symétrique de taille \(N\times N\) constituée de réels compris entre \(a\) et \(b\).

Show code cell source Hide code cell source

def mat_alea(N, a, b):
    """
    retourne une matrice aléatoire de taille NxN 
    dont les valeurs sont comprises entre a et b
    
    temps d'exécution pour N = 1000 : 
    5.44 ms ± 14.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    """
    return a + (b-a) * np.random.rand(N, N)

def mat_alea_loop(N, a, b):
    """
    retourne une matrice aléatoire de taille NxN 
    dont les valeurs sont comprises entre a et b
    
    temps d'exécution pour N = 1000 : 
    382 ms ± 2.94 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
    """
    A = np.empty((N, N))
    for i in range(N):
        for j in range(N):
            A[i, j] = a + (b-a) * np.random.rand()
    return A

a, b = -1, 1
N = 5
print(mat_alea(N, a, b))
# N = 1000
# %timeit mat_alea(N, a, b)
# %timeit mat_alea_loop(N, a, b)

[[-0.21120294 -0.50572425 -0.51560088  0.69948537  0.25184491]
 [ 0.64566727  0.92814932 -0.03742505 -0.86137614 -0.26908111]
 [ 0.54818252  0.87696546 -0.22646724 -0.87638037 -0.7954247 ]
 [ 0.71273054  0.73355715  0.91121208 -0.63565615  0.43731694]
 [ 0.32919775  0.89195629 -0.26698383  0.19107224 -0.74627066]]

Réponse

Sur mon ordinateur, voici à titre informatif les temps relevés pour N=1000:

• Solution 1 : 5.44 ms ± 14.9 µs per loop

• Solution 2 : 382 ms ± 2.94 ms per loop

Show code cell source Hide code cell source

def mat_alea_sym(N, a, b):
    """
    retourne une matrice aléatoire de taille NxN 
    dont les valeurs sont comprises entre a et b
    
    temps d'exécution pour N = 1000 : 
    7.18 ms ± 138 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    """
    A = a + (b-a) * np.random.rand(N, N)
    return .5*(A + A.T)

def mat_alea_sym_loop(N, a, b):
    """
    retourne une matrice aléatoire de taille NxN 
    dont les valeurs sont comprises entre a et b
    
    temps d'exécution pour N = 1000 : 
    268 ms ± 959 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
    """
    A = np.empty((N, N))
    for i in range(N):
        for j in range(i):
            A[i, j] = a + (b-a) * np.random.rand()
            A[j, i] = A[i, j]
    return A

a, b = -1, 1
N = 5
print(mat_alea_sym(N, a, b))
# N = 1000
# %timeit mat_alea_sym(N, a, b)
# %timeit mat_alea_sym_loop(N, a, b)

[[ 0.72925015  0.15903069 -0.43084435 -0.75469104 -0.05392134]
 [ 0.15903069  0.56643465  0.09314903  0.54806696 -0.06159308]
 [-0.43084435  0.09314903 -0.08034734  0.20382448 -0.1065672 ]
 [-0.75469104  0.54806696  0.20382448  0.74622941 -0.49541613]
 [-0.05392134 -0.06159308 -0.1065672  -0.49541613 -0.76771821]]

Réponse

Sur mon ordinateur, voici à titre informatif les temps relevés pour N=1000:

• Solution 1 : 7.18 ms ± 138 µs per loop

• Solution 2 : 268 ms ± 959 µs per loop

4.3. Exercice 3

Question

Créez la matrice suivante

\[\begin{split} \begin{pmatrix} 1 & 2 & 3 & \ldots & N-1 & N \\ 2 & 3 & 4 & \ldots & N & 1 \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ N-1 & N & 1 & \ldots & N-3 & N-2 \\ N & 1 & 2 & \ldots & N-2 & N-1 \end{pmatrix} \end{split}\]

Show code cell source Hide code cell source

def solution1(N):
    """Solution 1 (avec des boucles)"""
    A = np.empty((N, N))
    for i in range(N):
        for j in range(N):
            A[i, j] = 1 + ((i+j) % N)
    return A

def solution2(N):
    """Solution 2 (à l'aide d'une fonction)"""
    return np.fromfunction(lambda i, j: 1.+((i+j) % N), (N, N), dtype=int)

def solution3(N):
    """Solution 3 (en utilisant le caclul vectoriel)"""
    i, j = np.arange(N), np.arange(N)
    i.shape = (1, N)
    j.shape = (N, 1)
    return 1. + ((i+j) % N)

def solution4(N):
    """Solution 4 (calcul vectoriel et broadcast...)"""
    i = np.arange(N)
    return 1. + ((i[:, None] + i[None, :]) % N)

N = 10
print(solution3(N))

[[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
 [ 2.  3.  4.  5.  6.  7.  8.  9. 10.  1.]
 [ 3.  4.  5.  6.  7.  8.  9. 10.  1.  2.]
 [ 4.  5.  6.  7.  8.  9. 10.  1.  2.  3.]
 [ 5.  6.  7.  8.  9. 10.  1.  2.  3.  4.]
 [ 6.  7.  8.  9. 10.  1.  2.  3.  4.  5.]
 [ 7.  8.  9. 10.  1.  2.  3.  4.  5.  6.]
 [ 8.  9. 10.  1.  2.  3.  4.  5.  6.  7.]
 [ 9. 10.  1.  2.  3.  4.  5.  6.  7.  8.]
 [10.  1.  2.  3.  4.  5.  6.  7.  8.  9.]]

##############################################################################
############### désactiver pour ne pas alourdir la compilation ###############
# N = 1000
# for k, solk in enumerate([solution1, solution2, solution3, solution4]):
#     print(f"Solution {k+1} : ", end='')
#     %timeit A = solk(N)
##############################################################################

Réponse

Sur mon ordinateur, voici à titre informatif les temps relevés pour N=1000:

• Solution 1 : 146 ms ± 2.23 ms per loop

• Solution 2 : 4.18 ms ± 52 µs per loop

• Solution 3 : 3.81 ms ± 53.4 µs per loop

• Solution 4 : 3.79 ms ± 35.1 µs per loop

4.4. Exercice 4

Questions

1. Créez x un tableau 1D aléatoire composé de nombres réels entre \(-1\) et \(1\) de taille 6.

2. Testez les commandes suivantes et expliquez ce qu’elles font

i_p = x > 0
print(i_p)
print(x[i_p])

3. En utilisant la question 2, créez le tableau y tel que \(y_i = 1\) si \(x_i>0\) et \(y_i=-1\) si \(x_i<0\).

4. Pouvez-vous faire la même chose avec un tableau de dimension 2 ?

N = 6
x = -1 + 2 * np.random.rand(N)
print(x)
i_p = x > 0
print(i_p)
print(x[i_p])

[-0.46159711  0.43631082 -0.84922738  0.41184276  0.48133164  0.11474695]
[False  True False  True  True  True]
[0.43631082 0.41184276 0.48133164 0.11474695]

Show code cell source Hide code cell source

N = 6
for ndim in [1, 2]:
    shape = tuple([N]*ndim)
    print("*"*80)
    print(f"dimension {len(shape)}")
    x = -1 + 2 * np.random.rand(*shape)
    print(x)
    # Solution 1
    y = 1.*(x > 0) - 1.*(x < 0)
    print(y)
    # Solution 2
    y = np.zeros(x.shape)
    y[x > 0] = 1
    y[x < 0] = -1
    print(y)
    # Solution 3
    y = np.where(x > 0, 1., -1.)
    print(y)

********************************************************************************
dimension 1
[-0.58985031 -0.98908445  0.88726328 -0.20982467 -0.92480048 -0.7828979 ]
[-1. -1.  1. -1. -1. -1.]
[-1. -1.  1. -1. -1. -1.]
[-1. -1.  1. -1. -1. -1.]
********************************************************************************
dimension 2
[[-0.94968636 -0.20958602 -0.74424206 -0.28729228 -0.42325778 -0.54187424]
 [-0.60898661 -0.12812269 -0.71138921  0.06237934  0.80453012 -0.86208312]
 [-0.49987784 -0.49709624 -0.82497076  0.30560535 -0.80869656  0.61045513]
 [-0.64101306 -0.18911831 -0.19550509 -0.70315511 -0.16589988  0.08908123]
 [ 0.31617069 -0.66813723 -0.0959015  -0.33442626 -0.75599271  0.00893601]
 [ 0.67653458 -0.02964297  0.93796021  0.55892516 -0.95863191  0.12227249]]
[[-1. -1. -1. -1. -1. -1.]
 [-1. -1. -1.  1.  1. -1.]
 [-1. -1. -1.  1. -1.  1.]
 [-1. -1. -1. -1. -1.  1.]
 [ 1. -1. -1. -1. -1.  1.]
 [ 1. -1.  1.  1. -1.  1.]]
[[-1. -1. -1. -1. -1. -1.]
 [-1. -1. -1.  1.  1. -1.]
 [-1. -1. -1.  1. -1.  1.]
 [-1. -1. -1. -1. -1.  1.]
 [ 1. -1. -1. -1. -1.  1.]
 [ 1. -1.  1.  1. -1.  1.]]
[[-1. -1. -1. -1. -1. -1.]
 [-1. -1. -1.  1.  1. -1.]
 [-1. -1. -1.  1. -1.  1.]
 [-1. -1. -1. -1. -1.  1.]
 [ 1. -1. -1. -1. -1.  1.]
 [ 1. -1.  1.  1. -1.  1.]]

4.5. Exercice 5

Etant donné un vecteur \(x=(x_0, \ldots, x_{N+1})\), on appelle matrice de dérivation centrée la matrice \(A\) de taille \(N\times N\) définie par

\[\begin{split} A_{i, j} = \left\lbrace \begin{aligned} &\frac{1}{x_{i+1}-x_{i-1}}, && \text{si } 1\leq i=j-1 \leq N-1,\\ &-\frac{1}{x_{i+1}-x_{i-1}}, && \text{si } 2\leq i=j+1 \leq N,\\ &0 && \text{sinon}. \end{aligned} \right. \end{split}\]

Questions

1. Construisez un vecteur de points équi-répartis entre \(0\) et \(1\) (vous prendrez peu de points pour pouvoir afficher la matrice).

2. Perturbez les valeurs de ce vecteur par une variable aléatoire sans désordonner les points.

3. Construisez alors la matrice de dérivation centrée associée.

Pour vérifier votre construction, vous pourrez utilier ces valeurs :

x = [0.   0.19 0.4  0.58 0.82 1.  ]
A = [[ 0.          2.5         0.          0.        ]
     [-2.56410256  0.          2.56410256  0.        ]
     [ 0.         -2.38095238  0.          2.38095238]
     [ 0.          0.         -2.38095238  0.        ]]

Indication

Pour fabriquer un vecteur de points équi-répartis, pensez au choix aux commandes np.linspace, np.arange.

Show code cell source Hide code cell source

N = 5
x, dx = np.linspace(0, 1, N+1, retstep=True)
x[1:-1] += 0.25*dx*(-1+2*(10*np.random.rand(N-1) // 1)/10)
print(f"x = {x}")
xl, xm, xr = x[:-2], x[1:-1], x[2:]

# A = np.zeros((N-1, N-1))
# for i in range(N-2):
#     A[i, i+1] = 1/(xr[i]-xl[i])
# for i in range(1, N-1):
#     A[i, i-1] = -1/(xr[i]-xl[i])

A = np.diag(
    1/(xr[:-1]-xl[:-1]), k=1
) - np.diag(
    1/(xr[1:]-xl[1:]), k=-1
)
print(f"A = {A}")
print(A @ np.ones(A.shape[0]))

x = [0.   0.16 0.41 0.59 0.81 1.  ]
A = [[ 0.          2.43902439  0.          0.        ]
 [-2.3255814   0.          2.3255814   0.        ]
 [ 0.         -2.5         0.          2.5       ]
 [ 0.          0.         -2.43902439  0.        ]]
[ 2.43902439  0.          0.         -2.43902439]

4.6. Exercice 6

La matrice du Laplacien est une matrice essentielle en analyse numérique. En dimension 1 d’espace, elle est définie pour un vecteur \(x=(x_0, \ldots, x_{N+1})\) par

\[\begin{split} A_{i, j} = \left\lbrace \begin{aligned} & -\frac{2}{(x_{i+1}-x_i)(x_i-x_{i-1})}, && \text{si } 1\leq i=j \leq N,\\ & \frac{2}{(x_{i+1}-x_{i-1})(x_{i+1}-x_i)}, && \text{si } 1\leq i=j-1 \leq N-1,\\ & \frac{2}{(x_{i+1}-x_{i-1})(x_i-x_{i-1})}, && \text{si } 2\leq i=j+1 \leq N,\\ & 0, && \text{sinon}. \end{aligned} \right. \end{split}\]

Questions

1. Construisez un vecteur de points équi-répartis entre \(0\) et \(1\) (vous prendrez peu de points pour pouvoir afficher la matrice).

2. Perturbez les valeurs de ce vecteur par une variable aléatoire sans désordonner les points.

3. Construisez alors la matrice du Laplacien associée.

Pour vérifier votre construction, vous pourrez utiliser ces valeurs :

x = [0.   0.19 0.4  0.58 0.81 1.  ]
A = [[-50.12531328  23.80952381   0.           0.        ]
     [ 24.42002442 -52.91005291  28.49002849   0.        ]
     [  0.          27.100271   -48.30917874  21.20890774]
     [  0.           0.          20.70393375 -45.76659039]]

Show code cell source Hide code cell source

N = 5
x, dx = np.linspace(0, 1, N+1, retstep=True)
x[1:-1] += 0.25*dx*(-1+2*(10*np.random.rand(N-1) // 1)/10)
print(f"x = {x}")
xl, xm, xr = x[:-2], x[1:-1], x[2:]

# A = np.zeros((N-1, N-1))
# for i in range(N-1):
#     A[i, i] = -2/(xr[i]-xm[i])/(xm[i]-xl[i])
# for i in range(N-2):
#     A[i, i+1] = 2/(xr[i]-xl[i])/(xr[i]-xm[i])
# for i in range(1, N-1):
#     A[i, i-1] = 2/(xr[i]-xl[i])/(xm[i]-xl[i])

A = np.diag(
    -2/(xr-xm)/(xm-xl)
) + np.diag(
    2/(xr[:-1]-xl[:-1])/(xr[:-1]-xm[:-1]), k=1
) + np.diag(
    2/(xr[1:]-xl[1:])/(xm[1:]-xl[1:]), k=-1
)
print(f"A = {A}")
print(A @ np.ones(A.shape[0]))

x = [0.   0.24 0.37 0.63 0.79 1.  ]
A = [[-64.1025641   41.58004158   0.           0.        ]
 [ 39.44773176 -59.17159763  19.72386588   0.        ]
 [  0.          18.31501832 -48.07692308  29.76190476]
 [  0.           0.          33.78378378 -59.52380952]]
[-2.25225225e+01 -7.10542736e-15  0.00000000e+00 -2.57400257e+01]