diff --git a/1.1 Dichotomie.ipynb b/1.1 Dichotomie.ipynb index bb560f5..bece09a 100644 --- a/1.1 Dichotomie.ipynb +++ b/1.1 Dichotomie.ipynb @@ -1,272 +1,279 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Equations non-linéaires\n", "\n", "**Objectif :** trouver les zéros (ou racines) d'une fonction $f:[a,b] \\rightarrow \\mathbf R$ : \n", "$$\\alpha \\in [a,b] \\,:\\, f(\\alpha) = 0$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# importing libraries used in this book\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# defining the fonction\n", "def f(x):\n", " return x*np.sin(x*2.*np.pi) + 0.5*x - 0.25\n", "\n", "[a,b] = [-2,2]\n", "\n", "# points used to plot the graph \n", "z = np.linspace(a, b, 100)\n", "\n", "plt.plot(z, f(z),'-')\n", "\n", "# labels, title, legend\n", "plt.xlabel('x'); plt.ylabel('$f(x)$'); #plt.title('data')\n", "plt.legend(['$f(t)$'])\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Méthode de dichotomie ou bissection\n", "\n", "Si $f$ est continue et elle change de signe dans $[a,b]$, alors il existe au moins un $\\alpha$ tel que $f(\\alpha) = 0$.\n", "\n", "On peut alors définir l'algorithme suivant :\n", "\n", "$a^{(0)}=a$, $b^{(0)}=b$. Pour $k=0,1,...$\n", "\n", "1. $x^{(k)}=\\frac{a^{(k)}+b^{(k)}}{2}$\n", "2. si $f(x^{(k)})=0$, alors $x^{(k)}$ est le zéro cherché.\n", "\n", " Autrement:\n", "\n", " 1. soit $f(x^{(k)})f(a^{(k)})<0$, alors le zéro\n", " $\\alpha\\in[a^{(k)},x^{(k)}]$. \n", "\n", " On pose $a^{(k+1)}=a^{(k)}$ et $b^{(k+1)}=x^{(k)}$\n", "\n", " 2. soit $f(x^{(k)}f(b^{(k)})<0$, alors le zéro\n", " $\\alpha\\in[x^{(k)},b^{(k)}]$. \n", "\n", " On pose $a^{(k+1)}=x^{(k)}$ et $b^{(k+1)}=b^{(k)}$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercice\n", - "Ecrivez une fonction qui effectue l'algorithme de dichotomie ayant la structure suivante:\n", - "```python\n", - "def bisection(a,b,f,tolerance,maxIterations) :\n", - " # [a,b] interval of interest\n", - " # f function\n", - " # tolerance desired accuracy\n", - " # maxIterations : maximum number of iteration\n", - " # returns:\n", - " # zero, residual, number of iterations\n", - "```\n", + "Comprenez et completez la fonction suivante qui effectue l'algorithme de dichotomie \n", "\n", "Ensuite testez-la pour trouver la racine de la fonction $f(x) = x\\sin(2\\pi x) + \\frac12 x - \\frac14$ dans l'intervalle $[-1.5,1]$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def bisection(a,b,fun,tolerance,maxIterations) :\n", " # [a,b] interval of interest\n", " # fun function\n", " # tolerance desired accuracy\n", " # maxIterations : maximum number of iteration\n", " # returns:\n", " # zero, residual, number of iterations\n", " \n", " if (a >= b) :\n", " print(' b must be greater than a (b > a)')\n", " return 0,0,0\n", " \n", " # what we consider as \"zero\"\n", " eps = 1e-12\n", " \n", " # evaluate f at the endpoints\n", " fa = fun(a)\n", " fb = fun(b)\n", "\n", " if abs(fa) < eps : # a is the solution\n", " zero = a\n", " esterr = fa\n", " k = 0\n", " return zero, esterr, k\n", " \n", " if abs(fb) < eps : # b is the solution\n", " zero = b\n", " esterr = fb\n", " k = 0\n", " return zero, esterr, k\n", "\n", " if fa*fb > 0 :\n", " print(' The sign of FUN at the extrema of the interval must be different')\n", " return 0,0,0\n", "\n", "\n", "\n", " # We want the final error to be smaller than tol, \n", " # i.e. k > log( (b-a)/tol ) / log(2)\n", "\n", " nmax = int(np.ceil(np.log( (b-a)/tol ) / np.log(2)))\n", "\n", " \n", " # but nmax shall be smaller the the nmaximum iterations asked by the user\n", " if ( maxIterations < nmax ) :\n", " nmax = int(round(maxIterations))\n", " print('Warning: nmax is smaller than the minimum number of iterations necessary to reach the tolerance wished');\n", "\n", " # vector of intermadiate approximations etc\n", " x = np.zeros(nmax)\n", "\n", " # initial error is the length of the interval.\n", " esterr = (b - a)\n", "\n", " # do not need to store all the a^k and b^k, so I call them with a new variable name:\n", " ak = a\n", " bk = b\n", " # the values of f at those points are fa and fk\n", "\n", " for k in range(nmax) :\n", "\n", " # approximate solution is midpoint of current interval\n", - " x[k] = (ak + bk) / 2\n", - " fx = fun(x[k]);\n", + " # COMPLETE the code below\n", + " #> x[k] = \n", + " #> fx = \n", + " \n", " # error estimator is the half of the previous error\n", - " esterr = esterr / 2\n", + " #> esterr =\n", + "\n", + " # COMPLETE the code below\n", + "\n", + " # if we found the solution, stop the algorithm\n", + " if np.abs(fx) < eps :\n", + " # error is zero\n", + " #> zero = \n", + " #> esterr = \n", + " #> return \n", "\n", - " # (Complete the code below)\n", + " if fx*fa < 0 : # alpha is in (a,x)\n", + " #> bk = \n", + " elif fx*fb < 0 : # alpha is in (x,b)\n", + " #> ak = \n", + " else :\n", + " error('Algorithm not operating correctly')\n", "\n", "\n", "\n", " zero = x[k];\n", "\n", " if esterr > tol :\n", " print('Warning: bisection stopped without converging to the desired tolerance because the maximum number of iterations was reached');\n", "\n", " return zero, esterr, k, x \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[a,b] = [-.5,1]\n", "tol = 1e-10\n", "maxIter = 4\n", "zero, esterr, k, x = bisection(a,b,f,tol,maxIter)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plotBisectionIterations (a,b,f,x,N=200) :\n", " # plot the graphical interpretation of the Bisection method\n", " import matplotlib.pyplot as plt\n", "\n", " def putSign(y,f,text) :\n", " if (f(y) < 0) :\n", " plt.annotate(text, (y, 0.02*deltaAxis) )\n", " else :\n", " plt.annotate(text, (y, -0.05*deltaAxis) )\n", " \n", " \n", " z = np.linspace(a,b,N)\n", " plt.plot(z,f(z), 'b-')\n", " \n", " plt.plot([a,a], [0,f(a)], 'g:')\n", " plt.plot([b,b], [0,f(b)], 'g:')\n", "\n", " # For putting a sign at the initial point\n", " deltaAxis = plt.gca().axes.get_ylim()[1] - plt.gca().axes.get_ylim()[0]\n", "\n", " putSign(a,f,\"a\")\n", " putSign(b,f,\"b\")\n", " \n", " for k in range(x.size) :\n", " plt.plot([x[k],x[k]], [0, f(x[k])], 'g:')\n", " putSign(x[k],f,\"$x_\"+str(k)+\"$\")\n", "\n", " # Plot the x,y-axis \n", " plt.plot([a,b], [0,0], 'k-',linewidth=0.1)\n", " plt.plot([0,0], [np.min(f(z)),np.max(f(z))], 'k-',linewidth=0.1)\n", "\n", " plt.ylabel('$f$'); plt.xlabel('$x$');\n", "\n", " return" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams.update({'font.size': 16})\n", "plt.figure(figsize=(8, 5))\n", "\n", "plotBisectionIterations(a,b,f,x)\n", "\n", "# plt.savefig('Bisection-iterations.png', dpi=600)\n", "\n", "print(f'The estimated root is {zero:2.12f}, the estimated error {esterr:1.3e} and the residual is {f(zero):1.3e}, after {k} iterations')" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/1.2 NetwonRaphson.ipynb b/1.2 NetwonRaphson.ipynb index aedbcac..a1db454 100644 --- a/1.2 NetwonRaphson.ipynb +++ b/1.2 NetwonRaphson.ipynb @@ -1,457 +1,247 @@ { "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# importing libraries used in this book\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Méthode de Newton\n", "\n", "Soit $f:\\mathbb R \\rightarrow\\mathbb R$ une fonction différentiable.\n", "\n", "Soit $x^{(0)}$ un point donné. On considère l'équation de la droite $y(x)$ qui\n", "passe par le point $(x^{(k)},f(x^{(k)}))$ et qui a comme pente\n", "$f'(x^{(k)})$,\n", "\\begin{equation*}\n", " y(x)=f'(x^{(k)})(x-x^{(k)})+f(x^{(k)}).\n", "\\end{equation*}\n", "On définit $x^{(k+1)}$ comme étant le point où cette droite\n", "intersecte l'axe $x$, c'est-à-dire $y(x^{(k+1)})=0$. On en\n", "déduit que :\n", "$$ x^{(k+1)}=x^{(k)}-\\frac{f(x^{(k)})}{f'(x^{(k)})},\\,\\,\\,k=0,1,2\\ldots .$$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Exercice (4, Série 11)\n", + "#### Exercice (5, Série 1)\n", "\n", "On cherche les zéros de la fonction\n", "\n", "$$f(x) = \\dfrac{1}{2} \\sin \\left( \\dfrac{\\pi x}{2} \\right) + 1 - x \\; .$$\n", "\n", "\n", "1. Vérifiez qu'il y a au moins un zéro $\\alpha$ dans l'intervalle $[0,2]$.\n", "\n", "2. Ecrivez la méthode de Newton pour trouver le zéro $\\alpha$ de la fonction $f(x)$ et\n", "calculez la première itération à partir de la valeur initiale $x^{(0)}=1$.\n", "\n", "3. Calculez les zéros $\\alpha$ de la fonction\n", " $f$ avec la méthode de Newton (fonction `newton` que vous devrez écrire)\n", "\n", "```python\n", "def Newton( F, dF, x0, tol, nmax ) :\n", " # NEWTON Find the zeros of a nonlinear equations.\n", " # NEWTON(F,DF,X0,TOL,NMAX) tries to find the zero X of the \n", " # continuous and differentiable function F nearest to X0 using \n", " # the Newton method. DF is a function which take X and return the derivative of F.\n", " # If the search fails an error message is displayed.\n", " # \n", " # returns the value of the\n", " # residual R in X,the number of iterations N required for computing X and\n", " # INC the increments computed by Newton.\n", " \n", " return x, r, n, inc\n", "```\n", "\n", "Choisissez $x^{(0)} = 1$ comme point de départ pour la méthode\n", "et utilisez une tolérance $tol=10^{-4}$ sur la valeur absolue de l'incrément\n", "entre deux itérations\n", "successives $|x^{(k+1)}-x^{(k)}|$. \n", "\n", "*Dans le cas de la méthode de Newton,\n", "l'incrément est une bonne approximation de l'erreur.*\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def Newton( F, dF, x0, tol, nmax ) :\n", " '''\n", " NEWTON Find the zeros of a nonlinear equations.\n", " NEWTON(F,DF,X0,TOL,NMAX) tries to find the zero X of the \n", " continuous and differentiable function F nearest to X0 using \n", " the Newton method. DF is a function which take X and return the derivative of F.\n", " If the search fails an error message is displayed.\n", " \n", " Outputs : [x, r, n, inc, x_sequence]\n", " x : the approximated root of the function\n", " r : the absolute value of the residualin X\n", " n : the number of iterations N required for computing X and\n", " inc : the increments computed by Newton.\n", " x_sequence : the sequence computed by Newton\n", " '''\n", " \n", " # Initial values\n", " n = 0\n", " xk = x0\n", "\n", " # initialisation of loop components\n", " # increments (in abs value) at each iteration\n", " inc = []\n", " # in case we wish to plot the sequence \n", " x = [x0]\n", " \n", " # (Complete the code below)\n", " \n", " # Warning if not converged\n", " if n > nmax :\n", " print('Newton stopped without converging to the desired tolerance ')\n", " print('because the maximum number of iterations was reached')\n", " \n", " # (Complete the code below)\n", " \n", " \n", " return xk1, rk1, n, inc, np.array(x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = lambda x : 0.5*np.sin(np.pi*x/2)+1-x\n", "df = lambda x : 0.25*np.pi*np.cos(np.pi*x/2)-1\n", "\n", "x0 = 1\n", "tol = 1e-4\n", "nmax = 10\n", "\n", "zero, residual, niter, inc, x = Newton(f, df, x0, tol, nmax)\n", "\n", "print(f'The zero computed is {zero:1.4f}')\n", "print(f'Newton stoppedconverged in {niter} iterations'); \n", "print(f'with a residual of {residual:1.4e}.\\n');\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from NonLinearEquationsLib import plotNewtonIterations" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[a,b] = [0,3.5]\n", "x0 = 3\n", "zero, residual, niter, inc, x = Newton(f, df, x0, tol, nmax)\n", "\n", "# Rezise plots, which are usually too small\n", "plt.figure(figsize=(12, 4))\n", "plt.rcParams.update({'font.size': 12})\n", "\n", "# Subplot 1 over 2, 1st one\n", "plt.subplot(121)\n", "\n", "#plt.plot(range(MaxIterations), RelativeError, 'b:.')\n", "plt.plot(range(niter), inc, 'b:.')\n", " \n", "plt.xlabel('n'); plt.ylabel('$\\\\delta x$');\n", "plt.grid(True)\n", "#plt.xscale('log') \n", "plt.yscale('log')\n", "plt.legend(['$|\\\\delta x|$'])\n", "\n", "\n", "# Subplot 1 over 2, 2nd one\n", "plt.subplot(122)\n", "\n", "plotNewtonIterations (a,b,f,x,200)\n", "\n", "plt.show()\n", "\n", "# Rezise plots, which are usually too small\n", "plt.figure(figsize=(8, 4))\n", "plt.rcParams.update({'font.size': 16})\n", "\n", "plotNewtonIterations (a,b,f,x,200)\n", "# plt.savefig('Newton-iterations.png', dpi=600)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[a,b] = [-1,3.5]\n", "z = np.linspace(a,b,200)\n", "plt.plot(z,f(z), 'b-', x[6],f(x[6]), 'rx')\n", "plt.ylabel('$f(x)$'); plt.xlabel('$x$');\n", "\n", "# Plot the x,y-axis \n", "plt.plot([a,b], [0,0], 'k-',linewidth=0.1)\n", "plt.plot([0,0], [np.min(f(z)),np.max(f(z))], 'k-',linewidth=0.1)\n", "\n", "plt.legend(['$f$','$\\\\alpha$'])\n", "# plt.savefig('Newton-fx-alpha.png', dpi=600)\n", "\n", "plt.show()\n", "\n" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Exercice (2, Série 11)\n", - "\n", - "On considère les méthodes de point fixe $x^{(n+1)} = g_i(x^{(n)})\n", - "\\quad (i=1,2,3)$ avec:\n", - "\n", - "$$g_1(x^{(n)}) = \\dfrac{1}{2} e^{x^{(n)}/2} , \\qquad g_2(x^{(n)}) = -\\dfrac{1}{2}\n", - " e^{x^{(n)}/2}, \\qquad g_3(x^{(n)}) = 2 \\ln(2x^{(n)}),$$\n", - "dont les fonctions d'itération $g_i(x)$ sont visualisées sur la figure plus bas\n", - "\n", - "\n", - "1. Pour chaque point fixe $\\bar{x}$ de la fonction d'itération $g_i$ ($i=1,2,3$), on suppose d'avoir choisi une valeur \n", - "initiale $x^{(0)}$ proche de $\\bar{x}$. Etudiez si la méthode converge vers $\\bar{x}$.\n", - "\n", - "2. Pour chaque fonction d'itération $g_i$, déterminez graphiquement \n", - "pour quelles valeurs initiales $x^{(0)}$ la méthode de point fixe correspondante \n", - "converge et vers quel point fixe. \n", - "\n", - "3. Montrez que si $\\bar{x}$ est un point fixe de la fonction $g_i$ \n", - " ($i=1,2,3$), alors il est aussi un zéro de la fonction \n", - " $f(x)= e^x - 4x^2$ (dont le comportement est tracé sur la dernière figure). \n", - "\n", - "4. Comment peut-on calculer les zéros de $f$? \n", - "\n", - "*L'exercice 2 est à faire sur papier, ici une indication par ordinateur*\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from NonLinearEquationsLib import plotPhi, FixedPoint, plotPhiIterations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.rcParams['figure.figsize'] = [20, 5]\n", - "\n", - "plt.subplot(1,3,1)\n", - "[a,b] = [-2,5]\n", - "phi1 = lambda x : np.exp(x/2)/2\n", - "plotPhi(a,b,phi1,'$g_1$')\n", - "\n", - "plt.subplot(1,3,2)\n", - "[a,b] = [-2,5]\n", - "phi2 = lambda x : - np.exp(x/2)/2\n", - "plotPhi(a,b,phi2,'$g_2$')\n", - "\n", - "plt.subplot(1,3,3)\n", - "[a,b] = [1e-1,5]\n", - "phi3 = lambda x : 2*np.log(2*x)\n", - "plotPhi(a,b,phi3,'$g_3$')\n", - "\n", - "plt.show()\n", - "\n", - "# Graph of the fonction $f(x)=e^x-4x^2$\n", - "N = 100\n", - "z = np.linspace(a,b,N)\n", - "f = lambda x : np.exp(x) - 4*x*x\n", - "\n", - "plt.subplot(1,3,1)\n", - "plt.plot(z,f(z),'k-')\n", - "\n", - "plt.xlabel('x'); plt.ylabel('f(x)');\n", - "# Plot the x,y-axis \n", - "plt.plot([a,b], [0,0], 'k-',linewidth=0.1)\n", - "plt.plot([0,0], [np.min(f(z)),np.max(f(z))], 'k-',linewidth=0.1)\n", - "plt.legend(['f(x)'])\n", - "plt.title('Graph of $f(x)=e^x-4x^2$')\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Partie 1**\n", - "*Pour chaque point fixe $\\bar{x}$ de la fonction d'itération $g_i$ ($i=1,2,3$), on suppose choisir une valeur \n", - "initiale $x^{(0)}$ proche de $\\bar{x}$. Etudiez si la méthode converge vers $\\bar{x}$.*\n", - "\n", - "On va utiliser la fonction `FixedPoint` qui se trouve dans `NonLinearEquationsLib.py`\n", - "```python\n", - "def FixedPoint( phi, x0, a, b tol, nmax ) :\n", - " '''\n", - " FixedPoint Find the fixed point of a function by iterative iterations\n", - " FixedPoint( PHI,X0,a,b, TOL,NMAX) tries to find the fixedpoint X of the a\n", - " continuous function PHI nearest to X0 using \n", - " the fixed point iterations method. \n", - " [a,b] : if the iterations exit the interval, the method stops\n", - " If the search fails an error message is displayed.\n", - " \n", - " Outputs : [x, r, n, inc, x_sequence]\n", - " x : the approximated fixed point of the function\n", - " r : the absolute value of the residual in X : |phi(x) - x|\n", - " n : the number of iterations N required for computing X and\n", - " x_sequence : the sequence computed by Newton\n", - " \n", - " ...\n", - " return xk1, rk1, n, np.array(x) \n", - " '''\n", - "```" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tol = 1e-2\n", - "nmax = 10\n", - "\n", - "# Choose fonction phi\n", - "[a,b] = [-2,5]\n", - "phi = phi1\n", - "label = '$phi_1$'\n", - "\n", - "# Initial Point\n", - "x0 = 3\n", - "zero, residual, niter, x = FixedPoint(phi, x0, a,b, tol, nmax)\n", - "\n", - "plt.subplot(131)\n", - "plotPhi (a,b,phi,label)\n", - "plt.plot(x,phi(x), 'rx')\n", - "\n", - "# plot the graphical interpretation of the Fixed Point method\n", - "plotPhiIterations(x)\n", - "\n", - "plt.show()\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tol = 1e-2\n", - "nmax = 10\n", - "\n", - "# Choose fonction phi\n", - "intervals = [ [-2,5] , [-2,5] , [1e-2,5] ]\n", - "phiFunctions = [phi1, phi2, phi3]\n", - "labels = ['$phi_1$', '$phi_2$', '$phi_3$']\n", - "# Initial Points\n", - "initialPoints = [4.2,2,1]\n", - "\n", - "alpha = np.empty(3)\n", - "\n", - "for k in range(3) :\n", - " phi = phiFunctions[k]; x0 = initialPoints[k]\n", - " a = intervals[k][0]; b = intervals[k][1]\n", - " label = labels[k]\n", - " \n", - " alpha[k], residual, niter, x = FixedPoint(phi, x0, a,b, tol, nmax)\n", - "\n", - " plt.subplot(1,3,k+1)\n", - "\n", - " plotPhi (a,b,phi,label[k])\n", - " plt.plot(x,phi(x), 'rx')\n", - "\n", - " # plot the graphical interpretation of the Fixed Point method\n", - " plotPhiIterations(x)\n", - " \n", - " \n", - "\n", - "plt.show()\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Graph of the fonction $f(x)=e^x-4x^2$\n", - "N = 100\n", - "a,b = [-2,5]\n", - "z = np.linspace(a,b,N)\n", - "f = lambda x : np.exp(x) - 4*x*x\n", - "\n", - "plt.subplot(1,3,1)\n", - "plt.plot(z,f(z),'k-')\n", - "\n", - "# Solutions found:\n", - "plt.plot(alpha,f(alpha),'ro')\n", - "plt.annotate(\"$\\\\alpha_1$\", (alpha[0], 2))\n", - "plt.annotate(\"$\\\\alpha_2$\", (alpha[1], 2))\n", - "plt.annotate(\"$\\\\alpha_3$\", (alpha[2]+0.1, -4))\n", - "\n", - "plt.xlabel('x'); plt.ylabel('f(x)');\n", - "# Plot the x,y-axis \n", - "plt.plot([a,b], [0,0], 'k-',linewidth=0.1)\n", - "plt.plot([0,0], [np.min(f(z)),np.max(f(z))], 'k-',linewidth=0.1)\n", - "plt.legend(['f(x)','$\\\\alpha$'])\n", - "plt.title('Graph of $f(x)=e^x-4x^2$')\n", - "\n", - "plt.show()" - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "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.7.6" + "version": "3.7.7" } }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/2.1 Interpolation.ipynb b/2.1 Interpolation.ipynb index bc8f1e1..701f05e 100644 --- a/2.1 Interpolation.ipynb +++ b/2.1 Interpolation.ipynb @@ -1,347 +1,347 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 1 Interpolation et approximation de données\n", "\n", "## 1.1 Position du problème" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpolation de données\n", "\n", "Soit $n \\geq 0$ un nombre entier. Etant donnés $(n+1)$ noeuds \n", "distincts $x_0$, $x_1$,$\\dots$ $x_n$ et $(n+1)$ valeurs $y_0$,\n", "$y_1$,$\\dots$ $y_n$, on cherche un polynôme $p$\n", " de degré $n$, tel que\n", "\n", "$$p(x_j)=y_j \\qquad \\text{ pour } \\, 0\\leq j\\leq n.$$\n", "\n", "**Exemple** On cherche le polynôme $\\Pi_n$ de degré $n=4$ tel que $\\Pi_n(x_j) = y_j, j =1,...,5$ avec \n", "les données suivantes \n", "\n", "| $x_k$ | $y_k$ |\n", "| --- | --- |\n", "| 1 | 3 |\n", "| 1.5 | 4 |\n", "| 2 | 2 |\n", "| 2.5 | 5 |\n", "| 3 | 1 |\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "# importing libraries used in this book\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "# Some data given: x=1, 1.5, 2, 2.5, 3 and y = 3,4,2,5,1 \n", "x = np.linspace(1, 3, 5) # equivalent to np.array([ 1, 1.5, 2, 2.5, 3 ])\n", "y = np.array([3, 4, 2, 5, 1])\n", "\n", "# Plot the points using matplotlib\n", "plt.plot(x, y, 'ro')\n", "\n", "plt.xlabel('x'); plt.ylabel('y'); #plt.title('data')\n", "plt.legend(['data'])\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si ce polynôme existe, on le note $p=\\Pi_n$. On appelle $\\Pi_n$ le\n", "polynôme d'interpolation des valeurs $y_j$ aux noeuds $x_j,\\, j=0,\\dots,n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the interpolating function\n", "\n", "# Defining the polynomial function \n", "def p(x):\n", " # coefficients of the interpolating polynomial\n", " a = np.array([-140.0, 343.0, -872./3., 104.0, -40./3.])\n", " \n", " # value of the polynomial in all the points t\n", " return a[0] + a[1]*x + a[2]*(x**2) + a[3]*(x**3) + a[4]*(x**4)\n", "\n", "\n", "# points used to plot the graph \n", "z = np.linspace(1, 3, 100)\n", "\n", "plt.plot(x, y, 'ro', z, p(z))\n", "plt.xlabel('x'); plt.ylabel('y'); #plt.title('data')\n", "plt.legend(['data','$\\Pi_2(x)$'])\n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpolation de fonctions\n", "\n", "\n", "Soit $f\\in C^0(I)$ et $x_0,\\ldots,x_n\\in I$. \n", "Si on prend $$y_j=f(x_j),\\quad 0\\leq j\\leq n,$$ \n", "alors le polynôme d'interpolation\n", "$\\Pi_n(x)$ est noté $\\Pi_n f(x)$ et est appelé l'interpolant de $f$ aux\n", "noeuds $x_0$,$\\dots$ $x_n$.\n", "\n", "**Exemple** Soient $$x_1=1, x_2=1.75, x_3=2.5, x_4=3.25, x_5=4$$ les points d'interpolation et $$f(x) = x \\sin(2\\pi x).$$ On cherche l'interpolant $\\Pi_n f$ de degré $n=4$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# defining the fonction that we want to interpolate\n", "def f(x):\n", " return x*np.sin(x*2.*np.pi) \n", "\n", "# The interpolation must occour at points x=1, 1.75, 2.5, 3.25, 4\n", "x = np.linspace(1, 4, 5)\n", "\n", "# points used to plot the graph \n", "z = np.linspace(1, 4, 100)\n", "\n", "plt.plot(x, f(x), 'ro', z, f(z),':')\n", "\n", "# labels, title, legend\n", "plt.xlabel('x'); plt.ylabel('$f(x)$'); #plt.title('data')\n", "plt.legend(['$f(x_k)$','$f(x)$'])\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the interpolating function\n", "\n", "# Defining the polynomial function \n", "def p(x):\n", " # coefficients of the interpolating polynomial\n", " a = np.array([0, 7.9012, -13.037, 5.9259, -0.79012])\n", " \n", " # value of the polynomial in all the points x\n", " return a[0] + a[1]*x + a[2]*(x**2) + a[3]*(x**3) + a[4]*(x**4)\n", "\n", "\n", "# points where to evaluate the polynomial\n", "z = np.linspace(1, 4, 100)\n", "\n", "plt.plot(x, f(x), 'ro', z, f(z),':', z, p(z))\n", "plt.xlabel('$x$'); plt.ylabel('$f(x)$'); #plt.title('data')\n", "plt.legend(['$f(x_k)$','$f(x)$','$\\Pi_n(x)$'])\n", "\n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Matrice de Vandermonde\n", "\n", "Il est possible d'écrire un système d'équations et de trouver les coefficients de manière directe.\n", "Ce n'est pas toujours la meilleure solution.\n", "\n", "Nous cherchons les coefficients du polynôme $p(x) = a_0 + a_1 x + ... + a_n x^n$ qui satisfont les $(n+1)$ équations $p(x_k) = y_k, k=0,...,n$, c'est-à-dire\n", "\n", "$$a_0 + a_1 x_k + ... + a_n x_k^n = y_k, k=0,...,n$$\n", "\n", "Ce système s'écrit sous forme matricielle\n", "\n", "$$\\begin{pmatrix}\n", "1 & x_0 & x_0^2 & \\cdots & x_0^n \\\\\n", "\\vdots & & & & \\vdots \\\\\n", "1 & x_n & x_n^2 & \\cdots & x_n^n\n", "\\end{pmatrix}\n", "\\begin{pmatrix} a_0 \\\\ \\vdots \\\\ a_n \\end{pmatrix}\n", "=\\begin{pmatrix} y_0 \\\\ \\vdots \\\\ y_n \\end{pmatrix}$$\n", "\n", "Pour construire cette matrice, vous pouvez utiliser la fonction\n", "```python\n", "# Defining the mxn Vandermonde matrix \n", "def VandermondeMatrix(x):\n", " # Input\n", " # x : +1 array with interpolation nodes\n", " # Output\n", " # Matrix of Vandermonde of size n x n\n", "```\n", "que vous pouvez importer avec la commande\n", "```python\n", "from InterpolationLib import VandermondeMatrix \n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exemple** On cherche les coefficients du polynôme d'interpolation de degré $n=4$ \n", "des valeurs suivantes \n", "\n", "| $x_k$ | $y_k$ |\n", "| --- | --- |\n", "| 1 | 3 |\n", "| 1.5 | 4 |\n", "| 2 | 2 |\n", "| 2.5 | 5 |\n", "| 3 | 1 |\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from InterpolationLib import VandermondeMatrix \n", "\n", "# Some data given: x=1, 1.5, 2, 2.5, 3 and y = 3,4,2,5,1 \n", "x = np.linspace(1, 3, 5)\n", "y = np.array([3, 4, 2, 5, 1])\n", "n = x.size - 1\n", "\n", "A = VandermondeMatrix(x)\n", "# print(A)\n", "\n", "# compute coefficients\n", "# Resouds Ax = b avec b=y et rends x\n", "# >>> COMPLETE HERE <<<\n", "\n", "# print the coefficients on screen\n", "print('The coefficients a_0, ..., a_n are')\n", "print(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Now we can define the polynomial\n", "p = lambda x : a[0] + a[1]*x + a[2]*(x**2) + a[3]*(x**3) + a[4]*(x**4)\n", "\n", "# points used to plot the graph \n", "z = np.linspace(1, 3, 100)\n", "\n", "# Plot points and function\n", "# >>> COMPLETE HERE <<<\n", "plt.xlabel('x'); plt.ylabel('y'); #plt.title('data')\n", "plt.legend(['data','p(x)'])\n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Alternatives : polyfit et polyval\n", "\n", "Les fonctions `polyfit` et `polyval` de `numpy` font essentiellement la même chose que les paragraphes ci-dessous. Plus tard nous verrons des méthodes plus performantes.\n", "\n", "`a = numpy.polyfit(x, y, n, ... )` :\n", "\n", " * input : $x,y$ les données à interpoler, $n$ le degré du polynôme recherché\n", " * output : les coefficients du polynôme, _dans l'ordre inverse_ de ce que nous avons vu !\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Some data given: x=1, 1.5, 2, 2.5, 3 and y = 3,4,2,5,1 \n", "x = np.linspace(1, 3, 5)\n", "y = np.array([3, 4, 2, 5, 1])\n", "n = x.size - 1\n", "\n", "a = np.polyfit(x,y,n)\n", "\n", "# Now we can define the polynomial, with coeffs in the reverse order !\n", "p = lambda x : a[4] + a[3]*x + a[2]*(x**2) + a[1]*(x**3) + a[0]*(x**4)\n", "\n", "# We can also use polyval instead !\n", "# np.polyval(a,x)\n", "\n", "# points used to plot the graph \n", "z = np.linspace(1, 3, 100)\n", "\n", "plt.plot(x, y, 'ro', z, p(z), '.', z, np.polyval(a,z))\n", "plt.xlabel('x'); plt.ylabel('y'); #plt.title('data')\n", "plt.legend(['data','p(x)','polyval'])\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/2.2 Interpolation de Lagrange.ipynb b/2.2 Interpolation de Lagrange.ipynb index 9559b95..e8a3190 100644 --- a/2.2 Interpolation de Lagrange.ipynb +++ b/2.2 Interpolation de Lagrange.ipynb @@ -1,191 +1,191 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2 Interpolation de Lagrange\n", "\n", "### Base de Lagrange\n", "\n", "On considère les polynômes $\\varphi_k$, $k=0,\\dots, n$ de\n", "degré $n$ tels que\n", "\n", "$$\\varphi_k(x_j)=\\delta_{jk}, \\qquad k,j=0,\\dots, n,$$\n", "\n", "où $\\delta_{jk}=1$ si $j=k$ et $\\delta_{jk}=0$ si $j\\neq k$.\n", "Explicitement, on a\n", "\n", "$$\\varphi_k(x)=\\prod_{j=0,j\\ne k}^n\\dfrac{(x-x_j)}{(x_k-x_j)}.$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# importing libraries used in this book\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercice (théorique)** Vérifiez que\n", "\n", "1. $B = \\{\\varphi_k, k=0,\\dots, n\\}$ est une base de $P_n(\\mathbb R)$\n", "2. Chaque polynôme $\\varphi_k$ est de degré $n$\n", "3. $\\varphi_k(x_j)=\\delta_{jk}, \\qquad k,j=0,\\dots, n$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Defining the Lagrange basis functions\n", "def phi(x,k,z):\n", " # the input variables are:\n", " # x : the interpolatory points\n", " # k : which basis function\n", " # z : where to evaluate the function\n", " \n", "\n", " # careful, there are n+1 interpolation points!\n", " n = x.size - 1 \n", "\n", " # init result to one, of same type and size as z\n", " result = np.zeros_like(z) + 1\n", "\n", " # first few checks on k:\n", " if (type(k) != int) or (x.size < 1) or (k > n) or (k < 0):\n", " raise ValueError('Lagrange basis needs a positive integer k, smaller than the size of x')\n", " \n", " # loop on n to compute the product\n", " for j in range(0,n+1) :\n", " if (j == k) :\n", " continue\n", " if (x[k] == x[j]) :\n", " raise ValueError('Lagrange basis: all the interpolation points need to be distinct')\n", "\n", " # >>> COMPLETE HERE <<<\n", " result = \n", "\n", " return result\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exemple\n", "\n", "Pour $n=2$, $x_0=-1$, $x_1=0$, $x_2=1$, les polynômes de la base de Lagrange\n", "sont\n", "\n", "\\begin{equation*}\n", "\\begin{array}{lcl}\n", "\\varphi_0(x)&=&\\displaystyle\\dfrac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}=\n", "\\dfrac{1}{2}x(x-1),\\\\[2mm]\n", "\\varphi_1(x)&=&\\displaystyle\\dfrac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}=\n", "-(x+1)(x-1),\\\\[2mm]\n", "\\varphi_2(x)&=&\\displaystyle\\dfrac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}=\n", "\\dfrac{1}{2}x(x+1).\n", "\\end{array}\n", "\\end{equation*}\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the Lagrange Basis functions \n", "x = np.linspace(-1., 1, 3)\n", "z = np.linspace(-1.1, 1.1, 100)\n", "\n", "plt.plot(z, phi(x,0,z), 'g', z, phi(x,1,z), 'r', z, phi(x,2,z),':')\n", "\n", "plt.xlabel('x'); plt.ylabel('$\\\\varphi_{k}(x)$'); plt.title('Lagrange basis functions')\n", "plt.legend(['$\\\\varphi_{0}$','$\\\\varphi_{1}$','$\\\\varphi_{2}$'])\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercice\n", "\n", "Visualisez la base de Lagrange associée aux points $x_k = 1, 1.5, 2, 2.5, 3$.\n", "Evaluez sur le graphique les valeurs de $\\varphi_k(x_j)$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the Lagrange Basis functions \n", "x = np.linspace(1., 3, 5)\n", "z = np.linspace(0.9, 3.1, 100)\n", "\n", "plt.plot(z, phi(x,0,z), 'g', z, phi(x,1,z), 'r', z, phi(x,2,z),':', z, phi(x,3,z),':', z, phi(x,4,z),':')\n", "\n", "plt.xlabel('x'); plt.ylabel('phi(x)'); plt.title('Lagrange Basis functions')\n", "plt.legend(['$\\\\varphi_0$','$\\\\varphi_1$','$\\\\varphi_2$',\n", " '$\\\\varphi_3$','$\\\\varphi_4$'])\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Polynôme d'interpolation\n", "\n", "Le polynôme d'interpolation $\\Pi_n$ des\n", "valeurs $y_j$ aux noeuds $x_j$, $j=0,\\dots,n$, s'écrit\n", "\\begin{equation}\\label{e:int_lagr}\n", " \\Pi_n(x)=\\sum_{k=0}^n y_k \\varphi_k(x),\n", "\\end{equation}\n", "car il vérifie $\\Pi_n(x_j)=\\sum_{k=0}^n y_k \\varphi_k(x_j)=y_j$.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/2.3 Interpolation de fonction.ipynb b/2.3 Interpolation de fonction.ipynb index b2ce518..01a6abb 100644 --- a/2.3 Interpolation de fonction.ipynb +++ b/2.3 Interpolation de fonction.ipynb @@ -1,401 +1,401 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.3 Interpolation d'une fonction continue\n", "\n", "Soit $f: \\mathbf [a,b] \\rightarrow R$ continue et $x_0,\\ldots,x_n\\in [a,b]$ des noeuds distincts. Le polynôme d'interpolation\n", "$\\Pi_n(x)$ est noté $\\Pi_n f(x)$ et est appelé l'interpolant de $f$ aux\n", "noeuds $x_0,\\dots,x_n$.\n", "\n", "\n", "Si on prend $$y_k=f(x_k), k= 0,..., n,$$ \n", "alors on aura\n", "\\begin{equation*}\n", "\\Pi_nf(x)=\\sum_{k=0}^n\n", "f(x_k)\\varphi_k(x).\n", "\\end{equation*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercice\n", "\n", "Ecrivez une fonction Python qui a la définition suivante, en utilisant la fonction `phi` définie plus haut.\n", "Ecrivez aussi un petit test sur la base de l'exercice précédent.\n", "\n", "```python\n", "# Lagrange Interpolation of data (x,y), evaluated at ordinate(s) z\n", "def LagrangePi(x,y,z):\n", " # the input variables are:\n", " # x : the interpolatory points\n", " # y : the corresponding data at the points x\n", " # z : where to evaluate the function\n", " \n", "```\n", "\n", "Utilisez le fait que $\\{\\varphi_k, k = 0,...,n\\}$ est une base des polynômes de degré $\\leq n$\n", "et que le vecteur $y$ représente les coordonnées du polynôme d'interpolation recherché par rapport à cette base,\n", "c'est-à-dire\n", "$$\\Pi_n (z) = y_0 \\varphi_0 + ... + y_n \\varphi_n$$\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# importing libraries used in this book\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from InterpolationLib import LagrangeBasis as phi" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Lagrange Interpolation of data (x,y), evaluated at ordinate(s) z\n", "def LagrangePi(x,y,z):\n", " # the input variables are:\n", " # x : the interpolatory points\n", " # y : the corresponding data at the points x\n", " # z : where to evaluate the function\n", " \n", " # {phi(x,k,.), k=0,...,n} is a basis of the polynomials of degree n\n", " # y represents the coordinates of the interpolating polynomial with respect to this basis.\n", " # Therefore LagrangePi(x,y,.) = y[0] phi(x,0,.) + ... + y[n] phi(x,n,.)\n", "\n", " # >>> COMPLETE HERE <<<\n", "\n", " return result" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# vecteur des points d'interpolation\n", "x = np.linspace(0, 1, 5)\n", "\n", "# vecteur des valeurs\n", "y = np.array([3.38, 3.86, 3.85, 3.59, 3.49])\n", "\n", "z = np.linspace(-0.1, 1.1, 100)\n", "plt.plot(x, y, 'ro', z, LagrangePi(x,y,z))\n", "\n", "plt.xlabel('x'); plt.ylabel('y');\n", "plt.legend(['data','p(x)'])\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Erreur d'interpolation\n", "\n", "Soient $x_0$, $x_1$, $\\ldots$, $x_n$,\n", "$(n+1)$ nœuds *équirépartis* dans $I=[a,b]$ et soit $f\\in C^{n+1}(I)$.\n", "Alors\n", "\\begin{equation}\n", " \\max_{x\\in I} |f(x) - \\Pi_n f(x)| \\leq \\dfrac{1}{4(n+1)}\n", "\\left(\\frac{b-a}{n}\\right)^{n+1}\\max_{x\\in I}|f^{(n+1)}(x)|.\n", "\\end{equation}\n", "\n", "On remarque que l'erreur d'interpolation dépend de la dérivée\n", "$(n+1)$-ième de $f$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercice\n", "\n", "On considère les points d'interpolation $$x_0=1, x_1=1.75, x_2=2.5, x_3=3.25, x_4=4$$ \n", "et la fonction $$f(x) = x \\sin(2\\pi x)$$.\n", "\n", "1. Calculez la base de Lagrange associée à ces points. D'abord sur papier, ensuite utilisez Python pour en dessiner le graphique.\n", "\n", "2. Calculez le polynôme d'interpolation $\\Pi_n$ à l'aide de la base de Lagrange. D'abord sur papier, ensuite avec Python.\n", "\n", "4. Quelle est l'erreur théorique d'interpolation ? \n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Comportement pour $n$ grand: eg, la fonction de Runge.\n", "\n", "Le fait que\n", "$$\\lim_{n\\rightarrow\\infty} \\dfrac{1}{4(n+1)}\\left(\\dfrac{b-a}{n}\\right)^{n+1}=0$$\n", "n'implique pas forcément que $\\max_{x \\in I} |E_nf(x)|$ tende vers zéro quand $n\\rightarrow\\infty$.\n", "\n", "Soit $f(x)=\\dfrac{1}{1+x^2}$, $x\\in [-5,5]$. Si on l'interpole dans des\n", "noeuds équirépartis, l'interpolant présente des oscillations au voisinage des extrémités de l'intervalle.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Runge fonction\n", "f = lambda x : 1./(1+x**2)\n", "\n", "# Values of N to use\n", "Nrange = [3,5,10]\n", "# plotting points\n", "z = np.linspace(-5, 5, 100)\n", "\n", "for n in Nrange :\n", "\n", " # define x and y=f(x)\n", " # compute the Lagrange interpolant\n", " # plot it\n", "\n", " # >>> COMPLETE HERE <<<\n", "\n", "\n", "plt.plot(z,f(z), 'b')\n", "plt.xlabel('x'); plt.ylabel('y'); plt.title('Runge function')\n", "plt.legend(Nrange)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`sympy` est une librairie pour le calcul symbolique en Python. Nous allons l'utiliser pour étudier le comportement de l'erreur d'interpolation de la fonction de Runge." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Using symbolic python to compute derivatives\n", "import sympy as sp\n", "\n", "# define x as symbol\n", "x = sp.symbols('x')\n", "\n", "# define Runge function\n", "f = 1/(1+x**2)\n", "\n", "# pretty print the 5th derivative of f\n", "f5 = sp.diff(f, x,5)\n", "# display f5\n", "sp.init_printing(use_unicode=True)\n", "display(f5)\n", "\n", "\n", "# evalf can be used to compute the value of a function at a given point\n", "print('5th derivative evaluated at 3 :')\n", "print( f5.evalf(subs={x: 3}) )\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour définir une fonction qui accepte un array de valeur, il faut utiliser les lignes suivantes.\n", "\n", "Ensuite on peut aussi dessiner le graphe ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "# to evaluate a function at many given points, we need the following trick\n", "diff_f_func = lambda t: float(sp.diff(f,x,k).evalf(subs={x: t}))\n", "diff_f = np.vectorize(diff_f_func)\n", "\n", "# the derivative can be set with k (not very elegant...)\n", "k = 4\n", "print(diff_f(4.5))\n", "\n", "# plotting points\n", "z = np.linspace(-5, 5, 100)\n", "# plot the derivative between -5 and 5\n", "\n", "plt.plot(z,diff_f(z), 'b')\n", "plt.xlabel('t'); plt.ylabel('y'); plt.title('Derivatives of Runge function')\n", "plt.legend(Nrange)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... ou évaluer le maximum pour plusieurs $n$ et voir le comportement de $\\max |f^{(n)}|$ en fonction de $n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot max(abs(fn)) in the range -5,5\n", "z = np.linspace(-5, 5, 100)\n", "\n", "Nmax = 10\n", "maxValFn = np.zeros(Nmax)\n", "for k in range(Nmax):\n", " maxValFn[k] = np.max(np.abs(diff_f(z)))\n", " \n", "plt.plot(range(10), maxValFn)\n", "\n", "plt.yscale('log')\n", "plt.xlabel('n'); plt.ylabel('$\\max|\\partial f|$');\n", "plt.title('Max of $|\\partial^n f|$');\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpolation de Chebyshev\n", "\n", "\n", "Pour chaque entier positif $n\\geq 1$, pour $i=0,\\dots n$, on note\n", "$$\\hat{x}_i=-\\cos(\\pi i/n)\\in[-1,1]$$\n", "**les points de Chebyshev** et on définit\n", "\n", "$$x_i=\\dfrac{a+b}{2}+\\dfrac{b-a}{2}\\hat{x}_i\\in[a,b],$$\n", "\n", "pour un intervalle arbitraire $[a,b]$. Pour une fonction continue $f\\in\n", "C^1([a,b])$, le polynôme d'interpolation $\\Pi_n f$ de degré $n$ aux noeuds\n", "$\\{x_i,i=0,\\ldots,n\\}$ converge uniformément vers $f$ quand\n", "$n\\rightarrow \\infty$.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Chebichev points on the interval [-1,1]\n", "\n", "z = np.linspace(0,1, 100)\n", "plt.plot(np.cos(np.pi*z), np.sin(np.pi*z))\n", "\n", "n =5\n", "z = np.linspace(0,1, n+1)\n", "plt.plot(np.cos(np.pi*z), np.sin(np.pi*z), 'o')\n", "plt.plot(np.cos(np.pi*z), 0*z, 'x')\n", "\n", "for k in range(0,n+1) :\n", " plt.plot([np.cos(np.pi*z[k]),np.cos(np.pi*z[k])],[0,np.sin(np.pi*z[k])],':')\n", "\n", "plt.axis('equal')\n", "\n", "plt.xlabel('t'); \n", "plt.title('Chebyshev points')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exemple\n", "\n", "On reprend le\n", " même exemple mais on interpole la fonction de Runge dans les points de\n", "Chebyshev. La figure montre les polynômes de\n", "Chebyshev de degrés $n=5$ et $n=10$. On remarque que les oscillations\n", "diminuent lorsqu'on augmente le degré du polynôme.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Runge fonction\n", "f = lambda x : 1./(1+x**2)\n", "\n", "# Values of N to use\n", "Nrange = [3,5,10]\n", "# plotting points\n", "[a,b] = [-5,5]\n", "z = np.linspace(a,b, 100)\n", "\n", "for n in Nrange :\n", "\n", " # Chebyshev points on [-1,1]\n", " hx = -np.cos(np.pi*np.linspace(0,n,n+1)/n)\n", " # mapped to [a,b]\n", " x =(a+b)/2 + (b-a)/2*hx\n", "\n", " y = f(x);\n", "\n", " plt.plot(z, LagrangePi(x,y,z), ':')\n", "\n", "plt.plot(z,f(z), 'b')\n", "plt.xlabel('x'); plt.ylabel('y'); plt.title('Chebyshev interpolation')\n", "plt.legend(Nrange)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/2.4 Interpolation par intervalles.ipynb b/2.4 Interpolation par intervalles.ipynb index b5c8760..8ba8535 100644 --- a/2.4 Interpolation par intervalles.ipynb +++ b/2.4 Interpolation par intervalles.ipynb @@ -1,287 +1,287 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.4 Interpolation par intervalles\n", "\n", "### Interpolation linéaire par morceaux\n", "\n", "Soit $f: \\mathbf [a,b] \\rightarrow R$ continue et $a=x_0<\\ldots>> COMPLETE HERE <<<\n", "\n", "\n", "\n", "plt.xlabel('x'); plt.ylabel('y');\n", "plt.legend(['data','p(x)'])\n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercice\n", "\n", "Dessinez le graphe de la fonction de Runge et l'interpolateur linéaire par morceaux \n", "sur l'intervalle $[-5,5]$ pour $N=3,5,10$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# interval and function\n", "a = -5; b = 5\n", "f = lambda x : 1/(1+x**2)\n", "\n", "# >>> COMPLETE HERE <<<\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Erreur d'interpolation linéaire par morceaux \n", "\n", "#### Théorème \n", "\n", "Soient $f\\in C^{2}([a,b])$, $H = \\frac{b-a}{N}$, $x_i = a + iH$ pour $i=0,...,N$.\n", " \n", "Soit $E^H_1f(x) = f(x) - \\Pi_1^{H} f(x)$, alors\n", "$$\\max_{x\\in I}\\mid E^H_1f(x) \\mid\\leq \\dfrac{H^2}{8}\\max_{x \\in I} | f''(x)|.$$\n", "\n", "**Preuve**\n", "D'après la formule, sur chaque intervalle $I_i=[x_i,x_{i+1}]$, on a\n", " \n", "$$\\max_{x\\in[x_i,x_{i+1}]}\\mid E^H_1f(x)\\mid \\leq\n", " \\dfrac{H^2}{4(1+1)}\\max_{x\\in I_i}\\mid f''(x)\\mid.$$\n", "\n", "**Remarque**\n", "On peut aussi montrer que, si l'on utilise un polynôme\n", "de degré $n$ ($\\geq 1$) et si l'on dénote $E^H_nf(x) = f(x) - \\Pi_n^{H} f(x)$,\n", " dans chaque sous-intervalle $I_i$, on\n", "trouve\n", "\n", "$$\\max_{x\\in I}\\mid E^H_nf(x) \\mid \\leq \\dfrac{H^{n+1}}{4(n+1)} \\max_{x \\in I} |f^{(n+1)} (x)| \\, .$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpolation quadratique par morceaux\n", "\n", "Soit $f: \\mathbf [a,b] \\rightarrow R$ continue et $a=x_0<\\ldots$(n+1)$ points \n", "~~distincts~~ $x_0$, $x_1$,$\\dots$ $x_n$ et $(n+1)$ valeurs $y_0$,\n", "$y_1$,$\\dots$ $y_n$, on cherche un polynôme $p$\n", " de degré $m , tel que\n", "\n", "$$\\sum_{j=0}^n |y_j - p(x_j)|^2 \\qquad \\text{ soit le plus petit possible}.$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On appelle **polynôme aux moindres carrés de degré $m$**\n", "le polynôme $\\hat{p}_n(x)$ de degré $m$ tel que\n", "\\begin{equation}\n", " \\sum_{j=0}^n |y_j - \\hat{p}_n(x_j)|^2 \\leq \\sum_{j=0}^n |y_j - p_n(x_j)|^2\n", " \\qquad \\forall p_m(x) \\in \\mathbb{P}_m\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous cherchons les coefficients du polynôme $p(x) = a_0 + a_1 x + ... + a_m x^m$ qui satisfait *au mieux* les $(n+1)$ équations $p(x_k) = y_k, k=0,...,n$, c'est-à-dire\n", "\n", "$$a_0 + a_1 x_k + ... + a_m x_k^m = y_k, k=0,...,n$$\n", "\n", "Ce système s'écrit sous forme matricielle \n", "\n", "$$\\begin{pmatrix}\n", "1 & x_0 & \\cdots & x_0^m \\\\\n", "\\vdots & & & \\vdots \\\\\n", "1 & x_n & \\cdots & x_n^m\n", "\\end{pmatrix}\n", "\\begin{pmatrix} a_0 \\\\ \\vdots \\\\ a_m \\end{pmatrix}\n", "=\\begin{pmatrix} y_0 \\\\ \\vdots \\\\ y_n \\end{pmatrix}$$\n", "\n", "Puisque $m>> COMPLETE HERE <<<\n", "a = \n", "\n", "# print the coefficients on screen\n", "print('The coefficients a_0, ..., a_m are')\n", "print(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def polynomial(a,x):\n", " m = a.size-1\n", " # \\hat p = a_0 + a_1 x + ... + a_n x^m\n", " # is equal to the scalar product between the vectors a and (1, x, ..., x^m) :\n", " return np.power( np.tile(x, (m+1, 1)).T , np.linspace(0,m,m+1)).dot(a)\n", "\n", "\n", "# points used to plot the graph, slightly larger than data\n", "# >>> COMPLETE HERE <<<\n", "z = \n", "\n", "plt.plot(sigma, epsilon, 'ro', z, polynomial(a,z),'b')\n", "plt.xlabel('$\\sigma$'); plt.ylabel('$\\epsilon$');\n", "plt.legend(['data','$\\hat p_n$'])\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercice\n", "\n", "Depuis https://hsso.ch/fr/2012/b/14 on a téléchargé les données de la population suisse dans\n", "le fichier `Data/PopulationSuisse.csv`.\n", "\n", "Approximez l'évolution de la population avec un polynôme de degré $n=1,2,3,7$. \n", "\n", "Ensuite faites l'hypothèse de croissance exponentielle de la population, c'est-à-dire $p(x) = e^{a_1 x+a_0}$ où $a_0$ et $a_1$ sont des paramètres. Comment utiliser l'approximation polynômiale dans ce contexte ? " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Data of the population has been dowloaded from https://hsso.ch/fr/2012/b/14\n", "# into the file Data/PopulationSuisse.csv\n", "import pandas as pd \n", "\n", "# Read data from file 'PopulationSuisse.csv' \n", "data = pd.read_csv(\"Data/PopulationSuisse.csv\") \n", "# Preview the first 5 lines of the loaded data \n", "data.head()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = data['Année'].to_numpy()\n", "y = data['Population'].to_numpy()\n", "\n", "m = 2\n", "B = VandermondeMatrix(x,m)\n", "\n", "# >>> COMPLETE HERE <<<\n", "# ... and solve the exercise\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def polynomial(a,x):\n", " m = a.size-1\n", " # \\hat p = a_0 + a_1 x + ... + a_n x^m\n", " # is equal to the scalar product between the vectors a and (1, x, ..., x^m) :\n", " return np.power( np.tile(x, (m+1, 1)).T , np.linspace(0,m,m+1)).dot(a)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On assume une croissance exponentielle : $population (x) = C * e^{a_1x} = e^{a_o + a_1 x}$\n", "\n", "En d'autres termes: $\\log (population) (x) = a_o + a_1 x$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }