FEUILLE DE TP 5

Introduction aux schémas numériques pour f(x) = 0¶

Dans ce TP, nous allons implémenter les méthodes de dichotomie, de la sécante et de point fixe afin de résoudre des équations scalaires non linéaires $f(x)=0$.

In [24]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import numpy as np
import matplotlib.pyplot as plt

color_bleu = 'xkcd:dusky blue'
color_vert = 'xkcd:grassy green'
color_orange = 'xkcd:orangish'
color_rouge = 'xkcd:blood orange'
color_bleuclair = 'xkcd:sky blue'
color_vertclair = 'xkcd:grass green'
color_marron = 'xkcd:burnt sienna'
color_violet = 'xkcd:blue violet'

liste_color = [
    color_bleu,
    color_vert,
    color_orange,
    color_rouge,
    color_bleuclair,
    color_vertclair,
    color_marron,
    color_violet
]

Exercice 1 : Méthode de dichotomique¶

Pour cet exercice on fera les tests avec les fonctions suivantes sur l'intervalle $[-1,2]$ : $$ \begin{gathered} f_0(x) = e^x - 1 - \frac{x}{2} - \frac{x^2}{2},\\ f_1(x) = e^x - 1 - x,\\ f_2(x) = e^x - 1 - x - \frac{x^2}{2}. \end{gathered} $$

Question 1

  1. Commencez par implémenter ces 3 fonctions.
  2. Tracez les afin de vérifier qu'elles possèdent bien un unique 0 dans l'intervalle $[-1,2]$.
In [42]:
def f_0(x):
    return np.exp(x)-1-x*(1+x)/2

def f_1(x):
    return np.exp(x)-1-x

def f_2(x):
    return np.exp(x)-1-x*(1+x/2)

list_f = [f_0, f_1, f_2]
In [43]:
size = 5
fig = plt.figure(figsize=(size*len(list_f), size))
a, b = -1, 2
x = np.linspace(a, b, 1000)

for k, f in enumerate(list_f):
    ax = fig.add_subplot(1, len(list_f), k+1)
    ax.plot(x, f(x), color=color_rouge, linewidth=2)
    ax.set_title(f'${f.__name__}$', fontsize=20, color=color_bleu)
    ax.axhline(0, color='black', linestyle='dotted')

Question 2

Implémentez la méthode de Dichotomie à travers une fonction x, niter, aL, bL, cL = dichotomie(f, a, b, tol, itermax) qui prend en arguments

  • la fonction $f$ dont on cherche la racine,
  • l'intervalle $[a, b]$ avec $f(a) f(b) < 0$,
  • le paramètre $tol$ du test d'arrêt : on arrête si $|f(x)| < tol$
  • et le nombre $iterMax$ maximum d'itérations autorisées.

et qui retourne :

  • $x$ la solution approchée,
  • $niter$ le nombre d'itérations réalisées,
  • $aL$, $bL$ et $cL$ des ndarray contenant les termes des suites $(a_n)$, $(b_n)$ et $(c_n)$.
In [44]:
def dichotomie(f, a, b, tol=1e-6, itermax=500):
    """
    Approximation du zeros d'une fonction f passée en argument

    Parameters
    ----------

    f: function
        la fonction
    a: float
        bord gauche de l'intervalle
    b: float
        bord droit de l'intervalle
    tol: float (default 1.e-6)
        paramètre du test d'arrêt : arrêt si |f(x)| < tol ou si longueur de l'intervalle < tol
    itermax: int (default 1000)
        nombre maximal d'itération
        
    Returns
    -------
    
    x: float
        la solution approchée trouvée
    niter: int
        le nombre d'itérations
    an: ndarray
        les valeurs de la suite an
    bn: ndarray
        les valeurs de la suite bn
    cn: ndarray
        les valeurs de la suite cn
    """
    # METTRE LE CODE ICI
    an, bn = (a, b) if a < b else (b, a)
    fan, fbn = f(an), f(bn)
    la, lb, lc = [an], [bn], []
    if fan * fbn > 0:
        raise ValueError(f"Probleme dans dichotomie: {a}, {b}")
    niter, fcn = 0, 2*tol
    while bn-an > tol and abs(fcn) > tol and niter < itermax:
        cn = .5*(an+bn)
        fcn = f(cn)
        if fan*fcn <= 0:
            bn, fbn = cn, fcn
        if fbn*fcn <= 0:
            an, fan = cn, fcn
        la.append(an)
        lb.append(bn)
        lc.append(cn)
        niter += 1
    cn = .5*(an+bn)
    lc.append(cn)
    return cn, niter, np.asanyarray(la), np.asanyarray(lb), np.asanyarray(lc)

