CAMA : ma13 Système matriciel – Exercices

Cours du 04/05

Programmation vectorielle

Le but des exercices est

• d'avoir un programme qui donne la bonne réponse
• qui soit le plus rapide possible (et pour cela on utilise massivement Numpy)

En règle général si vous avez des for imbriqués c'est mauvais signe.

Méthode du pivot de Gauss partiel

L'ennoncé est dans le cours

Solution

def solve_gauss_partial(A, b):   # on prend le max dans la colonne i parmi les lignes en dessous (plus facile)
    for i in range(len(A)-1):
        pivot = np.argmax(np.abs(A[i:, i]))  # il n'y a que 3 lignes à ajouter pour échanger les lignes
        A[[i, pivot]] = A[[pivot, i]]
        b[[i, pivot]] = b[[pivot, i]]
        E = np.diag(np.array([1.,] * len(A), dtype=A.dtype))
        coefs = - A[i+1:,i] / A[i,i]
        E[i+1:,i] = coefs
        A[i:, i:] = E[i:,i:] @ A[i:,i:]
        b[i+1:] += coefs * b[i]   # multiplication terme à terme
    # A est maintenant triangulaire surpérieur
    res = np.zeros(len(b), dtype=b.dtype)
    res[-1] = b[-1] / A[-1,-1]
    for i in range(len(A)-1)[::-1]:
        res[i] = (b[i] - A[i,i+1:] @ res[i+1:]) / A[i,i]
    return res

e = 1E-6
A = np.array([[e, 1], [1, 2]], dtype='float32')
b = np.array([1., 3.], dtype='float32')
print(f"A\n {A} \nb\n {b}\n")
x = solve_gauss_partial(A, b)
print('solution : ',x)
print('vérification\n', A@x)

A
 [[0.000001 1.      ]
 [1.       2.      ]] 
b
 [1. 3.]

solution :  [1.0000019  0.99999905]
vérification
 [3.       0.999997]

Factorisation de Choleski

Il s'agit de décomposer A en

$A=B{B}^T$ avec B une matrice triangulaire inférieure. Cela n'est possible
que si la matrice A est symétrique et définie positive (c'est d'ailleurs un facon de vérifier qu'une
matrice est définie positive).

Écrire l'algorithme de Choleski qui prend A et retourne B (pour deviner l'algorithme, essayez de trouver les
coefficients de B à partir des coefficients de A sur une matrice A 4x4).

Solution

$$A=B{B}^T=\left[\begin{array}{cccc}{b}_{11}& 0& \dots & 0\\ b_{21}& b_{22}& \dots & 0\\ & ⋮& \\ b_{n1}& b_{n2}& \dots & b_{n,n}\end{array}\right]\left[\begin{array}{cccc}{b}_{11}& b_{21}& \dots & b_{n1}\\ 0& b_{22}& \dots & b_{n2}\\ & ⋮& \\ b_{n1}& b_{n2}& \dots & b_{n,n}\end{array}\right]=\left[\begin{array}{cccc}{b}_{11}^2& b_{11}{b}_{21}& \dots & b_{11}{b}_{n1}\\ x& \sum _{i=1}^2{b}_{2i}^2& \dots & \sum _{i=1}^2{b}_{2i}{b}_{ni}\\ & & ⋮& \\ x& x& \dots & \sum _{i=1}^2{b}_{n,i}^2\end{array}\right]$$
avec

$x$ la même valeur que de l'autre coté de la diagonale

On voit que

$b_{11}=\sqrt{a_{11}}$ et maintenant qu'on a

$b_{11}$ on peut trouver toute la première ligne de

$B^T$ :

$b_{j1}=a_{1j}/b_{11}$.

Une fois qu'on connait la première ligne de

$B^T$ , on s'attaque à la deuxième en commencant par trouver

$b_{22}$ puis ensuite tous les autres éléments de la ligne comme on a fait pour la première ligne.

On a donc dans le cas général :

• $b_{ii}=\sqrt{a_{ii}-\sum _{k=1}^{i-1}{b}_{ik}^2}$
• $b_{ji}=a_{ij}-\sum _{k=1}^{i-1}{b}_{ik}{b}_{jk}/b_{ii}=a_{ij}-\sum _{k=1}^{i-1}{b}_{ik}{b}_{kj}^T/b_{ii}\mathrm{\forall }j>i$

def Choleski(A):
    B = np.zeros(A.shape)
    for i in range(len(A)):
        B[i,i] = np.sqrt(A[i,i] - np.sum(np.square(B[i, :i])))         # garanti ok car A est def positive
        B[i+1:, i] = (A[i, i+1:] - B[i, :i] @ B.T[:i, i+1:]) / B[i,i]  # les ∑ sous forme de prod. scalaire
    return B

