LG

Outils pour utilisateurs

Outils du site

Piste : • lieur
math:matrices:systeme_lineaire:directe

Différences

Ci-dessous, les différences entre deux révisions de la page.

Lien vers cette vue comparative

Les deux révisions précédentesRévision précédente
Prochaine révision		Révision précédente
math:matrices:systeme_lineaire:directe [2019/02/16 06:51] – [QR Householder (A non symétrique)] : Calcul du Q rootmath:matrices:systeme_lineaire:directe [2020/04/28 23:00] (Version actuelle) – mhtml -> html root
Ligne 1: Ligne 1:
- 
 =====Systèmes linéaires===== =====Systèmes linéaires=====
 ====Résolution==== ====Résolution====
Ligne 29: Ligne 28:
  
 Puis on applique le même principe avec la matrice supérieure $U$ pour déterminer le vecteur $x$. 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$
 +<code python>
 +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]))
 +</code>
 +
 +Implémentation en python pour inverser une matrice. On applique le même algorithme ci-dessus mais avec B matrice identité.
 +<code python>
 +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)
 +</code>
 +
 +
 +===Pivot===
 +
 +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.
 +
 +<code python>
 +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
 +</code>
 +
 +Triangularisation :
 +<code python>
 +def gauss(A, b):
 +    pivotgauss(A, b)
 +    (x, sol) = bsub(A, b)
 +    return (x, sol)
 +</code>
 +
 +[[https://asi.insa-rouen.fr/enseignement/siteUV/ananum/04syslindirect.pdf|Résolution de systèmes linéaires par des méthodes directes : Gauss, LU]] {{ math:matrices:systeme_lineaire:directe:04syslindirect.pdf |Archive 14/10/2018}}
 +
 +Exemple :
 +
 +[[https://www2.mat.ulaval.ca/fileadmin/Cours/MAT-2910/A-11/Solution_3.pdf|Université Laval]] avec une matrice $4*4$ {{ math:matrices:systeme_lineaire:directe:solution_3.pdf |Archive 14/10/2018}}
  
 ====Méthode de Gauss / Crout : LU (A non symétrique)==== ====Méthode de Gauss / Crout : LU (A non symétrique)====
-Valable pour une matrice carrée définie positive.+A doit être inversible (pas de pivot nul).
  
 Étape 1 : Étape 1 :
Ligne 69: Ligne 212:
 Si l'un des $l_{jj}$ vaut 0, il faut appliquer un pivot. Si l'un des $l_{jj}$ vaut 0, il faut appliquer un pivot.
  
-[[https://asi.insa-rouen.fr/enseignement/siteUV/ananum/04syslindirect.pdf|Résolution de systèmes linéaires par des méthodes directes : Gauss, LU]] {{ math:matrices:systeme_lineaire:directe:04syslindirect.pdf |Archive 14/10/2018}}+Implémentation du code en python :
  
-Exemple :+<code 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) 
 +</code>
  
-[[https://www2.mat.ulaval.ca/fileadmin/Cours/MAT-2910/A-11/Solution_3.pdf|Université Laval]] avec une matrice $4*4$ {{ math:matrices:systeme_lineaire:directe:solution_3.pdf |Archive 14/10/2018}} 
 ====Méthode de Crout : CDtC (A symétrique)==== ====Méthode de Crout : CDtC (A symétrique)====
 Valable pour une matrice carrée symétrique définie positive. Valable pour une matrice carrée symétrique définie positive.
Ligne 99: Ligne 258:
 Méthode : Méthode :
  
-[[https://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky_algorithm|The Cholesky algorithm]] {{ math:matrices:systeme_lineaire:directe:cholesky_decomposition_-_wikipedia.mhtml |Archive 17/10/2018}}+[[https://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky_algorithm|The Cholesky algorithm]] {{ :math:matrices:systeme_lineaire:directe:cholesky_decomposition_-_wikipedia_2020-04-28_10_55_30_pm_.html |Archive du 13/04/2020 le 28/04/2020}}
  
 Soit Soit
Ligne 264: Ligne 423:
 Exemple : Exemple :
  
-$A = \begin{pmatrix} 4   &  12 & -16 \\ 12  &  37 & -43 \\ -16 & -43 &  98 \end{pmatrix}$+$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}$ $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}$
  
-====QR Householder (A non symétrique)====+Implémentation en python : 
 + 
 +<code 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())) 
 +</code> 
 + 
 +====QR (A non symétrique) méthode Householder====
 $QR = A$ avec R matrice diagonale supérieure et Q matrice orthogonale. $QR = A$ avec R matrice diagonale supérieure et Q matrice orthogonale.
  
Ligne 292: Ligne 493:
 $$Q = \prod_{k=1}^{N-1}Q_i$$ $$Q = \prod_{k=1}^{N-1}Q_i$$
  
-[[https://fr.wikipedia.org/wiki/D%C3%A9composition_QR|Décomposition QR]] {{ :math:matrices:systeme_lineaire:directe:decomposition_qr_wikipedia.mhtml |Archive du 06/02/2019}}+[[https://fr.wikipedia.org/wiki/D%C3%A9composition_QR|Décomposition QR]] {{ :math:matrices:systeme_lineaire:directe:decomposition_qr_wikipedia_2020-04-28_10_55_29_pm_.html |Archive du 16/05/2019 le 28/04/2020}}
  
 [[http://metronu.ulb.ac.be/MATH-H-202/an_ch3.pdf|Factorisation QR et systèmes surdéterminés]] {{ :math:matrices:systeme_lineaire:directe:an_ch3.pdf |Archive du 06/02/2019}} [[http://metronu.ulb.ac.be/MATH-H-202/an_ch3.pdf|Factorisation QR et systèmes surdéterminés]] {{ :math:matrices:systeme_lineaire:directe:an_ch3.pdf |Archive du 06/02/2019}}
Ligne 375: Ligne 576:
 2/7 & -6/35 & 33/35 \\ 2/7 & -6/35 & 33/35 \\
 \end{pmatrix}$$ \end{pmatrix}$$
 +
 +====QR (A non symétrique) méthode Gram-Schmidt====
 +$\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}
 +$
 +
 +[[https://www.math.ucla.edu/~yanovsky/Teaching/Math151B/handouts/GramSchmidt.pdf|QR Decomposition with Gram-Schmidt]] {{ :math:matrices:systeme_lineaire:directe:gramschmidt.pdf |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}
 +$$
 +
 +====QR (A non symétrique) méthode Givens====
 +
 +$\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)}$.
 +
 +[[http://www.cs.nthu.edu.tw/~cherung/teaching/2008cs3331/chap4%20example.pdf|An Example of QR Decomposition]] {{ :math:matrices:systeme_lineaire:directe:chap4_example.pdf |Archive du 17/02/2019}}
 +
 +[[http://www.cs.rpi.edu/~flaherje/pdf/lin13.pdf|The QR Reduction]] {{ :math:matrices:systeme_lineaire:directe:lin13.pdf |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}
 +$$
 +
 =====Bibliographie===== =====Bibliographie=====
 <blockquote>[1] : DESTUYNDER, Philippe. Méthodes numériques pour l'ingénieur. Lavoisier, 2010, 248 pages. <blockquote>[1] : DESTUYNDER, Philippe. Méthodes numériques pour l'ingénieur. Lavoisier, 2010, 248 pages.
math/matrices/systeme_lineaire/directe.1550296306.txt.gz · Dernière modification : 2019/02/16 06:51 de root