Question 3

  1. Reprenez les graphiques de la question 1.
  2. Calculez la solution approchée $x^\star$ donnée par la méthode de la dichotomie pour chacun des cas et représentez la à l'aide d'un scatter sur la figure.
  3. Sur une autre figure, tracez pour chaque fonction (sauf $f_1$ évidemment) les 3 nuages de points $(|a_n|, |a_{n+1}|)$, $(|b_n|, |b_{n+1}|)$ et $(|c_n|, |c_{n+1}|)$. Vous pouvez utiliser une échelle logarithmique ! Les courbes sont-elles cohérentes avec la théorie ? Quelle est l'ordre de la méthode ?
In [46]:
## METTRE LE CODE ICI

solution = 0
tol = 1.e-6
itermax = 5000
methode = dichotomie

size = 4
fig = plt.figure(figsize=(size*len(list_f), size*2))
fig.suptitle(f"{methode.__name__}", fontsize=20)
a, b = -1, 2
x = np.linspace(a, b, 1000)

for k, f in enumerate(list_f):
    ax = fig.add_subplot(2, len(list_f), k+1)
    ax.plot(x, f(x), color='black')
    ax.axhline(0, color='black', linestyle='dotted')
    title = f'${f.__name__}$'
    try:
        xstar, niter, aL, bL, cL = methode(f, a, b, tol=tol, itermax=itermax)
        title += f" (${niter}$ itérations)"
        dicoOK = True
    except:
        print(f"Echec de la méthode pour f_{k}")
        dicoOK = False
    ax.set_title(title, fontsize=16, color='black')
    if dicoOK:
        ax_e = fig.add_subplot(2, len(list_f), k+1+len(list_f))
        ax_e.set_xscale('log', base=2)
        ax_e.set_yscale('log', base=2)
        ax.scatter(xstar, f(xstar), s=20, color=color_rouge)
        ax_e.scatter(abs(cL[:-1]-solution), abs(cL[1:]-solution), color=color_orange, label=r'$c_n$')
        ax_e.scatter(abs(aL[:-1]-solution), abs(aL[1:]-solution), color=color_bleuclair, label=r'$a_n$')
        ax_e.scatter(abs(bL[:-1]-solution), abs(bL[1:]-solution), color=color_vertclair, label=r'$b_n$')
        ax_e.plot(abs(cL[:-1]-solution), abs(cL[:-1]-solution), color='black', alpha=0.5, label=r'pente $1$')
        ax_e.legend()
        ax_e.set_title("Erreur", fontsize=16, color='black')
        ax_e.set_xlabel(r"$|e_{n}|$", color='black')
        ax_e.set_ylabel(r"$|e_{n+1}|$", color='black')

fig.tight_layout()

Echec de la méthode pour f_1

Exercice 2 : Méthode de la sécante¶

Pour cet exercice on fera les tests avec les fonctions suivantes sur l'intervalle $[-1,1]$ : $$ \begin{gathered} f_0(x) = e^x - 1,\\ f_1(x) = e^{-x} - 1,\\ f_2(x) = 1 - e^x,\\ f_3(x) = 1 - e^{-x}. \end{gathered} $$

