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