{ "cells": [ { "cell_type": "markdown", "id": "e5c5ff2b", "metadata": {}, "source": [ "# Cours ENPC - Pratique du calcul scientifique" ] }, { "cell_type": "markdown", "id": "5d6c8224", "metadata": {}, "source": [ "## Systèmes linéaires" ] }, { "cell_type": "code", "execution_count": null, "id": "24baa415", "metadata": {}, "outputs": [], "source": [ "# Bibliothèques et configuration par défaut utiles pour le TD\n", "# Cela n'empêche nullement d'installer et de faire appel à d'autres bibliothèques si nécessaire\n", "\n", "using LinearAlgebra, Polynomials, Plots, LaTeXStrings\n", "\n", "Plots.default(fontfamily=\"Computer Modern\",\n", " titlefontsize=20,\n", " xlabelfontsize=20,\n", " ylabelfontsize=20,\n", " legendfontsize=12,\n", " xtickfontsize=12,\n", " ytickfontsize=12,\n", " framestyle=:box,\n", " linewidth=2,\n", " label=nothing,\n", " grid=true)" ] }, { "cell_type": "markdown", "id": "dd49f15f", "metadata": {}, "source": [ "### Exercice sur la factorisation de Cholesky\n", "\n", "L'objectif de cet exercice est de proposer un algorithme permettant de réaliser la factorisation de Cholesky d'une matrice réelle définie positive $\\mathsf{A}\\in\\mathbb{R}^{n×n}$ i.e. d'identifier la matrice réelle triangulaire inférieure $\\mathsf{C}$ telle que :\n", "\n", "$$\n", "\\tag{1}\n", "\\mathsf{A}=\\mathsf{C}\\mathsf{C}^T\n", "$$\n", "\n", "1. Construire une fonction qui prend comme argument une matrice `A` et qui renvoie la matrice `C`.\n", " Pour calculer $\\mathsf{C}$, s'appuyer sur une identification successive des colonnes de `C` en commençant par l'élément diagonal puis le reste de la colonne,\n", " en comparant les deux membres de (1).\n", "\n", "1. Construire une fonction permettant de générer une matrice définie positive aléatoire dans $\\mathbb{R}^{n×n}$. On pourra par exemple générer une matrice aléatoire `B` avec `rand(n,n)` puis calculer `B'B+I` où l'identité `I` est l'objet `UniformScaling` de la bibliothèque `LinearAlgebra`. (*L'ajout de l'identité permet de s'affranchir des éventuelles valeurs propres trop proches de `0`.*)\n", "\n", "1. Tester les algorithmes précédents par un code du type\n", " ```julia\n", " n = 1000\n", " A = generate_defpos_matrix(n)\n", " C = cholesky(A)\n", " norm(C*C' - A, Inf)\n", " ```\n", "\n", "1. On suppose maintenant que la matrice réelle définie positive `A` est à largeur de bande `b` supposée beaucoup plus petite que `n`. Réécrire la fonction de décomposition de Cholesky en exploitant la largeur de bande.\n", "\n", "1. Construire une fonction permettant de générer une matrice définie positive aléatoire à largeur de bande donnée `b`. On pourra commencer par générer une matrice aléatoire triangulaire inférieure `B` à largeur de bande `b` puis renvoyer `B'B+I`.\n", "\n", "1. Tester les algorithmes des deux questions précédentes. On prendra par exemple `n=1000` et `b=4`.\n", "\n", "1. *Optionnel :* réaliser une étude de la complexité des algorithmes non optimisé et optimisé pour matrices à largeur de bande donnée en traçant en échelle double log les temps de calcul pour différentes valeurs de `n`.\n", "\n", " *Suggestions :*\n", " - Prendre pour `n` des puissances de `2` de `2⁷` à `2¹⁰`\n", "\n", " - Utiliser la macro `@elapsed` pour enregistrer le temps de calcul d'un appel en faisant une moyenne sur un certain nombre de réalisations (par exemple `10`)\n", "\n", " - Exploiter la fonction `fit` de [`Polynomials.jl`](https://github.com/JuliaMath/Polynomials.jl) pour trouver la puissance de la complexité en la supposant de la forme $αn^\\beta$." ] }, { "cell_type": "code", "execution_count": null, "id": "ea4f4085", "metadata": {}, "outputs": [], "source": [ "function cholesky(A)\n", " # Your code comes here\n", " return C\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "ac9076f7", "metadata": {}, "outputs": [], "source": [ "function generate_defpos_matrix(n)\n", " B = rand(n,n)\n", " return B'B+I\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "25f0de98", "metadata": {}, "outputs": [], "source": [ "n = 1000\n", "A = generate_defpos_matrix(n)\n", "C = cholesky(A)\n", "norm(C*C' - A, Inf)" ] }, { "cell_type": "code", "execution_count": null, "id": "ceebf301", "metadata": {}, "outputs": [], "source": [ "function cholesky_banded(A, b)\n", " # Your code comes here\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "44b06979", "metadata": {}, "outputs": [], "source": [ "function generate_banded_matrix(n, b)\n", " C = [j≤i≤j+b ? rand() : zero(Float64) for i in 1:n, j in 1:n]\n", " return C*C'+I\n", "end\n", "\n", "n = 10 ; b = 4\n", "A = generate_banded_matrix(n,b)\n", "C = cholesky_banded(A,b)\n", "norm(C*C' - A, Inf)" ] }, { "cell_type": "code", "execution_count": null, "id": "d08e9477", "metadata": {}, "outputs": [], "source": [ "nb_samples = 10\n", "mean(x) = sum(x)/length(x)\n", "plot(xaxis=:log10, yaxis=:log10, xlabel=\"n\", ylabel=\"CPU time\", legend=:topleft)\n", "tf = Float64[] ; tb = Float64[]\n", "tn = 2 .^(7:10)\n", "for n in tn\n", " A = generate_banded_matrix(n, b)\n", " println(n)\n", " push!(tf, mean(@elapsed cholesky(A) for _ in 1:nb_samples))\n", " push!(tb, mean(@elapsed cholesky_banded(A,b) for _ in 1:nb_samples))\n", "end\n", "Pf = fit(log10.(tn),log10.(tf),1) ; af = round(coeffs(Pf)[2]; digits=2)\n", "plot!(tn, tf, marker=:o, label=\"Algo matrice pleine \"*latexstring(\"n^{$(af)}\"))\n", "ntn = 2 .^(11:12)\n", "for n in ntn\n", " A = generate_banded_matrix(n, b)\n", " println(n)\n", " push!(tb, mean(@elapsed cholesky_banded(A,b) for _ in 1:nb_samples))\n", "end\n", "append!(tn,ntn)\n", "Pb = fit(log10.(tn),log10.(tb),1) ; ab = round(coeffs(Pb)[2]; digits=2)\n", "plot!(tn, tb, marker=:diamond, label=\"Algo matrice bande \"*latexstring(\"n^{$(ab)}\"))" ] }, { "cell_type": "markdown", "id": "6b5584fc", "metadata": {}, "source": [ "### Exercice: itération de Richardson\n", "\n", "Considérer le système linéaire suivant:\n", "$$\n", " \\mathsf A \\mathbf x :=\n", " \\begin{pmatrix}\n", " 3 & 1 \\\\ 1 & 3\n", " \\end{pmatrix}\n", " \\begin{pmatrix}\n", " x_1 \\\\\n", " x_2\n", " \\end{pmatrix}\n", " =\n", " \\begin{pmatrix}\n", " 9 \\\\\n", " 11\n", " \\end{pmatrix} =: \\mathbf b.\n", "$$\n", "\n", " 1. Illustrer à l'aide de la fonction `Plots.contourf` les lignes de niveau de la fonction\n", " $$\n", " f(\\mathbf x) = \\frac{1}{2} \\mathbf x^T \\mathsf A \\mathbf x - \\mathbf b^T \\mathbf x.\n", " $$\n", "\n", " 1. Implémenter l'itération de Richardson avec $\\omega = 0.1$ pour résoudre le système,\n", " et illustrer les itérations au dessus des lignes de niveau de la fonction $f$.\n", " Comme critère d'arrêt, utiliser\n", " $$\n", " \\lVert \\mathsf A \\mathbf x - \\mathbf b \\lVert \\leq \\varepsilon \\lVert \\mathbf b \\lVert,\n", " $$\n", " evec $\\varepsilon$ suffisamment petit.\n", "\n", " 1. Faire un plot de l'erreur $e_k := \\lVert \\mathbf x^{(k)} - \\mathbf x_* \\rVert$ en fonction de $k$,\n", " en utilisant une échelle linéaire pour l'axe des abcisses et une échelle logarithmique pour l'axe des ordonnées,\n", " gràce à l'argument `yscale=:log` passé à la fonction `Plots.plot`.\n", "\n", " 1. En utilisant `Polynomials.fit`, calculer une approximation du type\n", " $$\n", " \\log(e_k) \\approx a + bk \\qquad \\Leftrightarrow \\qquad e_k \\approx \\exp(a) \\times \\exp(b)^k.\n", " $$\n", " Comparer la valeur de $\\exp(b)$ au rayon spectral $\\rho(\\mathsf A - \\omega \\mathsf I)$ et expliquer.\n", "\n", " 1. Calculer le paramètre $\\omega$ optimal et refaire le plot de la décroissance de l'erreur dans ce cas." ] }, { "cell_type": "code", "execution_count": null, "id": "843a92af", "metadata": {}, "outputs": [], "source": [ "using LaTeXStrings\n", "using LinearAlgebra\n", "using Plots\n", "using Polynomials\n", "\n", "A = [3. 1.; 1. 3.]\n", "b = [11.; 9.]\n", "sol = A\\b\n" ] }, { "cell_type": "markdown", "id": "5fe2c6db", "metadata": {}, "source": [ "### Exercice: itération de Gauss—Seidel\n", "\n", "On s'intéresse dans cet exercice à la résolution de l'équation de Poisson en dimension 1 avec conditions de Dirichlet homogènes:\n", "$$\n", "- u''(x) = b(x), \\qquad u(0) = u(1) = 0, \\qquad b(x) := \\exp(x).\n", "$$\n", "Une discrétisation de cette équation sur un maillage uniforme par la méthode des différences finies conduit au système linéaire\n", "$$\n", "\\frac{1}{h^2}\n", "\\begin{pmatrix}\n", " 2 & -1 \\\\\n", " -1 & 2 & -1 \\\\\n", " & -1 & 2 & -1 \\\\\n", " & & \\ddots & \\ddots & \\ddots & \\\\\n", " & & & -1 & 2 & -1 \\\\\n", " & & & & -1 & 2 \\\\\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", " u_1 \\\\\n", " u_2 \\\\\n", " u_3 \\\\\n", " \\vdots \\\\\n", " u_{n-2} \\\\\n", " u_{n-1}\n", "\\end{pmatrix}\n", "=\n", "\\begin{pmatrix}\n", " b(x_1) \\\\\n", " b(x_2) \\\\\n", " b(x_3) \\\\\n", " \\vdots \\\\\n", " b(x_{n-2}) \\\\\n", " b(x_{n-1})\n", "\\end{pmatrix},\n", "\\qquad\n", "h := \\frac{1}{n},\n", "\\qquad\n", "x_i := ih.\n", "$$\n", "où $h$ est l'espace entre les points de discrétisation et $(u_1, u_2, \\cdots, u_{n-1})$ sont les valeurs recherchées de la fonction inconnue $u$ aux points intérieurs du domaine $[0, 1]$.\n", "\n", "1. Calculer la solution du système pour $n = 50$ en utilisant la méthode `\\` de Julia.\n", "\n", " *Indication*: pour construire la matrice du système linéaire,\n", " on pourra utiliser `LinearAlgebra.SymTridiagonal` ainsi que la fonction `fill`.\n", "\n", "1. Implémenter la méthode de Gauss-Seidel (ou sa généralisation, la méthode de relaxation) afin de résoudre le système linéaire pour $n = 50$,\n", " en n'utilisant cette fois pas les fonctions `\\` et `inv` de Julia ni aucune bibliothèque logicielle.\n", " On initialisera l'itération à $\\mathbf x_0 = \\mathbf 0$.\n", "\n", "1. Vérifier sur un graphe que la solution approchée obtenue est proche de la solution exacte,\n", " qui est donnée par $$u(x) = \\exp(x) - 1 - (\\exp(1) - 1)x.$$" ] }, { "cell_type": "code", "execution_count": null, "id": "1531c439", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "8bf4c31e", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "e5eeb6a5", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "encoding": "# -*- coding: utf-8 -*-" }, "kernelspec": { "display_name": "Julia 1.8.5", "language": "julia", "name": "julia-1.8" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.5" } }, "nbformat": 4, "nbformat_minor": 5 }