{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Approcher une dérivée\n", "\n", "Marc Lorenzi\n", "\n", "25 avril 2019" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import math\n", "from sympy import *\n", "init_printing()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = (12, 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous allons dans ce notebook \n", "\n", "- Fabriquer des formules permettant d'approcher la dérivée d'une fonction $f$ en un point $x$.\n", "- Estimer l'erreur commise dans ces formules.\n", "- Minimiser l'erreur commise.\n", "\n", "Hormis dans la dernière partie, nous nous intéresserons uniquement à des dérivées d'ordre 1. Nous verrons tout à la fin que la généralisation aux dérivées d'ordre supérieur ne présente aucune difficulté.\n", "\n", "Au travail !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Taux d'accroissement" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 L'approximation \"naïve\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tout le monde le sait, pour approcher la dérivée de la fonction $f$ au point $x$ il n'y a qu'à prendre la quantité\n", "\n", "$$\\frac{f(x+h)-f(x)}{h}$$\n", "\n", "où $h$ est réel \"petit\" ? Ce que tout le monde ne sait peut-être pas, c'est que cette approximation n'est pas très bonne. Pourquoi ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici tout d'abord la fonction `approx_deriv`. Elle prend en paramètres une fonction $f$, un réel $x$ et un réel $h\\ne 0$ et renvoie le quotient dont nous venons de parler." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def approx_deriv(f, x, h):\n", " return (f(x + h) - f(x)) / h" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Remarque__ : Dans tout le notebook je considérerai $f(x)=e^x$, en $x=1$. Toutes les dérivées exactes de $f$ en 1 valent évidemment $e$. Il vous est fortement conseillé de tester sur d'autres fonctions et d'autres points !\n", "\n", "Ci-dessous, le graphe de l'erreur commise en fonction de $h$. La fonction `tracer_approx` prend en paramètres\n", "\n", "- une fonction $f$\n", "- Un réel $x$\n", "- Une fonction $f_1$ qui est censée être la fonction dérivée de $f$\n", "- Une fonction `fonc_approx` qui est du genre de notre fonction `approx_deriv`\n", "\n", "La fonction trace l'erreur commise en fonction de $h$, en coordonnées logarithmiques." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def tracer_approx(f, x, f1, fonc_approx):\n", " hs = [10. ** (-k / 50) for k in range(500)]\n", " ds = [abs(fonc_approx(f, x, h) - f1(x)) for h in hs]\n", " plt.loglog(hs, ds, 'k')\n", " plt.grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tracer_approx(math.exp, 1, math.exp, approx_deriv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bigre ! Que voit-on ? En lisant le dessin de droite à gauche, on voit que lorsque $h$ diminue, l'erreur commise diminue. Ceci n'est pas une surprise. Sauf que lorsque $h$ devient plus petit que, environ, $10^{-8}$, l'erreur se met à fluctuer de façon \"imprévisible\" et, globalement, à ré-augmenter. Que se passe-t-il ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 Taylor-Lagrange" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous allons utiliser dans tout ce document une version du théorème de Taylor-Lagrange que voici.\n", "\n", "__Théorème__ : Soit $x\\in\\mathbb R$. Soit $f$ une fonction définie dans un voisinage $V=[x-\\alpha,x+\\alpha]$ ($\\alpha>0$) de $x$. Soit $n\\in\\mathbb N$. On suppose que\n", "\n", "- $f$ est de classe $\\mathcal C^n$ sur $V$.\n", "- $f$ est $n+1$ fois dérivable sur $V$.\n", "\n", "Soit $h>0$ assez petit pour que $x+h\\in V$. Alors il existe $c\\in]x,x+h[$ tel que\n", "\n", "$$f(x+h)=\\sum_{k=0}^n \\frac{x^k}{k!}f^{(k)}(x)+\\frac{h^{n+1}}{(n+1)!}f^{(n+1)}(c)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $f$ une fonction de classe $\\mathcal C^2$ dans un voisinage du réel $x$. Soit $h>0$ assez petit pour que $x+h$ soit dans ce voisinage. La formule de Taylor-Lagrange appliquée avec $n=1$ nous dit que\n", "\n", "$$f(x+h)=f(x)+hf'(x)+\\frac {h^2}2f''(c)$$\n", "\n", "où $x0$. L'erreur commise est donc inférieure à \n", "\n", "$$\\varphi(h)=\\frac 1 2 M_2 h + 2\\frac{\\varepsilon} h$$\n", "\n", "On a\n", "\n", "$$\\varphi'(h)=\\frac 1 2 M_2 - 2\\frac{\\varepsilon} {h^2}$$\n", "\n", "et on constate que $\\varphi$ passe par un minimum pour\n", "\n", "$$h=\\sqrt{\\frac{4\\varepsilon}{M_2}}$$\n", "\n", "La valeur minimale de $\\varphi$ est alors, par un petit calcul,\n", "\n", "$$\\frac 3 2 \\sqrt{M_2\\varepsilon}$$\n", "\n", "__Exercice__ : Faites le petit calcul.\n", "\n", "__Conséquence__ : Sauf coup de chance (erreurs d'arrondi magiquement petites) on ne peut pas obtenir la valeur de $f'(x)$ avec une précision meilleure que quelque chose de l'ordre de $\\sqrt\\varepsilon$.\n", "\n", "__Remarque__ : La norme IEEE754 décrit la façon dont sont codés les différents formats de réels. Le langage Python utilise le format __double__, dans lequel la \"précision\" des nombres est de 52 bits (je n'entre pas dans les détails ici). Une bonne valeur pour $\\varepsilon$ est donc sans doute\n", "\n", "$$\\varepsilon = 2^{-52}>10^{-16}$$\n", "\n", "Comme $\\sqrt\\varepsilon>10^{-8}$, nous ne pouvons donc pas compter, sauf miracle, sur plus de 8 chiffres exacts après la virgule ! \n", "\n", "__Remarque__ : Pour ceux d'entre-vous qui voudraient en savoir plus, le module `sys` contient un objet `float_info` qui contient tous les renseignement que l'on peut désirer sur la structure des `float` en Python." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici donc le $\\epsilon$ que nous prendrons dans le notebook :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "2 ** (-52)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et voici la \"précision\" des flottants en Python :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "52 * math.log(2, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Qu'obtenons-nous comme valeurs numériques du $h$ optimal et comme majorant de l'erreur minimale ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def h_optimal(eps, M):\n", " return math.sqrt(4 * eps / M)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h0 = h_optimal(2 ** (-52), math.exp(1))\n", "print(h0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def erreur_minimale(eps, M):\n", " return 3 / 2 * math.sqrt(M * eps)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "erreur_minimale(2 ** (-52), math.exp(1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et voici l'approximation de la dérivée obtenue en prenant $h=h_0$, le $h$ théoriquement optimal." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "approx_deriv(math.exp, 1, h0) - math.exp(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.5 Tracés" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Traçons sur un même graphe l'erreur réelle et l'erreur théorique, pour la dérivée de l'exponentielle en 1.\n", "\n", "Remarquez la valeur $M_2=e$. En réalité $M_2=e^{1+h}>e$, il faudrait prendre un peu plus que $e$, mais allons-y comme ça !" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def erreur_theorique(h, eps, M):\n", " return h * M / 2 + 2 * eps / h" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def tracer_approx2(f, x, f1, fonc_approx, fonc_erreur, M):\n", " hs = [10. ** (-k / 50) for k in range(500)]\n", " ds = [abs(fonc_approx(f, x, h) - f1(x)) for h in hs]\n", " es1 = [fonc_erreur(h, 2 ** (-52), M) for h in hs]\n", " plt.loglog(hs, es1, 'r')\n", " plt.loglog(hs, ds, 'k')\n", " plt.grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tracer_approx2(math.exp, 1, math.exp, approx_deriv, erreur_theorique, math.exp(1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour $h$ petit, c'est le terme $\\frac{2\\varepsilon}{h}$ qui est prépondérant. Pour $h$ grand, au contraire, c'est le terme $\\frac 1 2 M_2 h$ qui domine. Remarquez, pour $h$ voisin de $1$, la courbe noire l'égèrement au-dessus de la courbe rouge. C'est dû à notre sous-estimation de la valeur de $M_2$ !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Reprenez la cellule ci-dessus avec\n", "\n", "- $f(x)=x^3$ en $x=1$\n", "- $f(x)=\\arctan x$ en $x=1$\n", "- $f(x)=\\ln x$ en $x=2$\n", "- $f(x)=\\sin x$ en $x=\\frac \\pi 4$\n", "\n", "Des surprises ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Une meilleure approximation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Combinons ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Peut-on obtenir de meilleures approximations de $f'(x)$ ? Oui, on peut. Supposons que $f$ est de classe $\\mathcal C^3$ au voisinage de $x$. Soit $h>0$ assez petit. On a, par la formule de Taylor-Lagrange,\n", "\n", "$$f(x+h)=f(x)+hf'(x)+\\frac {h^2} 2 f''(x)+\\frac{h^3}6 f'''(c_1)$$\n", "\n", "et \n", "\n", "$$f(x-h)=f(x)-hf'(x)+\\frac {h^2} 2 f''(x)-\\frac{h^3}6 f'''(c_2)$$\n", "\n", "où $x0$ est un majorant de $|\\varepsilon_1|$ et $|\\varepsilon_2|$. Posons\n", "\n", "$$\\varphi(h)=\\frac{h^2}{6}M_3+\\frac{\\varepsilon}{h}$$\n", "\n", "On a\n", "\n", "$$\\varphi'(h)=\\frac{h}{3}M_3-\\frac{\\varepsilon}{h^2}$$\n", "\n", "Cette fois-ci, le minimum de $\\varphi$ est atteint pour\n", "\n", "$$h=\\left(\\frac{3\\varepsilon}{M_3}\\right)^{1/3}$$\n", "\n", "et la valeur du minimum est\n", "\n", "$$\\frac 1 2 3^{2/3}M^{1/3}\\varepsilon^{2/3}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Tracés" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Traçons, comme dans le paragraphe précédent." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def erreur_theorique2(h, eps, M):\n", " return h ** 2 * M / 6 + eps / h" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tracer_approx2(math.exp, 1, math.exp, approx_deriv2, erreur_theorique2, math.exp(1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cette fois-ci, le $h$ optimal est de l'ordre de $10^{-5}$, et l'erreur minimale de l'ordre de $10^{-11}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def h_optimal2(eps, M):\n", " return (3 * eps / M) ** (1 / 3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(h_optimal2(2 ** (-52), math.exp(1)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def erreur_minimale2(eps, M):\n", " return 1/ 3 * 3 ** (2 / 3) * M ** (1 / 3) * eps ** (2 / 3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "erreur_minimale2(2 ** (-52), math.exp(1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons gagné 3 chiffres significatifs par rapport à la méthode naïve ! Peut-on faire mieux ? Oui." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Une infinité de formules d'approximation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mais pourquoi une infinité ? Deux ça ne suffit pas :-) ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Taylor-Lagrange, encore" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $f$ une fonction de classe $\\mathcal C^{n+1}$ au voisinage du réel $x$. Soit $h>0$. Évaluons $f(x +j_0h)$, $f(x +(j_0+1)h)$, $\\ldots$, $f(x +(j_0+n)h)$. Posons dans ce qui suit $n_j=j_0+j$.\n", "\n", "$$f(x+n_jh)=\\sum_{k=0}^n\\frac{n_j^kh^k}{k!}f^{(k)}(x)+\\frac{n_j^{n+1}h^{n+1}}{(n+1)!}f^{(n+1)}(c_j)$$\n", "\n", "Combinons ensuite toutes ces valeurs. Pour cela, donnons nous $n+1$ réels $\\lambda_j$, $j=0,\\ldots,n$. On a\n", "\n", "$$\\sum_{j=0}^n\\lambda_jf(x+n_jh)=\\sum_{k=0}^na_k\\frac{h^k}{k!}f^{(k)}(x)+\\frac{a_{n+1}h^{n+1}}{(n+1)!}$$\n", "\n", "où, pour $k=0,1,\\ldots,n$,\n", "\n", "$$a_k=\\sum_{j=0}^n n_j^k\\lambda_j$$\n", "\n", "et\n", "\n", "$$a_{n+1}=\\sum_{j=0}^nn_j^{n+1}\\lambda_jf^{(n+1)}(c_j)$$\n", "\n", "Est-il possible de choisir les $\\lambda_j$ pour que $a_0=0$, $a_1=1$, $a_2=\\ldots=a_n=0$ ? Admettons pour l'instant que cela est possible. On obtient alors\n", "\n", "$$\\sum_{j=0}^n\\lambda_jf(x+n_jh)=hf'(x)+\\frac{a_{n+1}h^{n+1}}{(n+1)!}$$\n", "\n", "Et, donc, la promesse d'une formule nous donnant une valeur approchée de $f'(x)$ !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Remarque__ : Dans tout ce qui suit je prendrai toujours en exemple des $n_j$ symétriques par rapport à 0. Par exemple :\n", "\n", "- $j_0=-1$, $n=2$ : les $n_j$ sont $-1,0,1$.\n", "- $j_0=-2$, $n=4$ : les $n_j$ sont $-2, -1,0,1,2$.\n", "- $j_0=-3$, $n=6$ : les $n_j$ sont $-3,-2, -1,0,1,2,3$.\n", "- etc.\n", "\n", "Cela n'a rien d'obligatoire, rien ne vous empêche de faire des tests différents de ceux-ci !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Vandermonde !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons repoussé à plus tard dans le paragraphe précédent le problème suivant : le système\n", "\n", "$$(S)\\quad a_k=\\sum_{j=0}^n n_j^k\\lambda_j$$\n", "\n", "où $0\\le k\\le n$ a-t-il une solution $(\\lambda_0,\\ldots,\\lambda_n)$ ? Nous avons là un système de $n+1$ équations à $n+1$ inconnues dont la matrice $A$ est la matrice dont les coefficients sont les $n_j^k$, $0\\le k\\le n, 0\\le j\\le n$ ($k$ est l'indice des lignes, $j$ est celui des colonnes). Cette matrice est bien connue, c'est une matrice de Vandermonde. Et comme les $n_j$ sont des nombres distincts, la théorie nous dit que $A$ est inversible. En d'autres termes $(S)$ est un système de Cramer. D'où la\n", "\n", "__Proposition__ : Le système $(S)$ a une unique solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.3 Résolvons le système" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `vandermonde` prend en paramètre une liste $s$ de réels et renvoie la matrice de Vandermonde associée." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def vandermonde(s):\n", " n = len(s)\n", " A = [[x ** k for x in s] for k in range(n)]\n", " return Matrix(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vandermonde([-2, -1, 0, 1, 2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il est maintenant facile de résoudre notre système, qui est de la forme $A\\Lambda=B$. Son unique solution est\n", "\n", "$$\\Lambda=A^{-1}B$$\n", "\n", "La fonction `coefs` prend en paramètres les entiers $j_0$ et $n$ dont nous parlons depuis un certain temps et renvoie la liste des coefficients $\\lambda_j,j=0,\\ldots,n$. Évidemment, `sympy` est notre ami :-). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def coefs(j0, n):\n", " s = [j0 + j for j in range(n + 1)]\n", " A = vandermonde(s)\n", " B = (n + 1) * [0]\n", " B[1] = 1\n", " B = Matrix(B)\n", " return A ** (-1) * B" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Remarque__ : le calcul des coefficients $\\lambda_j$ nécessite l'inversion _exacte_ d'une matrice de taille $n+1$. Il nous faudra donc être raisonnables sur la valeur de $n$. En pratique, les temps de calcul deviennent prohibitifs lorsque $n$ atteint 20." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.4 Exemples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Trois exemples ? Reprenons tout d'abord l'approximation naïve du début du notebook,\n", "\n", "$$f'(x)\\simeq\\frac{f(x+h)-f(x)}{h}$$\n", "\n", "Ici, $j_0=0$ et $n=1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = coefs(0, 1)\n", "C" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous trouvons effectivement $\\lambda_0=-1$ et $\\lambda_1=1$.\n", "\n", "Prenons maintenant l'approximation un peu moins naïve\n", "\n", "$$f'(x)\\simeq\\frac{f(x+h)-f(x-h)}{2h}$$\n", "\n", "Ici, $j_0=-1$ et $n=2$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = coefs(-1, 2)\n", "C" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Prenons maintenant $j_0=-2$ et $n=4$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = coefs(-2, 4)\n", "C" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Je vous laisse écrire la formule correspondante. Ou alors attendez un peu, Python va le faire pour nous." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.5 La \"formule\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici notre nouvelle fonction pour approcher une dérivée. La fonction `approx_deriv3` prend en paramètres\n", "\n", "- Une fonction $f$.\n", "- Un réel $x$.\n", "- Une liste $C$ de coefficients.\n", "- Un entier $j_0$.\n", "- Un réel $h>0$.\n", "\n", "Elle renvoie l'approximation\n", "\n", "$$f'(x)\\simeq\\frac 1 h\\sum_{j=0}^n\\lambda_j f(x+n_jh)$$\n", "\n", "où $n+1$ est la longueur de $C$ et $n_j=j_0+j$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def approx_deriv3(f, x, C, j0, h):\n", " s = 0\n", " n = len(C) - 1\n", " for j in range(0, n + 1):\n", " s += C[j] * f(x + h * (j0 + j))\n", " return s / h" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "approx_deriv3(math.exp, 1, coefs(-2, 4), -2, 1e-3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voulez-vous quelque chose de plus \"formel\" ? Voici la \"formule d'approximation\"." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = Function('f')\n", "x, h = symbols('x h')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "simplify(approx_deriv3(f, x, coefs(-2, 4), -2, h))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons donc\n", "\n", "$$f'(x)\\simeq\\frac{f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)}{12h}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.6 L'erreur commise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3.6.1 Quelques tests" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reprenons les représentations graphiques de l'erreur que nous avions faites au début du notebook. Si nous prenons le troisième exemple vu ci-dessous, voici ce que nous obtenons." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "j0 = -2\n", "n = 4\n", "C = coefs(j0, n)\n", "tracer_approx(math.exp, 1, math.exp, lambda f, x, h: approx_deriv3(f, x, C, j0, h))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'erreur minimale est obtenue pour $h\\simeq10^{-3}$ et elle est de l'ordre de $10^{-13}$. Encore deux chiffres significatifs gagnés ! Pleins d'espoir, tentons $j_0=-5$ et $n=10$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "j0 = -5\n", "n = 10\n", "C = coefs(j0, n)\n", "tracer_approx(math.exp, 1, math.exp, lambda f, x, h: approx_deriv3(f, x, C, j0, h))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Une erreur minimale inférieure à $10^{-14}$, obtenue pour $h\\simeq\\frac 1 {10}$ ... Sachant que Python travaille avec un peu plus de 15 chiffres significatifs, nous sommes à un chiffre du mieux possible !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3.6.2 La théorie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rappelons la formule trouvée plus haut :\n", "\n", "$$\\sum_{j=0}^n\\lambda_jf(x+n_jh)=hf'(x)+\\frac{a_{n+1}h^{n+1}}{(n+1)!}$$\n", "\n", "Si nous tenons compte des erreurs d'arrondis dues aux calculs en virgule flottante, il nous faut remplacer $f(x+n_jh)$ par $f(x+n_jh)+\\varepsilon_j$, où $\\varepsilon_j\\in\\mathbb R$ est la fameuse erreur. Ainsi,\n", "\n", "$$\\sum_{j=0}^n\\lambda_jf(x+n_jh)+ \\sum_{j=0}^n\\lambda_j\\varepsilon_j=hf'(x)+\\frac{a_{n+1}h^{n}}{(n+1)!}$$\n", "\n", "Si l'on suppose tous les $\\varepsilon_j$ inférieurs, en valeur absolue, à un réel $\\varepsilon > 0$, il vient\n", "\n", "$$\\left|f'(x)-\\frac 1 h\\sum_{j=0}^n\\lambda_jf(x+n_jh)\\right|\\le \\frac \\varepsilon h\\sum_{j=0}^n|\\lambda_j|+\\frac{|a_{n+1}|h^{n}}{(n+1)!}$$\n", "\n", "Rappelons-nous que\n", "\n", "$$a_{n+1}=\\sum_{j=0}^nn_j^{n+1}\\lambda_jf^{(n+1)}(c_j)$$\n", "\n", "On a donc \n", "\n", "$$|a_{n+1}|\\le M_{n+1}\\sum_{j=0}^n|n_j|^{n+1}|\\lambda_j|$$\n", "\n", "où $M_{n+1}$ est un majorant de $|f^{(n+1)}|$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Notation__ : Pour $k\\in\\mathbb N$, notons $S_k=\\sum_{j=0}^n|\\lambda_j||n_j|^k$. La fonction `sum_abs` ci-dessous calcule les sommes $S_k$. On lui passe évidemment en paramètre la liste $C$ des coefficients $\\lambda_j$, ainsi que l'entier $j_0$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def sum_abs(C, j0, k):\n", " s = 0\n", " n = len(C) - 1\n", " for j in range(0, n + 1):\n", " nj = j0 + j\n", " s += abs(C[j]) * abs(nj ** k)\n", " return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = coefs(-2, 4)\n", "C, [sum_abs(C, -2, k) for k in range(6)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Avec nos nouvelles notations nous avons donc\n", "\n", "$$\\left|f'(x)-\\frac 1 h\\sum_{j=0}^n\\lambda_jf(x+n_jh)\\right|\\le \\frac {\\varepsilon S_0} h+\\frac{M_{n+1}S_{n+1}h^{n}}{(n+1)!}$$\n", "\n", "Posons \n", "\n", "$$\\varphi(h)=\\frac {\\varepsilon S_0}h+\\frac{M_{n+1}S_{n+1}h^{n}}{(n+1)!}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Un peu de courage, nous y sommes presque (non, c'est pas vrai :-))." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def fact(n):\n", " p = 1\n", " for k in range(1, n + 1): p *= k\n", " return p" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def erreur_theorique3(C, j0, h, eps, M):\n", " n = len(C) - 1\n", " S0 = sum_abs(C, j0, 0)\n", " S1 = sum_abs(C, j0, n + 1)\n", " return eps * S0 / h + M * S1 * h ** n / fact(n + 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici le graphe de l'erreur théorique." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = coefs(-2, 4)\n", "j0 = -2\n", "hs = [10. ** (-k / 20) for k in range(250)]\n", "es1 = [erreur_theorique3(C, j0, h, 2 ** (-52), math.exp(1)) for h in hs]\n", "plt.loglog(hs, es1, 'r')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Une étude de $\\varphi$ comme celles que nous avons déjà faites plus haut montre que le minimum de $\\varphi$ est atteint pour\n", "\n", "$$h_0=\\left(\\frac{\\varepsilon S_0(n+1)!}{nM_{n+1}S_{n+1}}\\right)^{1/(n+1)}$$\n", "\n", "et que\n", "\n", "$$\\varphi(h_0)=\\varepsilon S_0\\left(1+\\frac 1 n\\right)\\left(\\frac{nM_{n+1}S_{n+1}}{\\varepsilon S_0(n+1)!}\\right)^{1/(n+1)}$$\n", "\n", "__Exercice__ : Faites les calculs ! Non, ce n'est pas difficile." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def h_optimal3(C, j0, eps, M):\n", " n = len(C) - 1\n", " S0 = sum_abs(C, j0, 0)\n", " S1 = sum_abs(C, j0, n + 1)\n", " return (eps * S0 * fact(n + 1) / (n * M * S1)) ** (1 / (n + 1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = coefs(-2, 4)\n", "h_optimal3(C, -2, 2 ** (-52), math.exp(1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def erreur_minimale3(C, j0, eps, M):\n", " n = len(C) - 1\n", " S0 = sum_abs(C, j0, 0)\n", " S1 = sum_abs(C, j0, n + 1)\n", " return eps * S0 * (1 + 1 / n) * ((n * M * S1) / (eps * S0 * fact(n + 1))) ** (1 / (n + 1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = coefs(-2, 4)\n", "erreur_minimale3(C, -2, 2 ** (-52), math.exp(1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Superposons comme d'habitude l'erreur réelle et son estimation théorique :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "j0 = -2\n", "n = 4\n", "C = coefs(j0, n)\n", "tracer_approx2(math.exp, 1, math.exp, \n", " lambda f, x, h: approx_deriv3(f, x, C, j0, h), \n", " lambda h, eps, M:erreur_theorique3(C, j0, h, eps, M), \n", " math.exp(1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et tentons pour finir une formule faisant intervenir 19 points !" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C = coefs(-9, 18)\n", "erreur_minimale3(C, -2, 2 ** (-52), math.exp(1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Juste pour rire, voici la formule : " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "simplify(approx_deriv3(f, x, C, -9, h))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "j0 = -9\n", "n = 18\n", "#C = coefs(j0, n)\n", "tracer_approx2(math.exp, 1, math.exp, \n", " lambda f, x, h: approx_deriv3(f, x, C, j0, h), \n", " lambda h, eps, M:erreur_theorique3(C, j0, h, eps, M), \n", " math.exp(1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Dérivées d'ordre supérieur" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour ne pas faire fuir un éventuel lecteur, je ne reprendrai pas dans cette section les calculs d'erreurs effectués dans les sections précédentes. Nous allons juste voir que notre méthode permet de calculer des valeurs approchées de dérivées d'ordre 2, 3, etc." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rappelons-nous la formule\n", "\n", "$$\\sum_{j=0}^n\\lambda_jf(x+n_jh)=\\sum_{k=0}^na_k\\frac{h^k}{k!}f^{(k)}(x)+\\frac{a_{n+1}h^{n+1}}{(n+1)!}$$\n", "\n", "où, pour $k=0,1,\\ldots,n$,\n", "\n", "$$a_k=\\sum_{j=0}^n n_j^k\\lambda_j$$\n", "\n", "et\n", "\n", "$$a_{n+1}=\\sum_{j=0}^nn_j^{n+1}\\lambda_jf^{(n+1)}(c_j)$$\n", "\n", "Nous avions alors décidé de choisir les $\\lambda_j$ pour que $a_1=1$ et $a_0=a_2=\\ldots=a_n=0$.\n", "\n", "Soit maintenant $p\\in[0,n]$. Nous pouvons choisir les $\\lambda_j$ pour que $a_p=1$ et tous les autres $a_k$ soient nuls, $0\\le k\\le n$. On obtient alors\n", "\n", "$$\\sum_{j=0}^n\\lambda_jf(x+n_jh)=\\frac{h^p}{p!}f^{(p)}(x)+\\frac{a_{n+1}h^{n+1}}{(n+1)!}$$\n", "\n", "Une valeur approchée de $f^{(p)}(x)$ est donc\n", "\n", "$$f^{(p)}(x)\\simeq \\frac {p!}{h^p}\\sum_{j=0}^n\\lambda_jf(x+n_jh)$$\n", "\n", "avec une erreur commise que nous n'étudierons pas. Pas parce que c'est difficile mais parce que c'est exactement comme nous avons déjà fait. Laissé au lecteur, sans état d'âme !" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def coefs2(j0, n, p):\n", " s = [j0 + j for j in range(n + 1)]\n", " A = vandermonde(s)\n", " B = (n + 1) * [0]\n", " B[p] = 1\n", " B = Matrix(B)\n", " return A ** (-1) * B" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def approx_deriv4(f, x, p, C, j0, h):\n", " s = 0\n", " n = len(C) - 1\n", " for j in range(0, n + 1):\n", " s += C[j] * f(x + h * (j0 + j))\n", " return s * fact(p) / h ** p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici par exemple une formule pour calculer une dérivée seconde." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "simplify(approx_deriv4(f, x, 2, coefs2(-1, 2, 2), -1, h))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "j0 = -1\n", "n = 2\n", "p = 2\n", "C = coefs2(j0, n, p)\n", "tracer_approx(math.exp, 1, math.exp, lambda f, x, h: approx_deriv4(f, x, p, C, j0, h))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'erreur minimale est de l'ordre de $10^{-7}$. Approcher une dérivée d'ordre supérieur est-il une tâche délicate ?\n", "\n", "__Exercice__ : Mettez ci-dessus $j_0$ à $-8$ et $n$ à $16$. Moralité : tout n'est pas perdu :-)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On tente une dérivée d'ordre 3 ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "factor(approx_deriv4(f, x, 3, coefs2(-2, 4, 3), -2, h))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "j0 = -2\n", "n = 4\n", "p = 3\n", "C = coefs2(j0, n, p)\n", "hs = [10. ** (-k / 40) for k in range(300)]\n", "ds = [abs(approx_deriv4(math.exp, 1, 3, C, j0, h) - math.exp(1)) for h in hs]\n", "plt.loglog(hs, ds, 'k')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Et maintenant ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Maintenant une question essentielle se pose : avez-vous appris quelque chose en lisant ce notebook ? Il y a une excellente façon d'en être sûrs :\n", "\n", "1. Calculez un majorant de l'erreur théorique, $\\varphi(h)$, commise dans le paragraphe 4 lors du calcul de $f^{(p)}(x)$. Cette erreur dépend de paramètres que nous appellerons sans plus de détails $M_{n+1}$, $\\varepsilon$ et $p$. Sans compter les coefficients $\\lambda_k$, bien entendu. Pour cela inspirez-vous du paragraphe 3. Les calculs sont totalement analogues. Allez, dans un élan de bonté, je vous donne la réponse. Vous devriez trouver\n", "\n", "$$\\varphi(h)=\\frac{p!S_{n+1}M_{n+1}}{(n+1)!}h^{n+1-p}+\\frac{p!S_0\\varepsilon}{h^p}$$\n", "\n", "2. Calculez la valeur $h_0$ du paramètre $h$ qui minimise $\\varphi$. Vous trouvez\n", "\n", "$$h_0=\\left(\\frac{p(n+1)!S_0\\varepsilon}{(n+1-p)S_{n+1}M_{n+1}}\\right)^{\\frac 1 {n+1}}$$\n", "\n", "3. Calculez $\\varphi(h_0)$. Si vous avez été soigneux, vous avez évidemment trouvé\n", "\n", "$$\\frac{p!S_0\\varepsilon(n+1)}{n+1-p}\\left(\\frac{(n+1-p)S_{n+1}M_{n+1}}{p(n+1)!S_0\\varepsilon}\\right)^{\\frac p {n+1}}$$\n", "\n", "4. Tracez sur un même graphe l'erreur théorique en rouge et l'erreur réelle en noir, et admirez.\n", "\n", "Avoir la possibilité de __visualiser__ des résultats est un avantage immense. La confrontation entre théorie et pratique permet de se convaincre que nos raisonnements, nos calculs, ne sont pas entâchés (de trop) d'erreurs, que nos théories sont pertinentes. Si vous voulez un petit encouragement pour faire le travail demandé, peut-être faudra-t-il que je vous avoue l'inavouable ? Ce que le scientifique ne dira pas (car, c'est bien connu, le scientifique ne se trompe jamais) ? \n", "\n", "__Avant d'obtenir les splendides et parfaits graphiques de ce notebook j'ai dû recommencer plusieurs fois tous mes calculs__ :-).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }