LG

Soit le système linéaire $AX=b$ avec $A$ matrice carrée $N \times N$ et $b$ et $X$ vecteurs de $\mathbb R^N$

Peut se résoudre si $det A \neq 0$.

Soit on détermine $A^{-1}$, on obtient $X = A^{-1}b$.

$$A^{-1}=\frac 1 {\det A} \,^{t}\!{\rm {com} A}$$

Dans la pratique, cette méthode est trop longue et inappliquée.

$$A = LU$$

• $L$ matrice triangulaire inférieure,
• $U$ matrice triangulaire supérieure avec des 1 sur la diagonale.

On détermine $LY=b$ puis $UX=Y$.

Lorsque L est déterminée, on ne peut calculer $y_i$ que si les $y_i$ inférieurs sont calculés :

$$y_i = \frac{b_i - \sum_{j=1}^{i-1}l_{i,j}y_j}{l_{ii}}$$

Si $l_{i,i} = 0$, deux possibilités : le numérateur de $y_i$ est différent de 0 alors, il n'y a aucune solution. Si le numérateur vaut 0, alors, il y a une infinité de résultats et on peut mettre n'importe quelle valeur à $y_i$.

Complexité : $N^2/2$ multiplications et additions, $N$ divisions.

Puis on applique le même principe avec la matrice supérieure $U$ pour déterminer le vecteur $x$.

Implémentation en python pour la résolution d'un système $Ax=b$

import numpy as np
 
def fsub(A,b):
    if len(A) != len(A[0]):
        raise ValueError('La matrice n''est pas carrée.')
    if len(A) != len(b):
        raise ValueError('La matrice A et le vecteur b n''ont pas la même dimension.')
 
    N = len(A)
    x = np.empty(N)
    sol = 1
    for i in range(0, N):
        somme = b[i] - np.sum(A[i][0:i] * x[0:i])
        if somme != 0 and A[i][i] == 0:
            return 0, np.nan
        if somme == 0 and A[i][i] == 0:
            x[i] = 0
            sol = -1
        else:
            x[i] = somme / A[i][i]
 
    return sol, x
 
def bsub(A,b):
    if len(A) != len(A[0]):
        raise ValueError('La matrice n''est pas carrée.')
    if len(A) != len(b):
        raise ValueError('La matrice A et le vecteur b n''ont pas la même dimension.')
 
    N = len(A)
    x = np.empty(N)
    sol = 1
    for i in range(N - 1, -1, -1):
        if A[i][i] == 0:
            determinant_null = True
        somme = b[i] - np.sum(A[i][i+1:N] * x[i+1:N])
# Pas de solution
        if somme != 0 and A[i][i] == 0:
            return 0, np.nan
# Infinité de solutions
        if somme == 0 and A[i][i] == 0:
            x[i] = 0
            sol = -1
        else:
            x[i] = somme / A[i][i]
 
    return sol, x
 
fsub(np.array([[1, 0, 0, 0],[1, 1, 0, 0],[1, 1, 1, 0],[1, 1, 1, 1]]), np.array([1, 2, 3, 4]))
fsub(np.array([[1, 0, 0, 0],[1, 0, 0, 0],[1, 1, 1, 0],[1, 1, 1, 1]]), np.array([1, 2, 3, 4]))
fsub(np.array([[1, 0, 0, 0],[2, 0, 0, 0],[1, 1, 1, 0],[1, 1, 1, 1]]), np.array([1, 2, 3, 4]))

Implémentation en python pour inverser une matrice. On applique le même algorithme ci-dessus mais avec B matrice identité.

import numpy as np
import math
import cmath
 
def pivotgauss(A):
    if len(A) != len(A[0]):
        raise ValueError('La matrice n''est pas carrée.')
 
    N = len(A)
    B = np.identity(N, dtype=np.float64)
    for i in range(0, N):
 
        indice = np.argmax(np.abs(A[i:,i]))
 
        if (A[i+indice, i+indice] == 0):
            A.fill(np.NaN)
