diff --git a/5 Equations non-lineaires.ipynb b/5 Equations non-lineaires.ipynb new file mode 100644 index 0000000..0449555 --- /dev/null +++ b/5 Equations non-lineaires.ipynb @@ -0,0 +1,627 @@ +{ + "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 \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "[a,b] = [-1.5,1]\n", + "tol = 1e-10\n", + "maxIter = 1000\n", + "zero, esterr, k = bisection(a,b,f,tol,maxIter)\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érifier qu'il y a au moins un zéro $\\alpha$ dans l'intervalle $[0,2]$.\n", + "\n", + "2. Ecrire la méthode de Newton pour trouver le zéro $\\alpha$ de la fonction $f(x)$ et\n", + "calculer la première itération à partir de la valeur initiale $x^{(0)}=1$.\n", + "\n", + "3. Calculer les zéro $\\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", + "Choisir $x^{(0)} = 1$ comme points de départ pour la méthode\n", + "et utiliser 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 approche assez bien l'erreur commise.*\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" + ] + }, + { + "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 HERE\n", + " \n", + " \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", + " 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", + "\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()" + ] + }, + { + "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 choisir une valeur \n", + "initiale $x^{(0)}$ proche de $\\bar{x}$. Etudier si la méthode converge vers $\\bar{x}$.\n", + "\n", + "2. Pour chaque fonction d'itération $g_i$, déterminer 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. Montrer 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}$. Etudier 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.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/NonLinearEquationsLib.py b/NonLinearEquationsLib.py new file mode 100644 index 0000000..0569538 --- /dev/null +++ b/NonLinearEquationsLib.py @@ -0,0 +1,229 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Jan 30 + +@author: Simone Deparis +""" + +import numpy as np + +def Bisection(a,b,fun,tolerance,maxIterations) : + # [a,b] interval of interest + # fun function + # tolerance desired accuracy + # maxIterations : maximum number of iteration + # returns: + # zero, residual, number of iterations + + if (a >= b) : + print(' b must be greater than a (b > a)') + return 0,0,0 + + # what we consider as "zero" + eps = 1e-12 + + # evaluate f at the endpoints + fa = fun(a) + fb = fun(b) + + if abs(fa) < eps : # a is the solution + zero = a + esterr = fa + k = 0 + return zero, esterr, k + + if abs(fb) < eps : # b is the solution + zero = b + esterr = fb + k = 0 + return zero, esterr, k + + if fa*fb > 0 : + print(' The sign of FUN at the extrema of the interval must be different') + return 0,0,0 + + + + # We want the final error to be smaller than tol, + # i.e. k > log( (b-a)/tol ) / log(2) + + nmax = int(np.ceil(np.log( (b-a)/tol ) / np.log(2))) + + + # but nmax shall be smaller the the nmaximum iterations asked by the user + if ( maxIterations < nmax ) : + nmax = int(round(maxIterations)) + print('Warning: nmax is smaller than the minimum number of iterations necessary to reach the tolerance wished'); + + # vector of intermadiate approximations etc + x = np.zeros(nmax) + + # initial error is the length of the interval. + esterr = (b - a) + + # do not need to store all the a^k and b^k, so I call them with a new variable name: + ak = a + bk = b + # the values of f at those points are fa and fk + + for k in range(nmax) : + + # approximate solution is midpoint of current interval + x[k] = (ak + bk) / 2 + fx = fun(x[k]); + # error estimator is the half of the previous error + esterr = esterr / 2 + + # if we found the solution, stop the algorithm + if np.abs(fx) < eps : + # error is zero + zero = x[k] + esterr = 0; + return zero, esterr, k + + if fx*fa < 0 : # alpha is in (a,x) + bk = x[k] + elif fx*fb < 0 : # alpha is in (x,b) + ak = x[k] + else : + error('Algorithm not operating correctly') + + + zero = x[k]; + + if esterr > tol : + print('Warning: bisection stopped without converging to the desired tolerance because the maximum number of iterations was reached'); + + return zero, esterr, k + + + + +def FixedPoint( phi, x0, a, b, tol, nmax ) : + ''' + FixedPoint Find the fixed point of a function by iterative iterations + FixedPoint( PHI,X0,a,b, TOL,NMAX) tries to find the fixedpoint X of the a + continuous function PHI nearest to X0 using + the fixed point iterations method. + [a,b] : if the iterations exit the interval, the method stops + If the search fails an error message is displayed. + + Outputs : [x, r, n, inc, x_sequence] + x : the approximated fixed point of the function + r : the absolute value of the residual in X : |phi(x) - x| + n : the number of iterations N required for computing X and + x_sequence : the sequence computed by Newton + ''' + + # Initial values + n = 0 + xk = x0 + + # initialisation of loop components + # increments (in abs value) at each iteration + inc = [] + # in case we wish to plot the sequence + x = [x0] + + # diff : last increment, + diff = tol + 1 # initially set larger then tolerance + + # Loop until tolerance is reached + while ( diff >= tol and n <= nmax and (xk >= a) and (xk <= b) ) : + # call phi + xk1 = phi(xk) + + # increments + diff = np.abs(xk1 - xk) + + # prepare the next loop + n = n + 1 + xk = xk1 + x.append(xk) + + # Final residual + rk1 = np.abs(xk1 - phi(xk1)) + + # Warning if not converged + if n > nmax : + print('FixexPoint stopped without converging to the desired tolerance ') + print('because the maximum number of iterations was reached') + + return xk1, rk1, n, np.array(x) + + +def plotPhi (a,b,phi,label='$\phi$', N = 100) : + ''' + simple plot of fonction with bisectrice of first quadrant + useful for preparing grpahics of fixed point iterations + [a,b] : interval, used for both x and y axis + phi : funciton to plot + label : label for the fonction, usually '$\phi$' + N : number of points for plotting + ''' + import matplotlib.pyplot as plt + + z = np.linspace(a,b,N) + plt.plot(z,phi(z),'k-') + plt.plot([a,b],[a,b],':', linewidth=0.5) + + plt.xlabel('x'); plt.ylabel(label); + # Plot the x,y-axis + plt.plot([a,b], [0,0], 'k-',linewidth=0.1) + plt.plot([0,0], [a,b], 'k-',linewidth=0.1) + plt.legend([label, 'y=x']) + plt.title('Graph de la fonction ' + label) + + +def plotPhiIterations (x) : + # plot the graphical interpretation of the Fixed Point method + import matplotlib.pyplot as plt + + plt.plot([x[0],x[0]], [0,x[1]], 'g:') + + for k in range(x.size-2) : + plt.plot([x[k],x[k+1]], [x[k+1],x[k+1]], 'g:') + plt.plot([x[k+1],x[k+1]], [x[k+1],x[k+2]], 'g:') + k=x.size-2 + plt.plot([x[k],x[k+1]], [x[k+1],x[k+1]], 'g:') + + # Putting a sign at the initial point + deltaAxis = plt.gca().axes.get_ylim()[1] - plt.gca().axes.get_ylim()[0] + if (x[1] < 0) : + plt.annotate("$x_0$", (x[0], 0.02*deltaAxis) ) + else : + plt.annotate("$x_0$", (x[0], -0.05*deltaAxis) ) + + +def plotNewtonIterations (a,b,f,x,N=200) : + # plot the graphical interpretation of the Newton method + import matplotlib.pyplot as plt + + z = np.linspace(a,b,N) + plt.plot(z,f(z), 'b-', x,f(x), 'rx') + + # Putting a sign at the initial point + deltaAxis = plt.gca().axes.get_ylim()[1] - plt.gca().axes.get_ylim()[0] + + # plot the graphical interpretation of the Newton method + plt.plot([x[0],x[0]], [0,f(x[0])], 'g:') + for k in range(x.size-1) : + plt.plot([x[k],x[k+1]], [f(x[k]),0], 'g-') + plt.plot([x[k+1],x[k+1]], [0,f(x[k+1])], 'g:') + + # Putting a sign at the initial point + if (f(x[k]) < 0) : + plt.annotate("$x_"+str(k)+"$", (x[k], 0.02*deltaAxis) ) + else : + plt.annotate("$x_"+str(k)+"$", (x[k], -0.05*deltaAxis) ) + + + plt.xlabel('f'); plt.ylabel('x'); + + # Plot the x,y-axis + plt.plot([a,b], [0,0], 'k-',linewidth=0.1) + plt.plot([0,0], [np.min(f(z)),np.max(f(z))], 'k-',linewidth=0.1) + + plt.legend(['f','(xk,fk)']) +