{
"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
}