#           C'est inacceptable car on inverse la matrice (cf. tp2).
            return (0, A)
 
        A[[i,i+indice]] = A[[i+indice,i]]
        B[[i,i+indice]] = B[[i+indice,i]]
 
        for j in range(i+1, N):
            multiplicateur = A[j][i]/A[i][i]
            B[j,:] -= B[i,:]*multiplicateur
            A[j][i:] -= A[i][i:]*multiplicateur
 
        B[i,:] /= A[i,i]
        A[i,:] /= A[i,i]
 
    for i in range(N-1, -1, -1):
        for j in range(0, i):
            B[j,:] -= A[j, i]*B[i,:]
    return (1, B)

La méthode du pivot a deux intérêts : éviter une division par 0 (voir ci-dessus) et réduire les erreurs numériques en triant les pivots par ordre décroissant.

Ci-dessous l'algorithme applique la méthode du pivot en même temps que la triangularisation.

Pour chaque colonne, on cherche la plus grande valeur absolue dans le triangle qu'on veut annuler et on inclut la diagonale. On permute la ligne puis on applique des permutations de lignes pour annuler toutes les valeurs de la colonne et la diagonale.

import numpy as np
 
def pivotgauss(A, b):
    if len(A) != len(A[0]):
        raise ValueError('La matrice n''est pas carrée.')
    if len(A) != len(b):
        raise ValueError('La matrice A et le vecteur b n''ont pas la même dimension.')
 
    N = len(A)
    for i in range(0, N):
 
        indice = np.argmax(np.abs(A[i:,i]))
 
        if (A[i+indice, i+indice] == 0):
#           C'est acceptable car c'est bsub qui déterminera
#           Si c'est aucune solution ou une infinité (cf. tp3).
            continue
 
        A[[i,i+indice]] = A[[i+indice,i]]
        b[i], b[i+indice] = b[i+indice], b[i]
 
        for j in range(i+1, N):
            multiplicateur = A[j][i]/A[i][i]
            b[j] = b[j] - b[i]*multiplicateur
            A[j][i:] = A[j][i:] - A[i][i:]*multiplicateur

Triangularisation :

def gauss(A, b):
    pivotgauss(A, b)
    (x, sol) = bsub(A, b)
    return (x, sol)

Résolution de systèmes linéaires par des méthodes directes : Gauss, LU Archive 14/10/2018

Exemple :

Université Laval avec une matrice $4*4$ Archive 14/10/2018

A doit être inversible (pas de pivot nul).

Étape 1 : $$ \begin{align*} &j = 1 \\ &i = 1 \quad &l_{jj} &= a_{jj} \\ &\forall i \in [2,N] \quad &u_{ji} &= \frac{a_{ji}}{l_{jj}} \end{align*} $$ Étape 2 : $$ \begin{align} \forall j &\in [2,N-1] : \\ i &= 1 \quad & l_{ji} &= a_{ji} \\ \forall &i \in [2, j] \quad & l_{ji} &= a_{ji} - \sum_{k=1}^{i-1}l_{jk}u_{ki} \\ \forall &i \in [j+1, N] \quad & u_{ji} &= \frac{a_{ji} - \sum_{k=1}^{j-1}l_{jk}u_{ki}}{l_{jj}} \end{align} $$ Étape 3: $$ \begin{align} \forall i &\in [1,N] : \\ j &= N \\ l_{ji} &= a_{ji} - \sum_{k=1}^{i-1}l_{jk}u_{ki} \\ \end{align} $$ Complexité : $\frac{N^3}{3}$ multiplications, $N$ divisions.

[1] §2.2 L'algorithme de Gauss.

À propos des méthodes de décomposition de type GAUSS Archive 14/​10/​2018

Si l'un des $l_{jj}$ vaut 0, il faut appliquer un pivot.

Implémentation du code en python :

def factorlu(A):
    L = np.identity(len(A))
    U = A.copy()
 
    for j in range(len(A)):
 
        if (U[j][j] == 0):
            L.fill(np.nan)
            U.fill(np.nan)
            return (L,U,0)
 
        for i in range(j+1,len(A)):
            L[i,j]=U[i][j]/U[j][j]
            U[i,j:] -=L[i,j]*U[j,j:]
 
    return (L,U,1)

Valable pour une matrice carrée symétrique définie positive.

$$A = C D {}^t C$$

• $C$ matrice triangulaire inférieure ayant des 1 sur la diagonale.
• $D$ matrice diagonale