Question 1

  1. Commencez par implémenter ces 4 fonctions.
  2. Tracez les afin de vérifier qu'elles possèdent bien un unique 0 dans l'intervalle $[-1,1]$.
In [47]:
def f_0(x):
    return np.exp(x)-1

def f_1(x):
    return np.exp(-x)-1

def f_2(x):
    return 1-np.exp(x)

def f_3(x):
    return 1-np.exp(-x)

list_f = [f_0, f_1, f_2, f_3]
In [48]:
size = 4
fig = plt.figure(figsize=(size*len(list_f), size))
a, b = -1, 1
x = np.linspace(a, b, 1000)

for k, f in enumerate(list_f):
    ax = fig.add_subplot(1, len(list_f), k+1)
    ax.plot(x, f(x), color=color_rouge, linewidth=2)
    ax.set_title(f'${f.__name__}$', fontsize=20, color=color_bleu)
    ax.axhline(0, color='black', linestyle='dotted')

Question 2

Implémentez la méthode de la sécante à travers une fonction x, niter, aL, bL, cL = secante(f, a, b, tol, itermax) qui prend en arguments

  • la fonction $f$ dont on cherche la racine,
  • l'intervalle $[a, b]$ avec $f(a) f(b) < 0$,
  • le paramètre $tol$ du test d'arrêt : on arrête si $|f(x)| < tol$
  • et le nombre $iterMax$ maximum d'itérations autorisées.

et qui retourne :

  • $x$ la solution approchée,
  • $niter$ le nombre d'itérations réalisées,
  • $aL$, $bL$ et $cL$ des ndarray contenant les termes des suites $(a_n)$, $(b_n)$ et $(c_n)$.
In [39]:
def secante(f, a, b, tol=1e-6, itermax=500):
    """
    Approximation du zeros d'une fonction f passée en argument

    Parameters
    ----------

    f: function
        la fonction
    a: float
        bord gauche de l'intervalle
    b: float
        bord droit de l'intervalle
    tol: float (default 1.e-6)
        paramètre du test d'arrêt : arrêt si |f(x)| < tol
    itermax: int (default 1000)
        nombre maximal d'itération
        
    Returns
    -------
    
    x: float
        la solution approchée trouvée
    niter: int
        le nombre d'itérations
    an: ndarray
        les valeurs de la suite an
    bn: ndarray
    """
    # METTRE LE CODE ICI
    an, bn = (a, b) if a < b else (b, a)
    fan, fbn = f(an), f(bn)
    la, lb, lc = [an], [bn], []
    if fan * fbn > 0:
        raise ValueError(f"Probleme dans fausse position: {a}, {b}")
    niter, fcn = 0, 2*tol
    while abs(fcn) > tol and niter < itermax:
        cn = (an*fbn - bn*fan) / (fbn - fan)
        fcn = f(cn)
        if fan*fcn <= 0:
            bn, fbn = cn, fcn
        if fbn*fcn <= 0:
            an, fan = cn, fcn
        la.append(an)
        lb.append(bn)
        lc.append(cn)
        niter += 1
    cn = (an*fbn - bn*fan) / (fbn - fan)
    lc.append(cn)
    return cn, niter, np.asanyarray(la), np.asanyarray(lb), np.asanyarray(lc)

Question 3

  1. Reprenez les graphiques de la question 1.
  2. Calculez la solution approchée $x^\star$ donnée par la méthode de la fausse position pour chacun des cas et représentez la à l'aide d'un scatter sur la figure.
  3. Sur une autre figure, tracez pour chaque fonction les 3 nuages de points $(|a_n|, |a_{n+1}|)$, $(|b_n|, |b_{n+1}|)$ et $(|c_n|, |c_{n+1}|)$. Vous pouvez utiliser une échelle logarithmique ! Les courbes sont-elles cohérentes avec la théorie ? Quelle est l'ordre de la méthode ?
