{ "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": [ "## Interpolation et approximation\n", "\n", "Les graphiques peuvent être tracés grâce à la bibliothèque [`Plots.jl`](https://github.com/JuliaPlots/Plots.jl)\n", "\n", "- [dépôt GitHub](https://github.com/JuliaPlots/Plots.jl)\n", "\n", "- [documentation](https://docs.juliaplots.org/stable/)\n", "\n", "- [tutoriel](https://docs.juliaplots.org/stable/tutorial/)\n", "\n", "- Il est possible de définir des options de tracé par défaut (qui peuvent être redéfinies ponctuellement au besoin) par\n", "\n", " ```julia\n", " Plots.default(fontfamily=\"Computer Modern\",\n", " titlefontsize=20,\n", " xlabelfontsize=20,\n", " ylabelfontsize=20,\n", " legendfontsize=16,\n", " xtickfontsize=16,\n", " ytickfontsize=16,\n", " linewidth=2,\n", " markersize=7,\n", " framestyle=:box,\n", " label=nothing,\n", " grid=true)\n", " ```" ] }, { "cell_type": "code", "execution_count": null, "id": "16be1d44", "metadata": {}, "outputs": [], "source": [ "# Quelques bibliothèques utiles. Vous pouvez modifier les paramètres à votre guise et ajouter des bibliothèques\n", "# si besoin (sauf s'il explicitement mentionné dans un exercice de ne pas utiliser de bibliothèque externe)\n", "\n", "using Polynomials, LinearAlgebra, 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", " label=nothing,\n", " grid=true)\n", "\n", "# For information on any of these, type\n", "# plotattr(\"framestyle\")" ] }, { "cell_type": "markdown", "id": "7d4cb731", "metadata": {}, "source": [ "### Exercice sur l'interpolation de Lagrange en arthmétique inexacte\n", "\n", "On étudie la performance de l’interpolation de Lagrange en arithmétique inexacte.\n", "\n", "1. Créer une fonction `get_lagrange_interpolant` en `Julia` qui accepte comme arguments un vecteur `x` de nœuds équidistants et un vecteur `u` de valeurs d'une fonction $u$, que l'on veut interpoler en ces nœuds, et qui construit et renvoie l'interpolant de Lagrange, $\\hat{u}$ comme la fonction `interpolant` en `Julia`. Assurez-vous que dans l'intérieur de votre fonction l'arithmétique est faite dans la même précision que l'arguement qu'elle reçoit. Indice: Regardez la documentation des fonctions `one` et `zero`.\n", "\n", "1. Tester votre fonction `get_lagrange_interpolant` sur la fonction $f:[-1, 1] \\to \\mathbb{R}$, $f(x) = 1 $ pour tout $x \\in [-1, 1]$ pour plusieurs $n$ et avec des nœuds et valeurs de type `Float64`, `Float32`, et `Float16`.\n", "\n", " Commenter ce qui se passe lorsque l'on utilise plus de nœuds et les types de moindre précision. Que peut expliquer ce comportement ? Quel est l'interpolant à $n$ nœuds équidistants de $f$ en arithmétique exacte ?\n", "\n", "1. *Facultatif :* examiner le cas d'une précision à 128 bits (type `Float128` disponible avec la bibliothèque [`Quadmath.jl`](https://github.com/JuliaMath/Quadmath.jl)" ] }, { "cell_type": "code", "execution_count": null, "id": "9cd79533", "metadata": {}, "outputs": [], "source": [ "function get_lagrange_interpolant(x, u)\n", " \"\"\"\n", " Définir une fonction qui renvoie une fonction `interpolant` qui \n", " est une interpolation de Lagrange de noeuds `x` et valeurs `u`\n", " \"\"\"\n", " function interpolant(y)\n", " \"\"\"\n", " Interpolant û comme il est défini dans le poly\n", " \"\"\"\n", " # Ajouter votre code ici : la fonction doit retourner la valeur de û au point y\n", " s = zero(eltype(u))\n", "\n", " n = length(x)\n", " for i=1:n\n", " s += prod((y-x[j])/(x[i]-x[j]) for j=1:n if j!=i)\n", " end\n", "\n", " return s\n", " end\n", " \n", " return interpolant\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "0fe9f855", "metadata": {}, "outputs": [], "source": [ "# Pour tester votre fonction, vous pouvez utiliser la fonction suivante en tournant les cellules qui suivent\n", "\n", "using Plots\n", "function lagrange_tester(n, m, num_type, fun)\n", " \"\"\"\n", " Tester `get_lagrange_interpolant` sur la fonction `fun`\n", " pour tout x avec des noeud d'une précision donné par `num_type`.\n", " Appliquer l'interpolant renvoyé par `get_lagrange_interpolant`\n", " à `m` point équidistant\n", " \n", " Arguments\n", " ---------\n", " n : nombre des noeuds\n", " \n", " m : nombre des points à évaluer l'interpolant\n", " \n", " num_type : type de noeuds et valeurs d'interpolation, ex. Float64\n", " \n", " fun: fonction à interpoler\n", " \n", " Retour\n", " ------\n", " \n", " p1 : plot de l'interpolation evalué à n+1 points equidistants et la vraie fonction\n", " \"\"\"\n", " # créer n noeuds equidistants et n valeurs de type num_type\n", " x = LinRange{num_type}(-1, 1, n)\n", " # idem avec m noeuds pour le tracé\n", " x_plot = LinRange{num_type}(-1, 1, m)\n", " û = get_lagrange_interpolant(x, fun.(x))\n", " \n", " # Dessiner l'interpolation avec m points equidistant\n", " x_plot = LinRange{num_type}(-1, 1, m)\n", " p1 = plot(x_plot, û.(x_plot), label = \"Interpolation\")\n", " \n", " # Dessiner la vraie fonction \n", " plot!(x_plot, fun.(x_plot), label = \"Vraie Fonction\")\n", " \n", " print(\"Temps d'evaluation d'interpolant: \")\n", " @time û.(x_plot)\n", " \n", " print(\"Erreur d'interpolation: $(maximum(x -> abs(û(x)-fun(x)), x_plot))\") \n", " \n", " return p1\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "fb152788-e536-4f93-a37e-f47113168513", "metadata": {}, "outputs": [], "source": [ "m = 200\n", "fun(x) = 1" ] }, { "cell_type": "code", "execution_count": null, "id": "6699a5d5", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(31, m, Float64, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "7a923552", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(41, m, Float64, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "3117ce70", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(51, m, Float64, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "f958d7c1", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(61, m, Float64, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "0b513249", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(11, m, Float32, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "372c22e9", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(21, m, Float32, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "d1ae2483", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(11, m, Float16, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "696f36e2", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(21, m, Float16, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "45c6e875", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(101, m, Float64, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "51ad0cfc", "metadata": {}, "outputs": [], "source": [ "using Quadmath\n", "\n", "lagrange_tester(101, m, Float128, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "d35afe33", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(121, m, Float128, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "7094932b", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(131, m, Float128, fun)" ] }, { "cell_type": "code", "execution_count": null, "id": "9201e4b3", "metadata": {}, "outputs": [], "source": [ "lagrange_tester(141, mj, Float128, fun)" ] }, { "cell_type": "markdown", "id": "5c31501e", "metadata": {}, "source": [ "### Exercice sur les noeuds d'interpolation\n", "\n", "Écrire un code `Julia` pour interpoler la fonction suivante à l'aide d'un polynôme de degré 20 sur l'intervalle $[-1, 1]$.\n", "$$\n", " f(x) = \\tanh\\left(\\frac{x+1/2}{\\varepsilon}\\right) + \\tanh\\left(\\frac{x}{\\varepsilon}\\right) + \\tanh\\left(\\frac{x-1/2}{\\varepsilon}\\right),\n", " \\qquad \\varepsilon = .01\n", "$$\n", "Utiliser des noeuds équidistants puis des noeuds de Tchebychev et comparer les deux approches en termes de précision.\n", "Tracer la fonction $f$ ainsi que les polynômes d'interpolation.\n", "\n", "*Indications :*\n", "\n", "- Pour limiter les erreurs d'arrondi numérique, il est préférable que la fonction renvoie un type `BigFloat`, autrement dit\n", "\n", " ```julia\n", " f(x) = BigFloat(tanh((x+1/2)/ε) + tanh(x/ε) + tanh((x-1/2)/ε))\n", " ```\n", "\n", "- Pour calculer rapidement les noeuds de Tchebychev, on peut exploiter la macro `@.` (comme toujours, il est conseillé de se référer à la documentation d'une commande en tapant `?` puis la commande dans la console). Cette commande évite d'utiliser des `.` après chaque fonction ou avant chaque opérateur.\n", "\n", " ```julia\n", " x = @. -cos(π*((0:n-1)+1/2)/n)\n", " ```\n", "\n", "- Le calcul du polynôme interpolateur peut être obtenu par la fonction [`fit`](https://juliamath.github.io/Polynomials.jl/stable/#Fitting-arbitrary-data) de la bibliothèque [`Polynomials.jl`](https://github.com/JuliaMath/Polynomials.jl).\n", "\n", "- Il peut être utile pour comparer les deux interpolations de limiter les valeurs minimale et maximale sur l'axe `y` à l'aide de l'option `ylims = (ymin,ymax)` dans une fonction de tracé `plot`, `scatter`, ou leurs équivalents terminant par `!`. On rappelle que, par convention en `Julia` (et non par obligation), une fonction dont le nom se termine par `!` modifie ses arguments. Dans le cas d'un graphe, la première commande initiant le graphe ne doit pas comporter de `!` (`plot`, `scatter`, ...) tandis que les suivantes incrémentant le même graphe se terminent par `!` (`plot!`, `scatter!`, ...). Toute omission du `!` est considéré comme un *redémarrage à zéro* du tracé.\n", "\n", "- Pour calculer la norme infinie d'une fonction afin d'évaluer la précision de l'interpolation, on pourra exploiter la fonction [`norm(...,Inf)`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.norm) de la bibliothèque `LinearAlgebra` avec une échantillonnage suffisamment fin des valeurs de la fonction ou exploiter la fonction `maximum`. Noter que la conversion d'un nombre `y` de type `BigFloat` en `Float64` se fait par `convert(Float64,y)` ou plus simplement ici `Float64(y)`.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "6ee22503", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "e62b3987", "metadata": {}, "source": [ "### Exercice sur les noeuds d'interpolation : animation de graphes\n", "\n", "En exploitant la macro [`@animate`](https://docs.juliaplots.org/latest/animations/) de la bibliothèque `Plots.jl`, créer une animation permettant de voir l'évolution superposée des interpolations avec noeuds équidistants et de Tchebychev de la fonction de Runge $u(x)=\\frac{1}{1+25x^2}$ sur l'intervalle $[-1,1]$. On fera varier le nombre de noeuds par exemple de 2 à 50 et on pourra tenir compte des indications de l'exercice précédent pour limiter par exemple l'extension en `y`." ] }, { "cell_type": "code", "execution_count": 1, "id": "317cfc7b", "metadata": {}, "outputs": [], "source": [] }, { "attachments": {}, "cell_type": "markdown", "id": "c0f716ce", "metadata": {}, "source": [ "### Exercice sur l'équation de Laplace\n", "\n", "On se propose d'implémenter une méthode numérique pour résoudre approximativement l'équation de Laplace avec conditions au bord homogène de Dirichlet:\n", "\n", "$$ u\\in C^2([0,1]),\\quad\\left\\{\\begin{aligned} -u''(x) & = f(x) & \\forall\\, x\\in(0,1),\\\\ u(0) & = u(1) = 0. \\end{aligned}\\right.$$\n", "Pour cela, on approxime $f$ avec un polynome interpolateur $\\hat f$, puis on résoud exactement l'équation de Laplace associée. Implémenter cette méthode.\n", "On pourra tester le cas où\n", "$$f(x) = \\exp(\\sin(2\\pi x))(\\sin(2\\pi x)-\\cos(2\\pi x)^2),$$\n", "auquel cas la solution analytique est donnée par\n", "$$ u(x)=(2\\pi)^{-2}\\exp((\\sin(2\\pi x))-1).$$\n", "\n", "*Indications :*\n", "- On pourra utiliser la fonction `fit` de la bibliothèque `Polynomials.jl` pour obtenir le polynôme interpolateur:\n", " ```julia\n", " p = fit(x,y)\n", " ```\n", "où `x` sont les noeuds d'interpolation et `y` sont les valeurs de la fonction à interpoler.\n", "- Pour calculer la solution analytique avec membre de droite polynomial, on pourra remarquer que toutes les solutions sont des polynômes, et que, sans condition au bord, la solution est unique modulo $\\mathbf{P}_1$.\n", "- On pourra utiliser la fonction `integrate` de la bibliothèque `Polynomials.jl` qui calcule une primitive d'un polynôme:\n", " ```julia\n", " P = integrate(p)\n", " ```\n", "- Utiliser le format `BigFloat` pour limiter les erreurs d'arrondi" ] }, { "cell_type": "markdown", "id": "a67c5aa2", "metadata": {}, "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 }