arthur/ArtStudies
1
0
Fork 0
mirror of https://github.com/ArthurDanjou/ArtStudies.git synced 2026-01-14 15:54:13 +01:00
Code Issues Packages Projects Releases Wiki Activity

Move all files

Browse Source
This commit is contained in:
Arthur Danjou
2024-09-17 22:36:02 +02:00
parent 3012b8a505
commit 5093ee4d25
74 changed files with 85 additions and 1251 deletions
Download Patch File Download Diff File Expand all files Collapse all files

368
L3/Analyse Matricielle/TP1_Methode_de_Gauss.ipynb Normal file
View File

@@ -0,0 +1,368 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\\newcommand{\\nr}[1]{\\|#1\\|}\n",
"\\newcommand{\\RR}{\\mathbb{R}}\n",
"\\newcommand{\\N}{\\mathbb{N}}\n",
"$$\n",
"### MEU352 2023/2024 - Analyse numérique matricielle et optimisation\n",
"\n",
"# TP1 - Résolution de systèmes linéaires triangulaires. Méthode de Gauss.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercice 0. Manipulation de vecteurs et de matrices.\n",
"\n",
"On aura besoin des modules de python ``numpy`` et ``matplotlib.pyplot``. On peut les charger en exécutant les commandes\n",
"\n",
"``import numpy as np``\n",
"\n",
"``import matplotlib.pyplot as plt``\n",
"\n",
"(on désignera alors le module ``numpy`` par ``np`` et ``matplotlib.pyplot`` par ``plt``. \n",
"\n",
"**Q1.** Executez les commandes suivantes et affichez le résultat. Essayez de comprendre ce que vous avez obtenu.\n",
"\n",
"``\n",
"u = np.array([1,2,3,4,5])\n",
"v = np.array([[1,2,3,4,5]])\n",
"su=u.shape\n",
"sv=v.shape\n",
"ut = np.transpose(u)\n",
"vt = np.transpose(v)\n",
"vt2 = np.array([[1],[2],[3],[4],[5]])\n",
"A = np.array([[1,2,0,0,0],[0,0,2,3,1],[0,0,0,2,2],[0,0,0,0,1],[1,1,1,0,0]])\n",
"B = np.array([[1,2,3,4,5],[2,3,4,5,6],[3,4,5,6,7],[4,5,6,7,8],[5,6,7,8,9]])\n",
"d=np.diag(A)\n",
"dd=np.array([np.diag(A)])\n",
"dt=np.transpose(d)\n",
"ddt=np.transpose(dd)\n",
"Ad=np.diag(np.diag(A))``\n",
"\n",
"**Q2.** Même question pour les commandes suivantes.\n",
"\n",
"``u*v, u*vt, vt*u, u/v, u/vt, v/v, v/vt, np.vdot(u,v), np.vdot(u,vt)``\n",
"\n",
"``A*B, np.dot(A,B)``\n",
"\n",
"``np.dot(A,u), np.dot(A,v), np.dot(v,A), np.dot(A,vt), np.linalg.inv(A), np.dot(np.linalg(inv(A)),A))``"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1. 0. 0. 0. 0.]\n",
" [0. 1. 0. 0. 0.]\n",
" [0. 0. 1. 0. 0.]\n",
" [0. 0. 0. 1. 0.]\n",
" [0. 0. 0. 0. 1.]]\n"
]
}
],
"source": [
"import numpy as np\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"u = np.array([1,2,3,4,5])\n",
"v = np.array([[1,2,3,4,5]])\n",
"su=u.shape\n",
"sv=v.shape\n",
"ut = np.transpose(u)\n",
"vt = np.transpose(v)\n",
"vt2 = np.array([[1],[2],[3],[4],[5]])\n",
"A = np.array([[1,2,0,0,0],[0,2,0,0,0],[0,0,3,0,0],[0,0,0,4,0],[0,0,0,0,5]])\n",
"B = np.array([[1,2,3,4,5],[2,3,4,5,6],[3,4,5,6,7],[4,5,6,7,8],[5,6,7,8,9]])\n",
"d=np.diag(A)\n",
"dd=np.array([np.diag(A)])\n",
"dt=np.transpose(d)\n",
"ddt=np.transpose(dd)\n",
"Ad=np.diag(np.diag(A))\n",
"\n",
"print(np.dot(np.linalg.inv(A), A))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercice 1. Résolution d'un système linéaire triangulaire.\n",
"\n",
"Soit $A\\in\\mathcal{M}_n(\\RR)$ une matrice triangulaire inférieure inversible, de taille $n\\in\\N,$ et $b\\in\\RR^n$. Comme $A$ est triangulaire inférieure, on peut résoudre le système $Ax=b$ par une technique dite de *descente* : la solution $x=(x_1,\\dots,x_n)$ est obtenue en calculant successivement ses composantes $x_i$ par les formules\n",
"$$\n",
"\\begin{align}\n",
"x_1&=\\frac{b_1}{A_{11}}\\\\\n",
"x_2&=\\frac{b_2-A_{21}\\,x_1}{A_{22}}\\\\\n",
"&\\vdots\\\\\n",
"x_n&=\\frac{b_n-(A_{n1}\\,x_1+\\cdots +A_{n\\,n-1}\\,x_{n-1})}{A_{nn}}\\\\\n",
"&\\\\\n",
"&\\\\\n",
"\\bigg(\\,\\,x_i&=\\frac{b_i-(A_{i1}\\,x_1+\\cdots +A_{i\\,i-1}\\,x_{i-1})}{A_{ii}}\\,\\, \\bigg)\n",
"\\end{align}\n",
"$$\n",
"\n",
"**Q1.** Définir une fonction ``descente`` qui prend en argument une matrice $A$ triangulaire inférieure inversible et un vecteur $b$ et qui retourne la solution $x$ du système $Ax=b$. Tester votre fonction sur une matrice $A$ à coefficients aléatoires et un second membre $b$ tel que la solution $x$ de $Ax=b$ soit connue.\n",
"\n",
"**Q2.** Écrire la solution $x$ du système $Ax=b$ lorsque $A$ est cette fois-ci triangulaire supérieure, en fonction des coefficients de $A$ et de $b$, en résolvant successivement les équations depuis la dernière jusqu'à la première (on dit qu'on résout le système $Ax=b$ par *remontée*).\n",
"\n",
"**Q3.** Modifier votre fonction ``descente`` en une fonction que vous appelerez ``remonte_descente`` qui permet la résolution du système $Ax=b$ lorsque $A$ est triangulaire inférieure ou triangulaire supérieure. Votre fonction devra tester si la matrice $A$ est triangulaire supérieure ou inférieure.\n",
"\n",
"*Commandes python : essayez les commandes ``np.tril(A), np.triu(A), np.tril(A,k), np.triu(A,k)``, avec $k=1$ ou $k=-1$, et ``np.random.rand(n,n)``, avec $n\\in\\N$. La somme $(A_{i1}\\,x_1+\\cdots +A_{i\\,i-1}\\,x_{i-1})$ peut être vue comme un produit scalaire entre deux vecteurs, utiliser ``np.vdot`` pour le produit scalaire*. "
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"def remontee_descente(A, b):\n",
" x = 0 * b\n",
" n = len(b)\n",
" if np.allclose(A, np.triu(A)):\n",
" for i in range(n-1, -1, -1):\n",
" x[i] = (b[i] - np.dot(A[i,i+1:], x[i+1:])) / A[i,i]\n",
" elif np.allclose(A, np.tril(A)):\n",
" for i in range(n):\n",
" x[i] = (b[i] - np.dot(A[i,:i], x[:i])) / A[i,i]\n",
" else:\n",
" raise ValueError(\"A est ni triangulaire supérieure ni triangulaire inférieure\")\n",
" return x"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.1093356479670479e-31\n"
]
}
],
"source": [
"n = 5\n",
"A = np.random.rand(n, n)\n",
"A = np.tril(A) + np.eye(n) * np.linalg.norm(A)\n",
"xe = np.array([1] * n)\n",
"b = np.dot(A, xe)\n",
"x = remontee_descente(A, b)\n",
"\n",
"print(np.dot(x - xe, x-xe))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercice 2. La méthode de Gauss.\n",
"\n",
"\n",
"On suppose que $A$ est une matrice carrée inversible et qu'il est\n",
"possible d'appliquer la méthode de Gauss à la matrice $A$ et donc la transformer en une matrice triangulaire supérieure $U$\n",
"à coefficients diagonaux non nuls simplement en effectuant\n",
"successivement des opérations élémentaires sur les lignes du type\n",
"$L_i$ devient $L_i + \\beta L_j$. On suppose donc que les pivots de\n",
"la méthode de Gauss sont tous non nuls. \n",
"\n",
"**Q1.** Vérifier que l'algorithme suivant permet de transformer une matrice donnée $A$ en une matrice triangulaire supérieure $U$ par la méthode de Gauss :\n",
"\n",
"\n",
"```\n",
"U = A # on prend une copie qu'on écrasera\n",
"pour j = 0 à n-1\n",
" pour i = j + 1 à n - 1\n",
" beta = U(i,j)/U(j,j) # U(j,j) est le pivot\n",
" pour k = j à n -1\n",
" U(i, k) = U(i,k) - beta * U(j,k) # ligne i devient ligne i - beta * ligne j\n",
" fin k\n",
" fin i\n",
"fin j\n",
"retourner U\n",
"```\n",
"\n",
"**Q2.** Ecrire une fonction Python de la forme ```met_gauss(A)``` correspondant à cet algorithme.\n",
"\n",
"*Remarque : vous pouvez écrire les commandes $U(i, k) = U(i,k) - \\beta U(j,k)$, pour $k=j,\\dots,n-1$, sans utiliser de boucle sur $k$, en écrivant le vecteur $(U(i,j),\\dots,U(i,n-1))$ comme ``U[i,j:]``.*\n",
"\n",
"**Q3.** Appliquer cette fonction aux matrices\n",
"$$\n",
"A=\\left (\n",
"\\begin{array}{ccc}\n",
"9 & 8 & 6 \\\\\n",
"7 & 6 & 12 \\\\\n",
"9 & 3 & 9\n",
"\\end{array}\n",
"\\right )\n",
"\\qquad \\mbox{ et } \\qquad \n",
"B=\\left (\n",
"\\begin{array}{cccc}\n",
"11 & 8 & 3 & 13 \\\\\n",
" 2 & 12 & 7 & 10 \\\\\n",
" 3 & 3 & 17 & 13 \\\\\n",
" 11 & 2 & 12 & 7\n",
"\\end{array}\n",
"\\right )\n",
"$$\n",
"\n",
"Les réponses attendues sont respectivement\n",
"\\begin{equation*}\n",
"A=\\left (\n",
"\\begin{array}{ccc}\n",
"9 & 8 & 6 \\\\\n",
"0 & -0.2222222 & 7.3333333 \\\\\n",
"0 & 0 & -162\n",
"\\end{array}\n",
"\\right )\\qquad\n",
" \\mbox{ et }\\qquad \n",
"B=\\left (\n",
"\\begin{array}{cccc}\n",
"11 & 8 & 3 & 13 \\\\\n",
" 0 & 10.5454545 & 6.45454545 & 7.63636364\\\\\n",
" 0 & 0 & 15.6810345 & 8.86206897 \\\\\n",
" 0 & 0 & 0 & -8.81693238\n",
"\\end{array}\n",
"\\right )\n",
"\\end{equation*}\n",
"\n",
"**Q4.** Adapter votre fonction ```met_gauss(A)``` en une fonction ```met_gauss_sys(A,b)``` de façon à que l'on puisse l'utiliser pour résoudre un système $Ax=b$, avec $b\\in\\RR^n$ donné. Pour cela, il faut le long de la méthode de Gauss faire les mêmes opérations sur la matrice $A$ et sur le second membre $b$. Cette fonction retournera la solution $x$ du système $Ax=b$ en écrivant le système triangulaire équivalent obtenu par la méthode de Gauss, et en résolvant ce système triangulaire avec la fonction ```remonte_descente```. La tester avec une matrice $A$ et un vecteur $b$ aléatoires par exemple."
]
},
{
"cell_type": "code",
"execution_count": 124,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"def met_gauss(A):\n",
" U = A\n",
" n = len(A)\n",
" for j in range(n):\n",
" for i in range(j+1, n):\n",
" beta = U[i,j]/U[j,j]\n",
" U[i,j:] = U[i,j:] - beta * U[j, j:]\n",
" return U"
]
},
{
"cell_type": "code",
"execution_count": 161,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"def met_gauss_sys(A, b):\n",
" n, m = A.shape\n",
" if n != m:\n",
" raise ValueError(\"Erreur de dimension : A doit etre carré\")\n",
" if n != b.size:\n",
" raise valueError(\"Erreur de dimension : le nombre de lignes de A doit être égal au nombr ede colonnes de b\")\n",
" U = np.zeros((n, n+1))\n",
" U = A\n",
" V = b\n",
" for j in range(n):\n",
" for i in range(j+1, n):\n",
" beta = U[i,j]/U[j,j]\n",
" U[i,j:] = U[i,j:] - beta * U[j, j:]\n",
" V[i] = V[i] - beta * V[j]\n",
" return remontee_descente(U, V)"
]
},
{
"cell_type": "code",
"execution_count": 165,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.3096323621833204e-32\n"
]
}
],
"source": [
"n = 5\n",
"A = np.random.rand(n, n) + float(n) * np.eye(n)\n",
"b = np.random.rand(n)\n",
"x = met_gauss_sys(A, b)\n",
"\n",
"print(np.dot(b - A.dot(x), b - A.dot(x)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## La méthode de Gauss avec stratégie de pivot partiel.\n",
"\n",
"En pratique, pour des questions de stabilité numérique, on a intérêt à choisir à l'étape $j$ un pivot $A_{k,j}$, avec $k\\geq j$ tel que $|A_{k,j}|$ est maximal (car cela signifie diviser par la quantité la plus grande possible). À l'étape $j$ de la méthode de Gauss, on commence alors par choisir $p$ tel que $|A_{p,j}|=\\max_{k\\geq j}|A_{k,j}|$ et on échange les lignes $p$ et $j$ de $A$.\n",
"\n",
"**Q5.** Créer une fonction ```met_gauss_pivot(A,b)``` qui permet la résolution du système $Ax=b$ en utilisant cette stratégie de choix de pivot. La tester avec le même exemple que dans la question précédente.\n",
"\n",
"**Q6. (Comparaison des deux méthodes).** Pour $n=10,20,30,\\dots,200$ :\n",
"* construire une matrice $A\\in\\mathcal{M}_n(\\RR)$ aléatoire, un vecteur $x_{ex}\\in\\RR^n$ aléatoire et calculer $b=Ax_{ex}$ ;\n",
"* Résoudre le système $Ax=b$ (dont la solution est $x=x_{ex}$) par la méthode de Gauss avec et sans choix de pivot ;\n",
"* Calculer la norme $\\|x-x_{ex}\\|$ pour chacune des méthodes.\n",
"\n",
"Comparer les résultats obtenus pour les deux méthodes. Vous pouvez représenter $\\|x-x_{ex}\\|$, ou, ce qu'est mieux, $\\mathrm{log}(\\|x-x_{ex}\\|)$ en fonction de la taille $n$ de la matrice.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.8"
}
},
"nbformat": 4,
"nbformat_minor": 4
}