Soit $\forall i \in [1,N]$ $$ d_i = a_{ii} - \sum_{j=1}^{i-1} l_{ij}^2 d_j \\ \forall j \in [i+1,N]\\ l_{ji} = (a_{ji} - \sum_{k=1}^{i-1} l_{jk} l_{ik} d_k)/d_i $$

Exemple : Résolution de systèmes linéaires par des méthodes directes : LDL’ et Choleski en page 14 Archive 14/10/2018

Valable pour une matrice carrée symétrique définie positive.

$$A = H {}^t H$$

• $H$ matrice triangulaire inférieure.

Méthode :

The Cholesky algorithm Archive du 13/04/2020 le 28/04/2020

Soit $$A^i= \begin{pmatrix} I_{i-1} & 0 & 0 \\ 0 & a_{i,i} & b_i^T \\ 0 & b_i & B^i \end{pmatrix}$$

$$L_i:= \begin{pmatrix} I_{i-1} & 0 & 0 \\ 0 & \sqrt{a_{i,i}} & 0 \\ 0 & \frac{b_i}{\sqrt{a_{i,i}}} & I_{n-i} \end{pmatrix}$$

$$ A^{i+1}= \begin{pmatrix} I_{i-1} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & B^i - \frac{b_i b_i^T}{a_{i,i}} \end{pmatrix} $$

Exemple :

$$ A = A^1= \begin{pmatrix} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \end{pmatrix} $$

$$ L_1= \begin{pmatrix} \sqrt4 & 0 & 0 \\ \frac{12}{\sqrt4} & 1 & 0 \\ \frac{-16}{\sqrt4} & 0 & 1 \end{pmatrix}= \begin{pmatrix} 2 & 0 & 0 \\ 6 & 1 & 0 \\ -8 & 0 & 1 \end{pmatrix} $$

$$ A^2= \begin{pmatrix} 1 & 0 & 0 \\ 0 & B^1 - \frac{b_1 b_1^T}{a_{1,1}} & \cdots \\ 0 & \vdots & \ddots \end{pmatrix} \quad B^1= \begin{pmatrix} 37 & -43 \\ -43 & 98 \end{pmatrix} \quad b_1=\begin{pmatrix} 12 \\ -16 \end{pmatrix} \quad A^2= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 5 \\ 0 & 5 & 34 \end{pmatrix} $$

$$ L_2= \begin{pmatrix} 1 & 0 & 0 \\ 0 & \sqrt1 & 0 \\ 0 & \frac{5}{\sqrt1} & 1 \end{pmatrix}= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 5 & 1 \end{pmatrix} $$

$$ A^3= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & B^2 - b_2 b_2^T \end{pmatrix} \quad B^2= \begin{matrix} 34 \end{matrix} \quad b_2=\begin{matrix} 5 \end{matrix} \quad A^3= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 9 \end{pmatrix} $$

$$ L_3= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \sqrt9 \end{pmatrix}= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 3 \end{pmatrix}$$

$$ L = L_1 L_2 L_3 = \begin{pmatrix} 2 & 0 & 0 \\ 6 & 1 & 0 \\ -8 & 5 & 3 \end{pmatrix} $$

$$ \begin{align*} \forall i &\in [1,N] \\ j &= i \quad & L_{j,j} &= \sqrt{ A_{j,j} - \sum_{k=1}^{j-1} L_{j,k}^2 } \\ \forall j &\in [1,i-1] \quad & L_{i,j} &= \frac{1}{L_{j,j}} \left( A_{i,j} - \sum_{k=1}^{j-1} L_{i,k} L_{j,k} \right) \end{align*} $$ Exemple : $A = \begin{pmatrix} 4 & 12 & -16
12 & 37 & -43
-16 & -43 & 98
\end{pmatrix}$

$L = \begin{pmatrix} \sqrt 4 = 2 & 0 & 0 \\ \frac{12}{2} = 6 & \sqrt{37-6^2} = 1 & 0 \\ \frac{-16}{2} = -8 & \frac{-43-(-8)*6}{1} = 5 & \sqrt{98-8^2-5^2} = 3 \end{pmatrix}$

Implémentation en python :

def cholesky(A):
    n = len(A)
    L = np.zeros((n,n), dtype=np.complex64)
    L[0,0] = cmath.sqrt(A[0,0])
    L[1:,0] = A[1:,0]/L[0,0]
 
    for j in range(1,n):
        if (not np.isreal(A[j, j])):
