math:matrices:systeme_lineaire:directe
Différences
Ci-dessous, les différences entre deux révisions de la page.
Les deux révisions précédentesRévision précédenteProchaine révision | Révision précédente | ||
math:matrices:systeme_lineaire:directe [2019/03/17 19:05] – [Principe de la décomposition en deux matrices LU] : ajout du code python root | math: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 30: | Ligne 29: | ||
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 | + | Implémentation en python |
<code python> | <code python> | ||
import numpy as np | import numpy as np | ||
Ligne 85: | Ligne 83: | ||
fsub(np.array([[1, | fsub(np.array([[1, | ||
</ | </ | ||
+ | |||
+ | 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(' | ||
+ | | ||
+ | N = len(A) | ||
+ | B = np.identity(N, | ||
+ | for i in range(0, N): | ||
+ | | ||
+ | indice = np.argmax(np.abs(A[i:, | ||
+ | |||
+ | if (A[i+indice, | ||
+ | A.fill(np.NaN) | ||
+ | # | ||
+ | return (0, A) | ||
+ | | ||
+ | A[[i, | ||
+ | B[[i, | ||
+ | | ||
+ | for j in range(i+1, N): | ||
+ | multiplicateur = A[j][i]/ | ||
+ | B[j,:] -= B[i,: | ||
+ | A[j][i:] -= A[i][i: | ||
+ | | ||
+ | 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) | ||
+ | </ | ||
+ | |||
+ | |||
+ | ===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' | ||
+ | |||
+ | 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, | ||
+ | if len(A) != len(A[0]): | ||
+ | raise ValueError(' | ||
+ | if len(A) != len(b): | ||
+ | raise ValueError(' | ||
+ | | ||
+ | N = len(A) | ||
+ | for i in range(0, N): | ||
+ | | ||
+ | indice = np.argmax(np.abs(A[i:, | ||
+ | |||
+ | if (A[i+indice, | ||
+ | # | ||
+ | # Si c'est aucune solution ou une infinité (cf. tp3). | ||
+ | continue | ||
+ | | ||
+ | A[[i, | ||
+ | b[i], b[i+indice] = b[i+indice], | ||
+ | | ||
+ | for j in range(i+1, N): | ||
+ | multiplicateur = A[j][i]/ | ||
+ | b[j] = b[j] - b[i]*multiplicateur | ||
+ | A[j][i:] = A[j][i:] - A[i][i: | ||
+ | </ | ||
+ | |||
+ | Triangularisation : | ||
+ | <code python> | ||
+ | def gauss(A, b): | ||
+ | pivotgauss(A, | ||
+ | (x, sol) = bsub(A, b) | ||
+ | return (x, sol) | ||
+ | </ | ||
+ | |||
+ | [[https:// | ||
+ | |||
+ | Exemple : | ||
+ | |||
+ | [[https:// | ||
+ | |||
====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 124: | 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:// | + | 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, | ||
+ | L[i, | ||
+ | U[i,j:] -=L[i, | ||
+ | |||
+ | return (L,U,1) | ||
+ | </ | ||
- | [[https:// | ||
====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 154: | Ligne 258: | ||
Méthode : | Méthode : | ||
- | [[https:// | + | [[https:// |
Soit | Soit | ||
Ligne 322: | 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}$ | $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, | ||
+ | L[0,0] = cmath.sqrt(A[0, | ||
+ | L[1:,0] = A[1:, | ||
+ | |||
+ | for j in range(1,n): | ||
+ | if (not np.isreal(A[j, | ||
+ | # Matrice A non hermitienne | ||
+ | L.fill(np.nan) | ||
+ | return (L,-1) | ||
+ | |||
+ | carre = A[j,j] - np.dot(L[j,: | ||
+ | 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, | ||
+ | if (A[i,j] != np.conjugate(A[j, | ||
+ | # Matrice A non hermitienne | ||
+ | L.fill(np.nan) | ||
+ | return (L,0) | ||
+ | L[i, | ||
+ | |||
+ | return (L,1) | ||
+ | |||
+ | |||
+ | A4 = np.array( | ||
+ | [[ 4, 12, -16], | ||
+ | [ 12, 37, -43+2j], | ||
+ | [-16, -43-2j, | ||
+ | (L, sol) = cholesky(A4) | ||
+ | print (" | ||
+ | print (" | ||
+ | print (L.dot(L.transpose().conjugate())) | ||
+ | </ | ||
====QR (A non symétrique) méthode Householder==== | ====QR (A non symétrique) méthode Householder==== | ||
Ligne 347: | Ligne 493: | ||
$$Q = \prod_{k=1}^{N-1}Q_i$$ | $$Q = \prod_{k=1}^{N-1}Q_i$$ | ||
- | [[https:// | + | [[https:// |
[[http:// | [[http:// |
math/matrices/systeme_lineaire/directe.1552845949.txt.gz · Dernière modification : 2019/03/17 19:05 de root