In [49]:
## METTRE LE CODE ICI

solution = 0
tol = 1.e-6
itermax = 5000
methode = secante

size = 4
fig = plt.figure(figsize=(size*len(list_f), size*2))
fig.suptitle(f"{methode.__name__}", fontsize=20)
a, b = -1, 2
x = np.linspace(a, b, 1000)

for k, f in enumerate(list_f):
    ax = fig.add_subplot(2, len(list_f), k+1)
    ax.plot(x, f(x), color='black')
    ax.axhline(0, color='black', linestyle='dotted')
    title = f'${f.__name__}$'
    try:
        xstar, niter, aL, bL, cL = methode(f, a, b, tol=tol, itermax=itermax)
        title += f" (${niter}$ itérations)"
        dicoOK = True
    except:
        print(f"Echec de la méthode pour f_{k}")
        dicoOK = False
    ax.set_title(title, fontsize=16, color='black')
    if dicoOK:
        ax_e = fig.add_subplot(2, len(list_f), k+1+len(list_f))
        ax_e.set_xscale('log', base=2)
        ax_e.set_yscale('log', base=2)
        ax.scatter(xstar, f(xstar), s=20, color=color_rouge)
        ax_e.scatter(abs(cL[:-1]-solution), abs(cL[1:]-solution), color=color_orange, label=r'$c_n$')
        ax_e.scatter(abs(aL[:-1]-solution), abs(aL[1:]-solution), color=color_bleuclair, label=r'$a_n$')
        ax_e.scatter(abs(bL[:-1]-solution), abs(bL[1:]-solution), color=color_vertclair, label=r'$b_n$')
        ax_e.plot(abs(cL[:-1]-solution), abs(cL[:-1]-solution), color='black', alpha=0.5, label=r'pente $1$')
        ax_e.legend()
        ax_e.set_title("Erreur", fontsize=16, color='black')
        ax_e.set_xlabel(r"$|e_{n}|$", color='black')
        ax_e.set_ylabel(r"$|e_{n+1}|$", color='black')
        
fig.tight_layout()

Question 4 Comparaison

On pourra dans cette question se contenter de comparer le nombre d'itérations.

  • Pour le cas test $\quad f(x) = e^{x} - 1 - x - 0.5x^2 \quad \text{sur} \quad [-0.9, 1]$
  • Comparer la méthode de dichotomie et celle de la sécante.
  • Que remarquez-vous ?
  • Reprendre dans le cas $\quad f(x) = e^{x} - 1 - 0.5 x - 0.5x^2 \quad \text{sur} \quad [-0.9, 1]$
In [50]:
## METTRE LE CODE ICI
list_f = [
    lambda x: np.exp(x)-1-x*(1+x)/2,
    lambda x: np.exp(x)-1-x,
    lambda x: np.exp(x)-1-x*(1+x/2),
]

tol = 1.e-6
itermax = 5000
a, b = -0.9, 1

for k, f in enumerate(list_f):
    print("-"*80)
    print(f"| fonction f{k} : ")
    print("-"*80)
    for methode in [dichotomie, fausse_position]:
        try:
            niter = methode(f, a, b, tol=tol, itermax=itermax)[1]
        except:
            niter = None
        print(f"|     Pour la méthode {methode.__name__:16}: {niter} itérations")
    print("-"*80)

--------------------------------------------------------------------------------
| fonction f0 : 
--------------------------------------------------------------------------------
|     Pour la méthode dichotomie      : 19 itérations
|     Pour la méthode fausse_position : 10 itérations
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
| fonction f1 : 
--------------------------------------------------------------------------------
|     Pour la méthode dichotomie      : None itérations
|     Pour la méthode fausse_position : None itérations
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
| fonction f2 : 
--------------------------------------------------------------------------------
|     Pour la méthode dichotomie      : 5 itérations
|     Pour la méthode fausse_position : 1920 itérations
--------------------------------------------------------------------------------