# Matrice A non hermitienne
            L.fill(np.nan)
            return (L,-1)
 
        carre = A[j,j] - np.dot(L[j,:j],L[j,:j].transpose().conjugate())
        if (carre <= 0):
# Matrice non définie strictement positive
            L.fill(np.nan)
            return (L,0)
        L[j,j] = cmath.sqrt(carre)
 
        for i in range(j+1,n):
            if (A[i,j] != np.conjugate(A[j,i])):
# Matrice A non hermitienne
                L.fill(np.nan)
                return (L,0)
            L[i,j]=(A[i,j]-np.dot(L[i,:i],L[j,:i].transpose().conjugate()))/L[j,j]
 
    return (L,1)
 
 
A4 = np.array(
        [[  4,  12, -16],
         [ 12,  37, -43+2j],
         [-16, -43-2j,  98]], dtype=np.complex64)
(L, sol) = cholesky(A4)
print ("L",L)
print ("sol",sol)
print (L.dot(L.transpose().conjugate()))

$QR = A$ avec R matrice diagonale supérieure et Q matrice orthogonale.

Soit $A_0 = A$

$\alpha_i = A_{i-1}[1]$ avec $A[1]$ correspondant à la première colonne de $A$ et $i \in [1,N-1]$.

$v_i = \alpha_i \pm \|\alpha_i\| e_1 $ avec $e_1$ vecteur nul sauf la première ligne vaut 1. Le signe à prendre est le même que la première ligne de $A[1]$ pour éviter une valeur nulle.

$$Q_i = I - 2 \frac{v_i v_i^t }{\|v_i\|_2^2}$$

$$A_i = \prod_{k=i-1}^{1}Q_i A$$

$$A_i= \begin{pmatrix} U_{i-1} & B \\ 0 & A_{i+1} \\ \end{pmatrix}$$ avec $U_x$ : matrice triangulaire supérieure de dimension $x$, et $b$ vecteur.

$$R = A_{N-1}$$

$$Q = \prod_{k=1}^{N-1}Q_i$$

Décomposition QR Archive du 16/05/2019 le 28/04/2020

Factorisation QR et systèmes surdéterminés Archive du 06/02/2019

Exemple :

$$A=A_0= \begin{pmatrix} 12 & -51 & 4 \\ 6 & 167 & -68 \\ -4 & 24 & -41 \\ \end{pmatrix}$$

$$\alpha_1 = \begin{pmatrix} 12 \\ 6 \\ -4 \\ \end{pmatrix}$$

$$\|\alpha_1\| = 14$$

$$v_1 = \alpha_1 \pm \|\alpha_1\| e_1 = \begin{pmatrix} 26 \\ 6 \\ -4 \\ \end{pmatrix}$$

$$Q_1 = I - 2 \frac{v_1 v_1^T}{\|v_1\|_2^2} = \begin{pmatrix} -6/7 & -3/7 & 2/7 \\ -3/7 & 82/91 & 6/91 \\ 2/7 & 6/91 & 87/91 \\ \end{pmatrix}$$

$$A_1 = Q_1 A_0 = \begin{pmatrix} -14 & -21 & 14 \\ 0 & 2261/13 & -854/13 \\ 0 & 252/13 & -553/13 \\ \end{pmatrix}$$

$$\alpha_2 = \begin{pmatrix} 2261/13 \\ 252/13 \\ \end{pmatrix}$$

$$\|\alpha_2\| = 175$$

$$v_2 = \alpha_2 \pm \|\alpha_2\| e_1 = \begin{pmatrix} 4536/13 \\ 252/13 \\ \end{pmatrix}$$

$$Q_2 = I - 2 \frac{v_2 v_2^T}{\|v_2\|_2^2} = \begin{pmatrix} -323/325 & -36/325 \\ -36/325 & 323/325 \\ \end{pmatrix}$$

$$R = A_2 = Q_2 A_1 = \begin{pmatrix} -14 & -21 & 14 \\ 0 & -175 & 70 \\ 0 & 0 & -35 \\ \end{pmatrix}$$