Matrice symetrique

Rappel : pas de boucles for imbriquées mais des opérations vectorielles et matricielles (sur des sous-matrices).

Créer une matrice symétrique définie positive est vérifier que votre programme marche bien.

Solution

A = np.random.randint(10, size=(4,4))
A = A + A.T                             # symmétrique
A = A + np.diag(A.sum(axis=0))          # diagonale dominante
print('A:\n', A)
B = Choleski(A)
print('B\n', B)
print('vérification\n', B @ B.T)

A:
 [[55  8 18  5]
 [ 8 33  7 10]
 [18  7 54  9]
 [ 5 10  9 28]]
B
 [[7.4161984871 0.           0.           0.          ]
 [1.0787197799 5.6423721639 0.           0.          ]
 [2.4271195049 0.7765914857 6.8924593995 0.          ]
 [0.6741998625 1.6434093681 0.8831939788 4.9055711788]]
vérification
 [[55.  8. 18.  5.]
 [ 8. 33.  7. 10.]
 [18.  7. 54.  9.]
 [ 5. 10.  9. 28.]]

Améliorer Jacobi

Lorsqu'on écrit une itération de la méthode de Jacobi avec l'ensemble des coefficients, on constate que
si on calcule la nouvelle valeur de x élément par élement alors lorsqu'on veut mettre à jour x[1],
on connait déjà x[0]. Idem lorsqu'on met à jour x[2] on connait déjà x[0] et x[1], etc.

L'idée de la version amméliorée de Jacobi est d'utiliser la nouvelle valeur de x[0] et non pas l'ancienne
comme c'est le cas dans l'algorithme du cours. Ainsi en utilisant des valeurs mise à jour on peut espérer
converger plus vite.

Dans ce chaque itération demande un calcul ligne par ligne et donc une boucle for dans la boucle for sur
les itérations.

Test d'arrêt

On ajoutera un argument error à la fonction qui indique la précision désirée du résultat. Par
défaut sa valeur est de 1E-6 et pour offrir une bonne garantie on arrête l'algorithme lorsque

$||x_{t+1}-x_t||<\mathtt{\text{error}}/10$.

Solution

def Jacobi(A, b, error=1E-6, verbose=False):
    L = np.tril(A) 
    U = -np.triu(A, k=1) 
    if verbose:
        print(f"L:\n {L}\nU\n {U}\n")
    previous_x = np.zeros(len(b))
    x = np.random.random(len(b))
    err = (error / 10) ** 2
    while np.sum(np.square(x - previous_x))  > err:
        previous_x = x.copy()
        if verbose:
            print(f"x = {x}")
        # on résoud  L x = U x + b  avec L matrice triangulaire inférieure
        y = U @ x + b
        x[0] = y[0] / L[0,0]
        for i in range(1,len(L)):
            x[i] = (y[i] - L[i,:i] @ x[:i]) / L[i,i]
    return x

A = np.random.randint(10, size=(4,4))
A = A + np.diag(A.sum(axis=0))
b = A.sum(axis=1)                     # ainsi la solution est [1,1,1,1]
print('A:\n', A, "\nb:\n", b, "\n")

Jacobi(A,b, verbose=True)

A:
 [[24  2  1  7]
 [ 5 19  4  6]
 [ 9  2 20  9]
 [ 2  9  9 32]] 
b:
 [34 34 40 52] 

L:
 [[24  0  0  0]
 [ 5 19  0  0]
 [ 9  2 20  0]
 [ 2  9  9 32]]
U
 [[ 0 -2 -1 -7]
 [ 0  0 -4 -6]
 [ 0  0  0 -9]
 [ 0  0  0  0]]

x = [0.8870874823 0.8448958895 0.2146829205 0.8281640711]
x = [1.0957657001 1.194392389  1.0147923641 0.9351814319]
x = [1.0020897014 1.0168049182 1.0265474982 0.9876765266]
x = [1.0010877908 0.9980164155 1.0052544156 0.9990120918]
x = [1.0002345046 0.9991440665 1.000424625  1.000106649 ]
x = [1.0000225291 0.9998709979 0.9999547701 1.0000475947]
x = [0.999998753  0.9999948204 0.9999796615 1.0000072549]
x = [0.9999991631 1.000002211  0.9999968908 1.0000003049]
x = [0.9999998564 1.0000005961 0.9999998678 0.9999998785]
x = [0.9999999913 1.0000000685 1.0000000518 0.9999999667]

array([1.0000000018, 0.9999999991, 1.0000000142, 0.9999999961])