Exercice 3 : Méthode d'itération (ou de point fixe)¶

Dans cette exercice, on va utiliser les fonctions et les données initiales $x_0$ correspondantes

  • a) $\quad f(x) = e^{x} - 1 - x \quad \text{sur} \quad [-1, 1], \quad x_0 = -0.5$
  • b) $\quad f(x) = x - \sin(x) \quad \text{sur} \quad [-\frac{\pi}{2}, \frac{5\pi}{2}], \quad x_0 = 1$
  • c) $\quad f(x) = x + \sin(x) \quad \text{sur} \quad [-\frac{\pi}{2}, \frac{5\pi}{2}], \quad x_0 = 1$
  • d) $\quad f(x) = x + \cos(x) -1 \quad \text{sur} \quad [-\frac{\pi}{2}, \frac{5\pi}{2}], \quad x_0 = 1$
  • e) $\quad f(x) = x - \cos(x) +1 \quad \text{sur} \quad [-\frac{\pi}{2}, \frac{5\pi}{2}], \quad x_0 = 1$

Remarque:

  • Dans a) on cherche la racine $x^\star=0$ de $e^{x} - 1 - 2 x = 0$.
  • Dans les cas tests b), c), d), e), on cherche la racine simple $x^\star = 0$ de $\sin(x) = 0$.

Question 1

  1. Implémentez les différentes fonctions demandées dans cet exercice.
  2. Tracez dans des sous-figures séparées le graphe de chaque fonction et ajoutez la droite d'équation $y=x$.
  3. Déterminez visuellement les points fixes de chaque fonction dans l'intervalle considéré.
  4. Ajoutez un point (scatter) pour chaque point fixe.
In [51]:
list_f = [
    {
        'expr': lambda x: np.exp(x) - 1 - x,
        'a': -1, 'b': 1, 'x0': -0.5,
        'l': 0
    },
    {
        'expr': lambda x: x - np.sin(x),
        'a': -.5*np.pi, 'b': 2.5*np.pi, 'x0': 1,
        'l': 0
    },
    {
        'expr': lambda x: x + np.sin(x),
        'a': -.5*np.pi, 'b': 2.5*np.pi, 'x0': 1,
        'l': np.pi
    },
    {
        'expr': lambda x: x + np.cos(x) - 1,
        'a': -.5*np.pi, 'b': 2.5*np.pi, 'x0': 1,
        'l': 0
    },
    {
        'expr': lambda x: x - np.cos(x) + 1,
        'a': -.5*np.pi, 'b': 2.5*np.pi, 'x0': 1,
        'l': 2*np.pi
    }
]
In [52]:
size = 4
fig = plt.figure(figsize=(size*len(list_f), size))

for k, dico in enumerate(list_f):
    a, b = dico['a'], dico['b']
    x0 = dico['x0']
    x = np.linspace(a, b, 1000)
    f = dico['expr']
    ax = fig.add_subplot(1, len(list_f), k+1)
    ax.plot(x, f(x), color=color_bleu)
    ax.set_title(f'$f_{k}$', fontsize=20)
    ax.plot(x, x, color=color_vert, linestyle='dotted')
    ax.scatter([x0], [x0], s=20, color=color_bleu)
    ax.scatter([dico['l']], [dico['l']], s=20, color=color_rouge)

Question 2

Implémentez la méthode du point fixe à travers une fonction x, niter, xL = point_fixe(f, x0, tol, itermax). Les valeurs par défaut de tol et de itermax seront respectivement 1.e-6 et 5000.

In [53]:
## METTRE LE CODE ICI