$$Q = \begin{pmatrix} -6/7 & -3/7 & 2/7 \\ -3/7 & 82/91 & 6/91 \\ 2/7 & 6/91 & 87/91 \\ \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & -323/325 & -36/325 \\ 0 & -36/325 & 323/325 \\ \end{pmatrix} = \begin{pmatrix} -6/7 & 69/175 & 58/175 \\ -3/7 & -158/175 & -6/175 \\ 2/7 & -6/35 & 33/35 \\ \end{pmatrix}$$

$\forall i \in [1,N]$

$a_i$ les vecteurs colonnes de A.

$u_i = a_i - \sum_{j=1}^{i-1} (a_i \cdot e_j) e_j$

$e_i = \frac{u_i}{\|u_i\|}$

$Q = \left[e_1 | … | e_N \right]$

$ \begin{align} & \forall i > j & R_{i,j} = 0 \\ & \forall i \le j & R_{i,j} = a_j e_i \\ \end{align} $

QR Decomposition with Gram-Schmidt Archive du 16/02/2019

Exemple :

$$A= \begin{pmatrix} 12 & -51 & 4 \\ 6 & 167 & -68 \\ -4 & 24 & -41 \\ \end{pmatrix}$$

$$u_1 = \begin{pmatrix} 12 \\ 6 \\ -4 \\ \end{pmatrix}$$

$$\|u_1\| = 14$$

$$e_1 = \frac{u_1}{\|u_1\|} = \begin{pmatrix} 6/7 \\ 3/7 \\ -2/7 \\ \end{pmatrix}$$

$$a_2 \cdot e_1 = (-51)*6/7 + 167 * 3/7 + 24 * (-2) / 7 = 21$$

$$u_2 = a_2 - (a_2 \cdot e_1) e_1 = \begin{pmatrix} -51 - 21 * 6 / 7 \\ 167 - 21 * 3 / 7 \\ 24 - 21 * (-2) / 7 \\ \end{pmatrix} = \begin{pmatrix} -69 \\ 158 \\ 30 \\ \end{pmatrix} $$

$$\|u_2\| = 175$$

$$e_2 = \frac{u_2}{\|u_2\|}$$

$$a_3 \cdot e_1 = 4 * 6 / 7 + (-68) * 3 / 7 + (-41) * (-2) / 7 = -14$$

$$a_3 \cdot e_2 = 4 * (-69) / 175 + (-68) * 158 / 175 + (-41) * 30 / 175 = -70$$

$$u_3 = a_3 - (a_3 \cdot e_1) e_1 - (a_3 \cdot e_2) e_2 = \begin{pmatrix} 4 - (-14) * 6 / 7 - (-70) * (-69) / 175 \\ -68 - (-14) * 3 / 7 - (-70) * 158 / 175 \\ -41 - (-14) * (-2) / 7 - (-70) * 30 / 175 \\ \end{pmatrix} = \begin{pmatrix} -58 / 5 \\ 6/5 \\ -33 \\ \end{pmatrix} $$

$$\|u_3\| = 35$$

$$e_3 = \frac{u_3}{\|u_3\|}$$

$$Q = \begin{pmatrix} 6/7 & -69/175 & -58/175 \\ 3/7 & 158/175 & 6/175 \\ -2/7 & 6/ 35 & -33/ 35 \\ \end{pmatrix} $$

$$R = \begin{pmatrix} 14 & 21 & -14 \\ 0 & 175 & -70 \\ 0 & 0 & 35 \\ \end{pmatrix} $$

$\forall j \in [1,N-1]$ avec j désignant les colonnes.

$\forall i \in [j+1,N]$ par ordre décroissant avec désignant les lignes.

On cherche donc à annuler chaque terme sur le triangle inférieur.

Au départ $R_{1,N}$ vaut A.

$r_{i,j} = \sqrt{\left ( R^{ij}_{i,j} \right ) ^ 2 + \left ( R^{ij}_{i-1,j} \right ) ^ 2}$

On calcule la matrice de rotation $G^{i,j}$ : vaut la matrice identité sauf : $$G_{i,i} = \frac{R_{i-1,j}}{r_{i,j}}$$ $$G_{i-1,i-1} = \frac{R_{i-1,j}}{r_{i,j}}$$ $$G_{i,i-1} = -\frac{R_{i,j}}{r_{i,j}}$$ $$G_{i-1,i} = \frac{R_{i,j}}{r_{i,j}}$$

