{ "cells": [ { "cell_type": "markdown", "id": "7c8c115e-dc7b-41b3-907a-700b3503fc61", "metadata": {}, "source": [ "## MT09 - TP3 - Automne 2024\n", "### Factorisation de Choleski, méthodes itératives, méthode de Newton" ] }, { "cell_type": "markdown", "id": "9c071e68-a309-4a2c-96ad-ecc661d6891a", "metadata": {}, "source": [ "### 1. Programmation de la factorisation de Choleski sous forme vectorielle\n", "Pour une matrice $A$ symétrique définie positive, programmez la factorisation de Choleski \n", "\n", "\\begin{align*}\n", "& b_{11} = \\sqrt{a_{11}},\\quad b_{i1}=\\frac{a_{i1}}{b_{11}}, \\ i>1, \\\\[1.5ex]\n", "& b_{jj} = \\sqrt{a_{jj}-\\sum_{k=1}^{j-1} b_{jk}^2},\\ j>1,\\quad\n", "b_{ij} =\\frac{1}{b_{jj}} \\left( a_{ij}-\\sum_{k=1}^{j-1} b_{ik} b_{jk}\\right),\\ 10$ et consomme localement à taux\n", "$k_i x_i$, $k_i>0$. Par ailleurs, pour réguler les stocks, entre deux régions $i$ et $j$, il y a des échanges de $j$ vers $i$ selon le flux d'échange $\\phi_{ij}=-k(x_j-x_i)$. On considère un réseau d'échanges\n", "étoilé plein.\n", "\n", "Si le système atteint son régime d'équilibre, on a alors un bilan conservatif échange-production-consommation de matière $m$ qui décrit pour chaque région $i$, $i\\in\\{1,...,N\\}$:\n", "\n", "$$\n", " \\sum_{j\\neq i} -k (x_j-x_i) = p_i -k_i x_i.\n", "$$\n", "\n", "Pour l'illustration on considère $N=5$ régions. On veut déterminer les stocks\n", "$x_i$ de matière première pour chaque région numéro $i$. Ecrire le système linéaire d'inconnues\n", "$(x_1,...,x_5)$.\n" ] }, { "cell_type": "markdown", "id": "7a2a0e00-d071-4c40-90b9-fecd00405d29", "metadata": { "jupyter": { "source_hidden": true }, "tags": [] }, "source": [ "Réponse : \n", "\n", "$$\n", "k\n", "\\begin{pmatrix}\n", "(\\frac{k_1}{k}+4) & -1 & -1 & -1 & -1 \\\\\n", "-1 & (\\frac{k_2}{k}+4) & -1 & -1 & -1 \\\\\n", "-1 & -1 & (\\frac{k_3}{k}+4) & -1 & -1 \\\\\n", "-1 & -1 & -1 & (\\frac{k_4}{k}+4) & -1 \\\\\n", "-1 & -1 & -1 & -1 & (\\frac{k_5}{k}+4)\n", "\\end{pmatrix}\n", "\\begin{pmatrix}x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4 \\\\ x_5 \\end{pmatrix}\n", "= \n", "\\begin{pmatrix} p_1 \\\\ p_2 \\\\ p_3 \\\\ p_4 \\\\ p_5 \\end{pmatrix}\n", "$$" ] }, { "cell_type": "markdown", "id": "3ccfed36-fa0a-4b28-a1b8-b4c378c4ff79", "metadata": { "tags": [] }, "source": [ "Que pouvez-vous dire de la matrice $A$ du système linéaire ?\n", "\n", "Pour l'application numérique, utiliser les paramètres $k=1$, \n", "$k_1= 0.1$, $k_2=k_3=k_4=0.2$, $k_5=3$, $p_1=8$, $p_2=p_3=1$, $p_4=p_5=2$. \n", "Ecrire en python la matrice $A$ du système linéaire ainsi que le second membre $b$." ] }, { "cell_type": "code", "execution_count": null, "id": "90c4f711-6ded-465b-a0e1-3e222a6aab61", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "aa37e842-5869-4d0b-84bc-b1570087f65f", "metadata": {}, "source": [ "#### Algorithme de Gauss-Seidel\n", "\n", "Comme la matrice est à diagonale strictement dominante, on peut utiliser la méthode de Gauss-Seidel\n", "pour résolution itérative du système linéaire $Ax=b$ (la convergence est garantie dans ce cas). Ecrire une fonction python\n", "```\n", "x, k = GS(A, b, x0, kmax, tol)\n", "```\n", "qui implémente de Gauss-Seidel pour la matrice $A$ le second membre $b$, le vecteur initial $x^0$, un nombre d'itérations maximale $k_{max}$ et utilisant le test de convergence\n", "\n", "$$\n", "\\frac{\\| Ax^k - b \\|}{\\| A x^0 - b \\|} < \\text{tol}\n", "$$\n", "avec une tolérance *tol*.\n", "\n", "Vérifier votre fonction ```GS()``` sur le système test $2\\times 2$ avec\n", "\n", "$$\n", "A = \\begin{pmatrix} 2 & -1 \\\\ -1 & 2\\end{pmatrix}, \\quad\n", "b = \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix}.\n", "$$\n", "\n", "On prendre $k_{max}=100$, $x^{(0)}=(0,0)^T$, $tol=10^{-12}$." ] }, { "cell_type": "code", "execution_count": null, "id": "f3e28dce-0d50-4907-926b-99cbc96802d5", "metadata": { "tags": [] }, "outputs": [], "source": [ "import numpy.linalg as la\n", "\n", "def GS(A, b, x0, kmax, tol):\n", " k = 0\n", " n = x0.size\n", " xk = np.copy(x0)\n", " # A COMPLETER\n", " return xk, k\n", "\n", "# Verif\n", "print(\"Vérification Gauss-Seidel ...\\n\")\n", "# ..." ] }, { "cell_type": "markdown", "id": "d23a28f5-93c4-4bf3-a303-71e1cdf8caed", "metadata": {}, "source": [ "Résoudre le système linéaire $Ax=b$ du système production-consommation à l'aide de la fonction ```GS()```.\n", "On utilisera $x^0=(0,...,0)$, $k_{max}=100$, $tol=10^{-10}$." ] }, { "cell_type": "code", "execution_count": null, "id": "0f145e9b-cc57-4568-80c5-b353a8ae2104", "metadata": { "tags": [] }, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "e971acbf-213f-49d9-a934-8dc77b717e4d", "metadata": { "tags": [] }, "source": [ "### 3. Méthode de Newton-Raphson (optionnel et hors évaluation, pour les élèves intéressé(e)s ou finissant en avance)\n", "La méthode de Newton-Raphson permet de résoudre numériquement un système d'équations algébriques non linéaires écrit sous forme condensée vectorielle\n", "\n", "$$\n", "\\mathbf{F}(\\mathbf{x}) = \\mathbf{0}.\n", "$$\n", "\n", "La méthode itérative s'écrit\n", "\n", "$$\n", "\\mathbf{x}^{(k+1)} = \\mathbf{x}^{(k)} - [J(\\mathbf{x}^{(k)})]^{-1}\\, \\mathbf{F}(\\mathbf{x}^{(k)})\n", "$$\n", "\n", "en partant de $\\mathbf{x}^{(0)}$ donné (initialization). La matrice\n", "$J(\\mathbf{x}^{(k)})$ est la matrice jacobienne de $\\mathbf{F}$ au point \n", "$\\mathbf{x}^{(k)}$.\n", "\n", "En pratique, on procède en deux étapes: on résout d'abord le système linéaire\n", "\n", "$$\n", "J(\\mathbf{x}^{(k)}) \\, \\mathbf{w}^{(k)} = -\\mathbf{F}(\\mathbf{x}^{(k)}),\n", "$$\n", "\n", "puis on fait l'affectation\n", "\n", "$$\n", "\\mathbf{x}^{(k+1)} = \\mathbf{x}^{(k)} + \\mathbf{w}^{(k)}.\n", "$$\n", "\n", "Ecrire une fonction python\n", "```\n", "xout, kout, Fout = newtonRaphson(x0, foncjac, kmax, tol)\n", "```\n", "qui résout l'équation $F(x)=0$ avec la méthode de Newton-Raphson avec ```x0``` l'initial guess, la fonction python\n", "```\n", "F, J = foncjac(x)\n", "``` \n", "retourne la fonction $F$ ainsi que sa matrice jacobienne $J(\\mathbf{x})=DF(\\mathbf{x})$ au point trouvé, $k_{max}$ est le nombre max d'itérations de Newton, $tol$ est la tolérance de convergence : on fera le test\n", "\n", "$$\n", "\\frac{\\|\\mathbf{F}(\\mathbf{x}^{(k)})\\|}{\\|\\mathbf{F}(\\mathbf{x}^{(0)})\\|} > tol\n", "$$\n", "\n", "(si oui et si $k