{ "cells": [ { "cell_type": "markdown", "id": "657ee43e-abd8-4efe-b5e0-6a632114529a", "metadata": { "tags": [] }, "source": [ "# Cours ENPC - Pratique du calcul scientifique" ] }, { "cell_type": "markdown", "id": "a30eaeba", "metadata": {}, "source": [ "## Intégration numérique" ] }, { "cell_type": "markdown", "id": "532fe401", "metadata": {}, "source": [ "### Exercice théorique sur les polynômes de Legendre\n", "\n", "Pour $n∈\\mathbb{N}$, on définit le polynôme de Legendre $L_n$ par la formule de Rodrigues :\n", "$$L_n(x)=\\frac{1}{n!2^n}\\frac{\\mathrm{d}^n}{\\mathrm{d} x^n}\\left((x^2-1)^n\\right)$$\n", "\n", "On définit également le produit scalaire suivant dans l'espace des fonctions continues de $[-1,1]$ dans $\\mathbb{R}$ :\n", "\n", "\n", "$$\n", "\\tag{1}\n", "\\begin{array}{rccl}\n", "⟨•,•⟩ :&\n", "\\mathcal{C}([-1,1])×\\mathcal{C}([-1,1])&→&\\mathbb{R} \\\\\n", "&(f,g)&↦&⟨f,g⟩=\\int_{x=-1}^1 f(x) g(x) \\mathrm{d} x\n", "\\end{array}\n", "$$\n", "\n", "Démontrer les propriétés suivantes de la famille des polynômes $L_n$ :\n", "\n", "1. Le degré de $L_n$ est $n$\n", "\n", "1. $L_n$ a la même parité que $n$ : $L_n(-x)=(-1)^n L_n(x)$\n", "\n", "1. Le coefficient dominant de $L_n$ est $\\frac{(2n)!}{n!^2 2^n}=\\frac{{2n \\choose n}}{2^n}$\n", "\n", "1. $L_n$ peut se réécrire sous la forme suivante\n", " $$L_n(x)=\\frac{1}{2^n} \\sum_{k=0}^n {n\\choose k}^2 (x+1)^{n-k} (x-1)^k$$\n", " avec ${n\\choose k}=\\frac{n!}{k!(n-k)!}$\n", "\n", "1. $L_n(1)=1$ et $L_n(-1)=(-1)^n$\n", "\n", "1. $(L_k)_{0≤k≤n}$ forme une base de $\\mathbb{R}_n[X]$\n", "\n", "1. $L_n$ est orthogonal à tout polynôme de $\\mathbb{R}_{n-1}[X]$ au sens du produit scalaire (1) \n", " *Indication : utiliser des intégrations par parties.*\n", "\n", "1. $L_n$ est scindé à racines simples dans $]-1,1[$ \n", " *Indication : montrer par l'absurde en supposant que $L_n$ change $k$ fois de signe sur $]-1,1[$ avec $k\n", " $$\n", " \\tag{2}\n", " \\begin{align*}\n", " &L_0=1, \\quad L_1=X\\\\\n", " &(n+1)L_{n+1}=(2n+1)X L_n - n L_{n-1}\\;(n≥1)\n", " \\end{align*}\n", " $$\n", " *Indication : considérer une décomposition adéquate de $X L_n$, exploiter les propriétés d'orthogonalité, de parité et de termes dominants.*\n", "\n", "1. La base orthonormée au sens de (1) formée des polynômes $φ_n=\\frac{L_n}{\\lVert L_n \\rVert}=\\sqrt{\\frac{2n+1}{2}}L_n$ satisfait la rélation de récurrence suivante\n", " $$\n", " \\begin{align*}\n", " &φ_0=\\frac{1}{\\sqrt{2}}, \\quad φ_1=\\sqrt{\\frac{3}{2}}X\\\\\n", " &\\frac{n+1}{\\sqrt{2n+3}}\\,φ_{n+1}=\\sqrt{2n+1}\\,X \\,φ_n - \\frac{n}{\\sqrt{2n-1}}\\, φ_{n-1}\\;(n≥1)\n", " \\end{align*}\n", " $$\n", "\n", "1. En déduire que les polynômes $φ_n$ satisfont l'équation matricielle suivante pour $N∈\\mathbb{N}$ fixé\n", " \n", " $$\n", " \\tag{3}\n", " \\begin{pmatrix}\n", " x φ_0(x) \\\\\n", " x φ_1(x) \\\\\n", " x φ_2(x) \\\\\n", " \\vdots \\\\\n", " x φ_{N-1}(x) \\\\\n", " x φ_N(x)\n", " \\end{pmatrix}\n", " =\n", " \\underbrace{\n", " \\begin{pmatrix}\n", " 0 & α_1 \\\\\n", " α_1 & 0 & α_2 \\\\\n", " & α_2 & 0 & α_3 \\\\\n", " & & \\ddots & \\ddots & \\ddots \\\\\n", " & & & α_{N-1} & 0 & α_N \\\\\n", " & & & & α_N & 0\n", " \\end{pmatrix}\n", " }_{\\mathsf A}\n", " \\underbrace{\n", " \\begin{pmatrix}\n", " φ_0(x) \\\\\n", " φ_1(x) \\\\\n", " φ_2(x) \\\\\n", " \\vdots \\\\\n", " φ_{N-1}(x) \\\\\n", " φ_N(x)\n", " \\end{pmatrix}\n", " }_{\\mathbf{\\boldsymbol{Φ}}(x)}\n", " +\n", " \\begin{pmatrix}\n", " 0 \\\\\n", " 0 \\\\\n", " 0 \\\\\n", " \\vdots \\\\\n", " 0 \\\\\n", " α_{N+1} φ_{N+1}(x)\n", " \\end{pmatrix}\n", " $$\n", " avec $α_n=\\frac{n}{\\sqrt{(2n+1)(2n-1)}}=\\frac{n}{\\sqrt{4n^2-1}}$\n", "\n", "1. En exploitant les racines $(r_n)_{0≤n≤N}$ de $L_{N+1}$ dans (3) déterminer les valeurs propres et vecteurs propres de $\\mathsf A$. Préciser pourquoi les vecteurs propres forment une base orthogonale de $\\mathbb{R}^{N+1}$ (au sens du produit scalaire canonique de $\\mathbb{R}^{N+1}$).\n", "\n", "1. On forme la matrice suivante $\\mathsf{P}$ dans laquelle on place les vecteurs $Φ(r_n)_{0≤n≤N}$ selon les lignes\n", " $$\n", " \\mathsf{P} =\n", " \\begin{pmatrix}\n", " φ_0(r_0) & φ_1(r_0) & \\dots & φ_N(r_0) \\\\\n", " φ_0(r_1) & φ_1(r_1) & \\dots & φ_N(r_1) \\\\\n", " φ_0(r_2) & φ_1(r_2) & \\dots & φ_N(r_2) \\\\\n", " \\vdots & \\vdots & \\dots & \\vdots \\\\\n", " φ_0(r_N) & φ_1(r_N) & \\dots & φ_N(r_N)\n", " \\end{pmatrix}\n", " $$\n", " Montrer que $\\mathsf{P}\\mathsf{P}^T$ est une matrice diagonale que l'on note $\\mathsf{D}=\\mathsf{P}\\mathsf{P}^T$ avec\n", " $$\n", " d_{ii}=\\lVert Φ(r_i)\\rVert^2=\\sum_{n=0}^N φ_n(r_i)^2\n", " $$ \n", " **NB : il est très important de noter pour l'implémentation numérique ultérieure que les vecteurs $Φ(r_i)$ ne sont a priori pas unitaires au sens de la norme canonique de $\\mathbb{R}^{N+1}$ mais il est possible de retrouver $Φ(r_i)$ à partir d'un vecteur propre associé à $r_i$ par le fait que sa première composante vaut $φ_0(r_i)=\\frac{1}{\\sqrt{2}}$.**\n", "\n", "1. Déduire que les vecteurs définis par les **colonnes** de $\\mathsf{P}$, autrement dit les vecteurs\n", " $$\n", " \\mathbf{\\boldsymbol{Ψ}}_n =\n", " \\begin{pmatrix}\n", " φ_n(r_0) \\\\\n", " φ_n(r_1) \\\\\n", " φ_n(r_2) \\\\\n", " \\vdots \\\\\n", " φ_n(r_N)\n", " \\end{pmatrix}\n", " \\quad 0≤n≤N\n", " $$\n", " forment une base orthonormale au sens du produit scalaire de $\\mathbb{R}^{N+1}$ dont la matrice dans la base canonique est $\\mathsf{D}^{-1}$. \n", " *Indication : la question revient à montrer que $\\mathsf{P}^T \\mathsf{D}^{-1}\\mathsf{P}=\\mathsf{I}$ où $\\mathsf{P}$ joue le rôle d'une matrice de passage.* \n", "\n", "1. On introduit le produit scalaire sur l'espace des polynômes de degré inférieur ou égal à $N$\n", " $$\n", " \\begin{array}{rccl}\n", " ⟨•,•⟩_{N+1} :&\n", " \\mathbb{R}_N[X]×\\mathbb{R}_N[X]&→&\\mathbb{R} \\\\\n", " &(P,Q)&↦&⟨P,Q⟩_{N+1}=\\sum_{i=0}^N w_i P(r_i) Q(r_i)\n", " \\end{array}\n", " $$\n", " où l'on a noté\n", " \n", " $$\n", " \\tag{4}\n", " w_i=\\frac{1}{d_{ii}}=\\frac{1}{\\lVert Φ(r_i)\\rVert^2}=\\frac{1}{\\sum_{n=0}^N φ_n(r_i)^2}\n", " $$\n", " Montrer que $⟨•,•⟩_{N+1}$ définit bien un produit scalaire sur $\\mathbb{R}_N[X]$ et que l'isomorphisme suivant entre espaces euclidiens est une isométrie :\n", " $$\n", " \\begin{array}{ccl}\n", " (\\mathbb{R}_N[X],⟨•,•⟩)&→&(\\mathbb{R}_N[X],⟨•,•⟩_{N+1}) \\\\\n", " P&↦&P\n", " \\end{array}\n", " $$ \n", " *Indication : étudier les valeurs prises par les deux produits scalaires sur la base des polynômes $(φ_n)_{0≤n≤N}$.*\n", "\n", "1. Déduire que pour tout polynôme $R∈\\mathbb{R}_N[X]$ on a\n", " $$\n", " \\int_{x=-1}^1 R(x) \\mathrm{d} x = \\sum_{i=0}^N w_i R(r_i)\n", " $$ \n", " *Indication : considérer les polynômes $R$ et $1$.*\n", "\n", "1. Montrer que le résultat précédent se généralise à tout polynôme de degré inférieur ou égal à $2N+1$\n", " \n", " $$\n", " \\tag{5}\n", " ∀ P∈ \\mathbb{R}_{2N+1}[X],\\;\n", " \\int_{x=-1}^1 P(x) \\mathrm{d} x = \\sum_{i=0}^N w_i P(r_i)\n", " $$ \n", " *Indication : utiliser la division euclidienne de $P$ par $L_{N+1}$.*\n", "\n", "1. Montrer qu'une méthode alternative à (4) pour identifier les poids repose sur l'expression suivante\n", " \n", " $$\n", " \\tag{6}\n", " ∀ i ∈ [0,1,\\dotsc,N], \\;\n", " w_i=\\frac{2}{(1-r_i^2) \\left(L'_{N+1}(r_i)\\right)^2}\n", " $$\n", " *Indication : pour $i$ donné, commencer par poser $L_{N+1}=(X-r_i)Q$ où $Q$ est un polynôme de degré $N$ et écrire $⟨L'_{N+1},Q⟩$ de plusieurs façons.*" ] }, { "cell_type": "markdown", "id": "2ff5e83c", "metadata": {}, "source": [ "### Exercice numérique de construction des points et poids de Gauss\n", "\n", "1. Appliquer à la main les algorithmes exposés dans l'exercice théorique dans le cas $N=1$ soit $2$ points et poids de Gauss-Legendre\n", "\n", " - méthode principale s'appuyant sur la diagonalisation de la matrice $\\mathsf{A}$ en (3),\n", "\n", " - méthode alternative s'appuyant sur la construction explicite du polynôme de Legendre d'ordre $N+1=2$ par la récurrence (2) puis le calcul des racines de $L_2$ et enfin des poids par (6).\n", "\n", " *Réponse : on doit trouver dans les deux cas $r_1=-r_0=\\frac{1}{\\sqrt{3}}$ et $w_0=w_1=1$*\n", "\n", "1. Implémenter une fonction en `Julia` qui à tout entier `N` en entrée renvoie `N+1` points et poids de Gauss-Legendre selon la méthode principale exposée dans l'exercice théorique. \n", "\n", " *Indications :*\n", " \n", " - La matrice $\\mathsf{A}$ en (3) pourra être construite à l'aide de l'opérateur [`LinearAlgebra.SymTridiagonal`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.SymTridiagonal).\n", "\n", " - Les valeurs propres et vecteurs propres d'une matrice peuvent être calculées par [`LinearAlgebra.eigen`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.eigen). \n", " On notera d'une part que les valeurs propres sont triées par ordre croissant et sont renvoyées sous forme du vecteur `E.values` si `E` désigne le résultat de l'opérateur `LinearAlgebra.eigen` et d'autre part que les vecteurs propres sont normalisés et sont les colonnes de la matrice renvoyée par `E.vectors` si bien que le $i^\\textrm{ème}$ vecteur propre est donné par `E.vectors[:,i]` (attention la numérotation en `Julia` démarre à **1** comme en `Fortran` et non à **0** comme en `Python` ou en `C/C++`).\n", "\n", " - Attention à bien tenir compte du *Nota Bene* exprimé dans l'exercice théorique s'agissant de la normalisation des vecteurs $Φ(r_i)$ qu'il n'est en réalité pas utile de calculer en entier pour extraire les poids de Gauss-Legendre.\n", "\n", "1. Implémenter la méthode alternative en calculant $L_{N+1}$ à partir de la formule de récurrence (2) puis ses racines et les poids par (6).\n", "\n", " *Indications :*\n", "\n", " - On pourra s'appuyer sur la bibliothèque [`Polynomials.jl`](https://github.com/JuliaMath/Polynomials.jl) pour les calculs algébriques sur les polynômes.\n", "\n", " - Les racines sont obtenues par la fonction [`Polynomials.roots`](https://juliamath.github.io/Polynomials.jl/stable/reference/#Polynomials.roots). Attention à bien les réordonner (avec `sort`) dans l'objectif de les comparer à ceux obtenus par la première méthode.\n", "\n", " - La dérivée d'un polynôme peut être calculée par [`Polynomials.derivative`](https://juliamath.github.io/Polynomials.jl/stable/reference/#Polynomials.derivative). Alternativement, il est également possible de construire la dérivée par exploitation d'une relation de récurrence couplée à celle des polynômes $L_n$ (2).\n", "\n", "1. Comparer les résultats des deux méthodes sur plusieurs valeurs de `N` (par exemple `N ∈ 0:10`). (Code fourni.)\n", "\n", "1. Vérifier que les polynômes $X^d$ ($0≤d≤2N+1$) sont intégrés exactement. (Code fourni.)" ] }, { "cell_type": "code", "execution_count": 1, "id": "36dc5f39", "metadata": {}, "outputs": [], "source": [ "function gauss_legendre(N)\n", " # Votre code pour la question 2.\n", " return nodes, weights\n", "end" ] }, { "cell_type": "code", "execution_count": 2, "id": "80959572", "metadata": {}, "outputs": [], "source": [ "# Il peut être utile de définir une fonction auxiliaire ici\n", "\n", "function gauss_legendre_bis(N)\n", " # Votre code pour la question 3.\n", " return nodes, weights\n", "end" ] }, { "cell_type": "code", "execution_count": 5, "id": "ac96e1d2", "metadata": {}, "outputs": [], "source": [ "GaussLegendre = Tuple{Vector{Float64}, Vector{Float64}}[]\n", "for N ∈ 0:10\n", " x, w = gauss_legendre(N) ; push!(GaussLegendre,(x,w))\n", " x̃, w̃ = gauss_legendre_bis(N)\n", " println(\"$(N+1) point$(N>0 ? \"s\" : \"\") → $(x ≈ x̃ && w ≈ w̃ ? \"✅\" : \"❌\")\")\n", "end\n", "x, w = gauss_legendre(1)\n", "println(\"$(x) ≈ [-1/√3,1/√3] → $(x ≈ [-1/√3,1/√3] ? \"✅\" : \"❌\")\")\n", "println(\"$(w) ≈ [1,1] → $(w ≈ ones(2) ? \"✅\" : \"❌\")\")" ] }, { "cell_type": "code", "execution_count": 4, "id": "78591ba0", "metadata": {}, "outputs": [], "source": [ "for N ∈ 0:length(GaussLegendre)-1\n", " x, w = GaussLegendre[N+1]\n", " res=[]\n", " for d ∈ 0:2:2N+1\n", " xᵈ(x) = x^d\n", " push!(res,\"∫₋₁¹x^$(d)dx → $(sum(w .* xᵈ.(x)) ≈ 2/(d+1) ? \"✅\" : \"❌\") ; \")\n", " end\n", " println(res...)\n", "end" ] }, { "cell_type": "markdown", "id": "e2b18712", "metadata": {}, "source": [ "### Application numérique de l'intégration de Gauss-Legendre\n", "\n", "1. Construire une fonction qui prend en entrée une fonction `u`, des bornes `a` et `b` et un nombre `N` (définissant `N+1` points de Gauss selon la convention prise plus haut) renvoyant l'intégration numérique par Gauss-Legendre de la fonction.\n", "\n", "1. Construire une fonction prenant en entrée une fonction `u`, des bornes `a` et `b` et un pas `h` et renvoyant l'intégration de `u` entre `a` et `b` selon la méthode trapézoïdale composite de pas `h`. On commencera par transformer le pas `h` en un nombre d'intervalle permettant d'approcher au mieux `h` par la formule suivante `N = round(Int,(b-a)/h)`.\n", "\n", "1. Faire de même avec la méthode de Simpson.\n", "\n", "1. Intégrer la fonction de Runge $u(x)=1/(1+25x^2)$ sur $[-1,1]$ analytiquement puis tracer en fonction du nombre d'intervalles d'intégration les erreurs relatives correspondant aux différentes méthodes d'intégration numérique. On pourra prendre une progression logarithmique du nombre d'intervalles entre $2^3$ et $2^8$ (tableau construit par `2 .^(3:8)`) et tracer les erreurs en échelles log en x et y (options `xaxis=:log10` et `yaxis=:log10` dans un `plot`).\n", "\n", "1. Appliquer la méthode de Romberg à 1 et 2 itérations à la formule d'intégration trapézoïdale, i.e. implémenter $J_1(h)=\\frac{4J(h)-J(2h)}{3}$ et $J_2(h)=\\frac{16J_1(h)-J_1(2h)}{15}$ si $J(h)$ désigne le résultat de l'intégration trapézoïdale et ajouter les graphes des erreurs obtenues aux précédents.\n", "\n", " *Indications :*\n", "\n", " - si `U` et `V` sont des vecteurs colonnes, `[U ; V]` renvoie le vecteur colonne issu de la concaténation de `U` et `V`\n", "\n", " - `sum([U ; V])` renvoie donc la somme de tous les éléments de `U` et `V` donc équivalent à `sum(U)+sum(V)`\n", "\n", " - `U[1]` ou `U[begin]` renvoie le premier terme du vecteur `U`\n", "\n", " - `U[end]` renvoie le dernier terme du vecteur `U`\n", " \n", " - `U[end-1]` renvoie l'avant-dernier terme du vecteur `U` etc.\n", "\n", " - `U[2:2:end-1]` renvoie le vecteur extrait de `U` ne comportant que les indices pairs entre le deuxième et l'avant-dernier\n", "\n", " - `U[3:2:end-2]` renvoie le vecteur extrait de `U` ne comportant que les indices impairs entre le troisième et l'antépénultième" ] }, { "cell_type": "code", "execution_count": 6, "id": "c2cd9b59", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 7, "id": "dac31df6", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 8, "id": "791dd77a", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 9, "id": "7d079a4f", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 10, "id": "6ee22503", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "citation-manager": { "items": {} }, "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 }