Le $R^{(n+1)}$ suivant vaut $G*R^{(n)}$.

An Example of QR Decomposition Archive du 17/02/2019

The QR Reduction Archive du 17/02/2019

Exemple :

$$R^{(0)} = A= \begin{pmatrix} 12 & -51 & 4 \\ 6 & 167 & -68 \\ -4 & 24 & -41 \\ \end{pmatrix}$$

$i=3, j=1$

$r_{3,1} = \sqrt{(-4)^2+6^2} = 2 \sqrt{13}$

$$G^{(1)} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{G^{(0)}_{i-1,j}}{r_{3,1}} & \frac{G^{(0)_{i,j}}}{r_{3,1}} \\ 0 & -\frac{G^{(0)}_{i,j}}{r_{3,1}} & \frac{G^{(0)}_{i-1,j}}{r_{3,1}} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{6}{2 \sqrt{13}} & \frac{-4}{2 \sqrt{13}} \\ 0 & -\frac{-4}{2 \sqrt{13}} & \frac{6}{2 \sqrt{13}} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{3}{\sqrt{13}} & \frac{-2}{\sqrt{13}} \\ 0 & \frac{2}{\sqrt{13}} & \frac{3}{\sqrt{13}} \\ \end{pmatrix} $$

$$R^{(1)} = G^{(1)} R^{(0)} = \begin{pmatrix} 12 & -51 & 4 \\ 2 \sqrt{13} & 453/\sqrt{13} & -122/\sqrt{13} \\ 0 & 406/\sqrt{13} & -259/\sqrt{13} \\ \end{pmatrix} $$

$i=2, j=1$

$r_{2,1} = \sqrt{(12)^2+(2 \sqrt{13})^2} = 14$

$$G^{(2)} = \begin{pmatrix} \frac{G^{(0)}_{i-1,j}}{r_{2,1}} & \frac{G^{(0)_{i,j}}}{r_{2,1}} & 0 \\ -\frac{G^{(0)}_{i,j}}{r_{2,1}} & \frac{G^{(0)}_{i-1,j}}{r_{2,1}} & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} = \begin{pmatrix} \frac{12}{14} & \frac{2 \sqrt{13}}{14} & 0 \\ -\frac{2 \sqrt{13}}{14} & \frac{12}{14} & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} = \begin{pmatrix} \frac{6}{7} & \frac{\sqrt{13}}{7} & 0 \\ -\frac{\sqrt{13}}{7} & \frac{6}{7} & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} $$

$$R^{(2)} = G^{(2)} R^{(1)} = \begin{pmatrix} 14 & 21 & -14 \\ 0 & 483/\sqrt{13} & -122/\sqrt{13} \\ 0 & 406/\sqrt{13} & -259/\sqrt{13} \\ \end{pmatrix} $$

$i=3, j=2$

$r_{3,2} = \sqrt{(483/\sqrt{13})^2+(406/\sqrt{13})^2} = 175$

$$G^{(3)} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{G^{(0)}_{i-1,j}}{r_{3,2}} & \frac{G^{(0)_{i,j}}}{r_{3,2}} \\ 0 & -\frac{G^{(0)}_{i,j}}{r_{3,2}} & \frac{G^{(0)}_{i-1,j}}{r_{3,2}} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{483/\sqrt{13}}{175} & \frac{406/\sqrt{13}}{175} \\ 0 & -\frac{406/\sqrt{13}}{175} & \frac{483/\sqrt{13}}{175} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \frac{69}{25\sqrt{13}} & \frac{58}{25\sqrt{13}} \\ 0 & -\frac{58}{25\sqrt{13}} & \frac{69}{25\sqrt{13}} \\ \end{pmatrix} $$

$$R^{(2)} = G^{(2)} R^{(1)} = \begin{pmatrix} 14 & 21 & -14 \\ 0 & 175 & -70 \\ 0 & 0 & -35 \\ \end{pmatrix} $$

[1] : DESTUYNDER, Philippe. Méthodes numériques pour l'ingénieur. Lavoisier, 2010, 248 pages.

[2] : Timothy, Vismor. Matrix Algorithms. January 30, 2015, 69 pages. Archive 11/11/2018