mirror of
https://github.com/ArthurDanjou/ArtStudies.git
synced 2026-01-14 15:54:13 +01:00
Arthur DANJOU
d5a6bfd339
- Added missing commas in various print statements and function calls for better syntax. - Reformatted code to enhance clarity, including breaking long lines and aligning parameters. - Updated function signatures to use float type for sigma parameters instead of int for better precision. - Cleaned up comments and documentation strings for clarity and consistency. - Ensured consistent formatting in plotting functions and data handling.
387 lines
14 KiB
Plaintext
387 lines
14 KiB
Plaintext
{
|
|
"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",
|
|
"\n",
|
|
"%matplotlib inline\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(\n",
|
|
" [\n",
|
|
" [1, 2, 0, 0, 0],\n",
|
|
" [0, 2, 0, 0, 0],\n",
|
|
" [0, 0, 3, 0, 0],\n",
|
|
" [0, 0, 0, 4, 0],\n",
|
|
" [0, 0, 0, 0, 5],\n",
|
|
" ],\n",
|
|
")\n",
|
|
"B = np.array(\n",
|
|
" [\n",
|
|
" [1, 2, 3, 4, 5],\n",
|
|
" [2, 3, 4, 5, 6],\n",
|
|
" [3, 4, 5, 6, 7],\n",
|
|
" [4, 5, 6, 7, 8],\n",
|
|
" [5, 6, 7, 8, 9],\n",
|
|
" ],\n",
|
|
")\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(\n",
|
|
" \"Erreur de dimension : le nombre de lignes de A doit être égal au nombr ede colonnes de b\",\n",
|
|
" )\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
|
|
}
|