def point_fixe(f, x0, tol=1.e-6, itermax=5000):
    """
    Recherche de point fixe : méthode brute x_{n+1} = f(x_n)
    
    Parameters
    ----------
    
    f: function
        la fonction dont on cherche le point fixe
    x0: float
        initialisation de la suite
    tol: float (default 1.e-6)
        critère d'arrêt : |x_{n+1} - x_n| < tol
    itermax: int
        le nombre maximal d'itérations autorisées
        
    Returns
    -------
    
    x: float
        la valeur trouvée pour le point fixe
    niter: int
        le nombre d'itérations effectuées
    xL: ndarray
        la suite des itérés de la suite
    """
    x = x0
    xL = [x]
    xold = x - 2*tol
    niter = 0
    while abs(x-xold) > tol:
        xold, x = x, f(x)
        xL.append(x)
        niter += 1
    return x, niter, np.asanyarray(xL)

Question

Pour chacun des cas tests ci-dessus faire

  • Tracez dans un même graphique la courbe de la fonction $f$ et celle de la droite $x \mapsto x$.
  • Générez la suite des itérés $x_n$ de la méthode d'itération du point fixe.
  • Représentez sur le même graphique la solution numérique $x^\star$ obtenue en dessinant un point à la position $(x^\star, f(x^\star))$
  • Ajoutez en titre le nombre d'itérations effectuées pour atteindre le critère d'arrêt.
  • Représentez la ligne brisée joignant les points $(x_0,x_0)$, $(x_0,x_1)$, $(x_1, x_1)$, $(x_1, x_2)$, $\ldots$, $(x_n,x_n)$.
  • Que remarquez-vous ?

Question

Pour chacun des cas tests ci-dessus, faites une figure illustrant l'ordre de convergence de la suite vers le point fixe (attention, le point fixe obtenu n'est pas toujours le même).

  • Vous pourrez tracer par exemple le nuage de points $(e_n, e_{n+1})$ où $e_n$ est l'erreur $|x_n-\ell|$ à l'étape $n$ avec $\ell$ la valeur limite du point fixe.
  • Pour visualiser l'ordre, vous pourrez ajouter aux graphiques des droites de pente $1$, $2$ ou $3$.

Pouvez-vous expliquer les différents comportements ?

In [57]:
## METTRE LE CODE ICI
size = 4
fig = plt.figure(figsize=(size*len(list_f), 2*size))
fig.suptitle("Méthode du point fixe", fontsize=20)

for k, dico in enumerate(list_f):
    a, b = dico['a'], dico['b']
    x0 = dico['x0']
    x = np.linspace(a, b, 1000)
    f = dico['expr']
    xstar, niter, xL = point_fixe(f, x0)
    
    ax = fig.add_subplot(2, len(list_f), k+1)
    ax.plot(x, f(x), color=color_bleu)
    ax.plot(x, x, color=color_vert, linestyle='dotted')
    xx = np.repeat(xL[:, None], 2, axis=1).flatten()
    ax.plot(xx[:-1], xx[1:], color=color_rouge)
        
    ax_e = fig.add_subplot(2, len(list_f), k+1+len(list_f))
    ax_e.set_xscale('log', base=2)
    ax_e.set_yscale('log', base=2)
    ax.scatter(xstar, f(xstar), s=20, color=color_rouge)
    ax.set_title(f'$f_{k}$ (${niter}$ itérations)', fontsize=16)
    e = abs(xL - dico['l'])
    ax_e.scatter(e[:-1], e[1:], color=color_orange, alpha=0.5)
    for p in [1, 2, 3]:
        ax_e.plot(e[:-1], e[:-1]**p, alpha=.5, label=f'pente ${p}$', color=liste_color[3+p])
    ax_e.legend()
    ax_e.set_title("Erreur", fontsize=16)
    ax_e.set_xlabel(r"$|e_{n}|$")
    ax_e.set_ylabel(r"$|e_{n+1}|$")
    ax.set_title(f'$f_{k}$ ({niter} itérations)', fontsize=16)
    
fig.tight_layout()
In [ ]: