{
"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": [
"## 1.1.1 Méthode de la bissection\n",
"\n",
""
]
},
{
"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": [
"## 1.1.2 Bissection, critères d'arrêts\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"## 1.1.3 Bissection, exemple\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercice\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) - 1\n",
"\n",
" nmax = int(np.ceil(np.log( (b-a)/tol ) / np.log(2))) - 1\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+1)\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+1) :\n",
"\n",
" # approximate solution is midpoint of current interval\n",
" # COMPLETE the code below\n",
" #> x[k] = \n",
" #> fx = \n",
" \n",
" # error estimator is the half of the previous error\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",
" 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] = [-1.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": [
"from NonLinearEquationsLib import plotBisectionIterations"
]
},
{
"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 (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
}