Différences

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

--- math:matrices:systeme_lineaire:directe [2019/03/17 20:54] – [Méthode de Gauss / Crout : LU (A non symétrique)] : ajout pivot root
+++ math:matrices:systeme_lineaire:directe [2020/04/28 23:00] (Version actuelle) – mhtml -> html root
@@ Ligne 1: / Ligne 1: @@
 =====Systèmes linéaires=====
 ====Résolution====
@@ Ligne 30: / Ligne 29: @@
 Puis on applique le même principe avec la matrice supérieure $U$ pour déterminer le vecteur $x$.
-Implémentation en python :
+Implémentation en python pour la résolution d'un système $Ax=b$
 <code python>
 import numpy as np
@@ Ligne 85: / Ligne 83: @@
 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>
-====Méthode de Gauss / Crout : LU (A non symétrique)====
-Valable pour une matrice carrée définie positive.
-Étape 1 :
+Implémentation en python pour inverser une matrice. On applique le même algorithme ci-dessus mais avec B matrice identité.
-$$
+<code python>
-\begin{align*}
+import numpy as np
-&j = 1 \\
+import math
-&i = 1 \quad &l_{jj} &= a_{jj} \\
+import cmath
-&\forall i \in [2,N] \quad &u_{ji} &= \frac{a_{ji}}{l_{jj}}
-\end{align*}
-$$
-Étape 2 :
+def pivotgauss(A):
-$$
+    if len(A) != len(A[0]):
-\begin{align}
+        raise ValueError('La matrice n''est pas carrée.')
-\forall j &\in [2,N-1] : \\
-i &= 1 \quad & l_{ji} &= a_{ji} \\
+    N = len(A)
-\forall &i \in [2, j] \quad & l_{ji} &= a_{ji} - \sum_{k=1}^{i-1}l_{jk}u_{ki} \\
+    B = np.identity(N, dtype=np.float64)
-\forall &i \in [j+1, N] \quad & u_{ji} &= \frac{a_{ji} - \sum_{k=1}^{j-1}l_{jk}u_{ki}}{l_{jj}}
+    for i in range(0, N):
-\end{align}
-$$
+        indice = np.argmax(np.abs(A[i:,i]))
-Étape 3:
+        if (A[i+indice, i+indice] == 0):
-$$
+            A.fill(np.NaN)
-\begin{align}
+#           C'est inacceptable car on inverse la matrice (cf. tp2).
-\forall i &\in [1,N] : \\
+            return (0, A)
-j &= N \\
-l_{ji} &= a_{ji} - \sum_{k=1}^{i-1}l_{jk}u_{ki} \\
+        A[[i,i+indice]] = A[[i+indice,i]]
-\end{align}
+        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>
-Complexité : $\frac{N^3}{3}$ multiplications, $N$ divisions.
-[[math:matrices:systeme_lineaire:directe#bibliographie|[1]]] §2.2 L'algorithme de Gauss.
-[[https://www.code-aster.org/V2/doc/v11/fr/man_r/r6/r6.02.01.pdf|À propos des méthodes de décomposition de type GAUSS]] {{ math:matrices:systeme_lineaire:directe:r6.02.01.pdf |Archive 14/10/2018}}
-Si l'un des $l_{jj}$ vaut 0, il faut appliquer un pivot.
 ===Pivot===
@@ Ligne 173: / Ligne 172: @@
 [[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)====
+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.
+[[math:matrices:systeme_lineaire:directe#bibliographie|[1]]] §2.2 L'algorithme de Gauss.
+[[https://www.code-aster.org/V2/doc/v11/fr/man_r/r6/r6.02.01.pdf|À propos des méthodes de décomposition de type GAUSS]] {{ math:matrices:systeme_lineaire:directe:r6.02.01.pdf |Archive 14/10/2018}}
+Si l'un des $l_{jj}$ vaut 0, il faut appliquer un pivot.
+Implémentation du code en python :
+<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>
 ====Méthode de Crout : CDtC (A symétrique)====
 Valable pour une matrice carrée symétrique définie positive.
@@ Ligne 198: / Ligne 258: @@
 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
@@ Ligne 366: / Ligne 426: @@
 $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 :
+<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====
@@ Ligne 391: / Ligne 493: @@
 $$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}}