diff --git a/1 Equations non-lineaires.ipynb b/1 Equations non-lineaires.ipynb index 9b3300b..c7ba387 100644 --- a/1 Equations non-lineaires.ipynb +++ b/1 Equations non-lineaires.ipynb @@ -1,654 +1,657 @@ { "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 that we want to interpolate\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)})= 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", " # error estimator is the half of the previous error\n", " esterr = esterr / 2\n", " \n", " # if we found the solution, stop the algorithm\n", " if np.abs(fx) < eps :\n", " # error is zero\n", " zero = x[k]\n", " esterr = 0;\n", " return zero, esterr, k \n", "\n", " if fx*fa < 0 : # alpha is in (a,x)\n", " bk = x[k]\n", " elif fx*fb < 0 : # alpha is in (x,b)\n", " ak = x[k] \n", " else :\n", " error('Algorithm not operating correctly')\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": [ "from NonLinearEquationsLib import plotBisectionIterations\n", "\n", "[a,b] = [-.5,1]\n", "tol = 1e-10\n", "maxIter = 4\n", "zero, esterr, k, x = bisection(a,b,f,tol,maxIter)\n", "\n", "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')\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", "\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", - " \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", "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 }