TangenteX.com

Poser le problème

La recherche des zéros d'une fonction réelle à variable réelle est une occupation fort prenante en analyse numérique et donc en physique numérique. Rappelons que l'on appelle "zéro" d'une fonction, ou encore "racine", la(les) solution(s) de l'équation f(x) = 0. En anglais, cela se nomme "root finding", terme que vous trouverez fréquemment dans la littérature.

Dans le cas d'une fonction linéaire, le problème est plutôt trivial : le zéro d'une fonction de type y = ax + b est assez simple à trouver! Il n'en est pas de même pour le reste des fonctions, les plus nombreuses! Prenons par exemple, la fonction carré: c'est encore simple (!) . Pour une fonction polynôme d'ordre 2, ça le fait aussi... ax2 +bx + c = 0, c'est encore faisable facilement... Mais pour les fonctions un peu plus compliquées, polynomiales d'ordre supérieur à 2, trigo ou exponentielle, le problème devient vite délicat !

Nous allons aborder deux méthodes très classiques dans cette page, avec des exemples de programmation en Python, Scilab et Maple. Ces programmes sont très facilement portables en n'importe quel langage: C, VB ou autres. En passant, je vous indiquerai les fonctions existantes dans ces langages pour traiter le problème.

Dans la suite, je ne traiterai que des fonctions réelles à une variable réelle, sachant que la transposition à des fonctions réelles à plusieurs variables réelles est relativement simple. Le cas des fonctions complexes sort du scope de notre site. Nous considérerons d'abord le cas d'une fonction simplement continue puis ensuite le cas d'une fonction continue et dérivable.

Pour illustrer le principe, j'ai choisi une fonction polynomiale classique f(x) = x3 + 4x2 + 5. Elle est continue et bien sur dérivable. Mais vous pourrez vous faire la main sur une autre fonction!

La méthode par dichotomie

C'est une méthode assez intuitive, qui a le bon goût d'être peu exigeante. Il suffit que la fonction f dont on cherche le ou les zéros soit continue! Point n'est besoin qu'elle soit dérivable comme pour la méthode de Newton que nous allons voir plus loin.

Le principe est simple: on choisit un intervalle [a,b] dans lequel la fonction f s'annule. Pour simplifier, on supposera que la fonction ne s'annule qu'une fois sur cet intervalle Cet intervalle peut être aussi grand que l'on veut, mais bon, on essayera de le limiter.. Pour s'assurer que la fonction s'annule dans l'intervalle, il suffit de vérifier le signe du produit f(a).f(b) : s'il est négatif, c'est bon !

On découpe notre intervalle [a,b] en deux sous-intervalles égaux, dont le nom de la méthode (dichotomie ou bisection). On vérifie l'intervalle dans lequel la fonction s'annule, toujours par le même moyen, en calculant le signe de f(an)f(bn), si an et bn sont les bornes du sous-intervalle.

Et on recommence: on découpe ce sous-intervalle en deux, et on teste le signe du produit dans les sous-sous-intervalles. Et on recommence encore...

Mais combien de fois ? Une petite étude de convergence de la suite que constitue la dichotomie montre que si on appelle ? \( \epsilon \) l'erreur qu'on accepte de commettre, alors le nombre k d'itérations minima est donné par \( k \gt \left(\dfrac{log_2(b - a)}{\epsilon} - 1 \right) \). En fait, j'utiliserai un critère d'arrêt basé sur la largeur de l'intervalle: la recherche s'arrêtera lorsque l'intervalle sera plus petit qu'une valeur fixée (la tolérance..). Cette valeur déterminera la précision du calcul mais aussi sa durée...

Voici un petit programme Scilab, téléchargeable ici, qui utilise cette méthode:

//*******************************************************************
// recherche des zéros d'une fonction par la méthode de la dichotomie
// Dominique Lefebvre octobre 2012
// www.tangenteX.com
//*******************************************************************
funcprot(0); // pour éviter un message de warning (redef des fonctions)
  //Définition de la fonction à traiter
  function [y] = f(x)
  y = x^3 + 4*x^2 + 5;
endfunction
// Fonction de tracé de la fonction
function plot_fonction(x1,x2,fn,pas)
  x = x1:pas:x2;
  y = fn(x);
  plot2d(x,y,style = 2,axesflag = 5);
endfunction
// Méthode de la dichotomie - La fonction retourne la valeur approchée
// à epsilon près du zéro dans l'intervalle [a,b] et le nombre d'itérations
function [x0,nbiter] = dichotomie(a,b,fonction,eps)
  n = 0;
  Ya = f(a);
  while (abs(b - a) > 2*eps)
    n = n + 1;
    // split de l'intervalle [a,b]
    X = 0.5*( a + b);
    Y = f(X);
    // détermination du sous-intervalle contenant le zéro
    if (Ya*Y) < 0 then
      b = X
    else
      Ya = Y;
      a = X;
    end;
   end;
   // retour des résultats
   x0 = X;
   nbiter = n;
endfunction
//*************************************************************************
// Programme principal
//*************************************************************************
// Définition des paramètres
x1 = -5.0;
x2 = 0.0;
pas = 0.1;
epsilon = 1.0e-6;
// tracé de la fonction
plot_fonction(x1,x2,f,pas);
// résolution de f(x) = 0 par dichotomie
x0 = 0;
nbiter = 0;
[x0,nbiter] = dichotomie(x1,x2,f,epsilon);
printf("Zero : %f nb iterations: %d\n",x0,nbiter);

Ce code reprend l'algorithme que j'ai décris plus haut. Notez quand même que je prends la précaution de tracer le graphe de la fonction, histoire de vérifier la valeur approximative de la racine. On peut peut être aussi remarquer l'usage du nom d'une fonction en paramètre d'une autre fonction. C'est une propriété intéressante que ne possèdent pas tous les langages de programmation.

Si vous lancez l'exécution de ce programme dans Scilab, vous obtiendrez la valeur suivante pour la racine : -4,273750 en 22 itérations, avec une tolérance epsilon égale à 1.0-6. 22 itérations, c'est beaucoup, la convergence de la méthode est assez lente!

L'usage de cette méthode présente aussi un autre inconvénient : il faut connaitre un intervalle approximatif qui renferme la racine. On le choisit généralement en traçant le graphe de la fonction... Avantage de la méthode : elle est simple et applicable à une fonction non dérivable, ce qui peut être intéressant dans certains cas.

Si vous voulez calculer la racine d'une autre équation, rien de plus simple: changez la définition dans la fonction f(x) et adaptez selon vos besoins les valeurs de x1 et x2 du programme principal, qui définissent l'intervalle qui contient la racine.

Une remarque avant de passer à Newton : vous aurez observé, car vous êtes attentifs, que je n'ai pas contrôlé le nombre maximum d'itérations, ce qui n'est pas très prudent. En effet, si la dichotomie diverge, mon programme va boucler indéfiniment ! Mauvais plan pour une utilisation pratique! Il faudrait donc rajouter un test sur un nombre max d'itérations, par exemple 100, et sortir si la convergence en epsilon n'est pas obtenue au bout de ce nombre d'itérations. A vos claviers...

La méthode de Newton-Raphson

Le principe

La méthode de Newton-Raphson nécessite la dérivabilité de la fonction, au moins au premier ordre. C'est généralement le cas pour les fonctions courantes, mais c'est quand même une contrainte qu'il faut vérifier. De plus, il faudra donc calculer la dérivée. Si ce n'est pas faisable analytiquement, il faudra coupler une routine de calcul de dérivée au programme de calcul des racines. Ou alors utiliser le module sympy de Python ou Maple ou un autre logiciel de calcul formel (Octave par exemple).

Le principe de la méthode de Newton-Raphson est relativement simple: on suppose qu'en un point xk, le zéro de la fonction f(x) est "proche" du zéro de la tangente de f(x). Or je sais que l'équation de la tangente en xk est y = f(xk) + f'(xk)(x - xk). Si je pose x = xk+1 et que j'identifie à 0, j'obtiens f(xk) + f'(xk)(xk+1 - xk) = 0 d'où la suite récurrente xk+1 = xk - f(xk)/ f'(xk).

En initiant la suite avec une valeur xi "proche" de la racine, on converge rapidement vers le résultat recherché. Petit inconvénient, il faut trouver une valeur approximativement proche de la racine. Comme pour la dichotomie, on regarde le graphe de la fonction...

Vous constaterez que la méthode de Newton converge bien plus vite que la dichotomie, à condition que x0 soit bien choisi (pas trop distant!) et que le zéro recherché soit simple (i.e. f'(zéro) <> 0). Je ne rentrerai pas dans les détails, si cela vous intéresse vous pourrez consulter votre cours de calcul numérique...

Implémentation avec Scilab

Je vous propose pour commencer une implémentation simple de la méthode de Newton-Raphson avec Scilab.

function [x0,nbiter] = Newton(xi,f,df,eps)
  // Initialisation
  n = 0;
  X = xi;
  // Boucle de calcul
  while (abs(f(X)) > eps)
    n = n + 1;
    X = X - f(X)/df(X);  
  end;
  // Retour des résultats
  x0 = X;
  nbiter = n;   
endfunction

Je passe à la fonction Newton() le point xi d'initialisation de la recherche, la fonction f(x) et sa dérivée f'(x) et la précision souhaitée. la fonction me retourne la valeur approchée de la racine de la fonction dans l'intervalle et le nombre d'itérations de calcul.

Voyons maintenant le programme entier. Il faut d'abord définir la fonction f(x) et sa dérivée f'(x), ce que je fais ainsi, en prenant la fonction de mon exemple ci-dessus :

function [y] = f(x)
  y = x^3 + 4*x^2 + 5;
endfunction

function [y] = df(x)
  y = 3*x^2 + 8*x;
endfunction

Le code du programme principal est :

// Bornes de l'intervalle de recherche
x1 = -5.0;
x2 = 0.0;
// Pas de tracé
pas = 0.1;
// Précision souhaitée
epsilon = 1.0e-6;

// Tracé de la fonction
plot_fonction(x1,x2,f,pas);

// Résolution de f(x) = 0 par la méthode de Newton-Raphson
xi = -4;
nbiter = 0;
[x0,nbiter] = Newton(xi,f,df,epsilon);
printf("Zero : %f  nb iterations: %d\n",x0,nbiter);

Vous pourrez télécharge le programme Scilab de recherche des racines par la méthode de Newton, téléchargeable ici.

Si vous lancez l'exécution de ce programme dans Scilab, vous obtiendrez la valeur suivante pour la racine : -4,273749 en 4 itérations, avec une précision epsilon égale à 1.0-6. La convergence est bien plus rapide que celle de la dichotomie ! Vous pouvez essayer avec d'autres valeurs de xi, en obtenant des résultats semblables.
Vous noterez que les valeurs obtenues pour le zéro sont très similaires, à 10-6 près ! Même remarque à propos de la limitation du nombre d'itérations que pour la dichotomie...

Implémentation avec Python

La fonction Newton() proposée ci-dessus souffre d'un défaut qui peut être rédhibitoire. Il se niche dans un détail, plus précisément dans l'instruction X = X - f(X)/df(X). Imaginez que lors du calcul, la valeur de la dérivée df(x) devienne nulle, notre ordinateur n'apprécierait pas du tout : plantage assuré !

Elle présente un autre défaut, qui n'aboutirait pas à un plantage mais à un bouclage infini, ce qui n'est pas mieux. Nous avons posé l'hypothèse, non démontrée, que notre suite récurrente était convergente. Et si elle ne l'était pas ? Notre fonction bouclerait infiniment...

Je vous propose donc de remédier à ces deux inconvénients désagréables dans la fonction NewtonRaphson() en Python que voici :

def NewtonRaphson(f,df,x,eps,NbMaxIter):
    X = f(x)
    n = 0
    while abs(X) > eps and n <= NbMaxIter:
        Xpoint = df(x)
        try:
            x = x - X/Xpoint
            n += 1
            X = f(x)
        except ZeroDivisionError:
            x = 0.0
            n =0
    return x,n

Le nombre d'itérations est limité par la borne NbMaxIter définie par le programmeur. En cas de division par zéro, une exception sera levée et le script s'arrêtera normalement sans plantage. A part ces deux modifications, le code est identique à la fonction Newton() de Scilab.

Le script Python est très similaire au programme Scilab :

# paramètres de calcul
xi = -4.0
epsilon = 1.0e-6
NbIter = 100

x0,n = NewtonRaphson(f,df,xi,epsilon,NbIter)

print'racine : ',x0, ' nb iterations : ',n

Remarque : dans la fonction NewtonRaphson(), je ne fais aucun test sur les variables d'entrée. Les règles de l'art voudraient que je teste que les bornes de l'intervalle de recherche ne soient pas nulles et de signe opposé.

En lançant le script, nous obtenons le résultat suivant, heureusement identique à celui obtenu avec Scilab :

racine :  -4.27374868631  nb iterations :  4

Le script Python est téléchargeable ici.

Les commandes Python, Scilab et Maple

Root finding avec Python

Python dispose dans le package scipy.optimize d'une fonction fsolve(), qui permet la recherche des racines d'une fonction. Sa documentation complète est disponible ici.

En voici une implémantation simple : 

# importation des librairies
from scipy.optimize import fsolve

# Fonction dont on cherche les racines
def f(x):
    return x**3 + 4*x**2 + 5
   
# programme principal
xi = -4.0
x0 = fsolve(f,xi)
print'racine : ',x0

Le résultat obtenu est sans surprise :

racine :  [-4.27374869]

Root finding avec Scilab

Scilab permet de rechercher assez simplement la (les) racines(s) d'une fonction. Voyons comment procéder dans le cas de notre exemple. Après avoir lancé Scilab, vous obtenez le prompt de commande:

Scilab restitue le résultat. En l'occurence, il indique que la racine vaut -4,276096. Le second nombre que vous lisez (0,0007554) est la valeur de f(x0), x0 étant la valeur de la racine estimée.
Vous constatez qu'il y a un écart avec les résultats que nous obtenons par la dichotomie et par Newton. Il est aussi différent de celui obtenu par fsolve() de Python. C'est assez curieux...

Root finding avec Maple

 Vous lancez Maple pour obtenir le prompt. Puis vous tapez la commande fsolve(x^3 + 4x^2 + 5,x) . C'est tout ! Maple vous retourne immédiatement la solution : -4,273749 . C'est précisément le résultat obtenu par la méthode de Newton en 4 itérations...

Pour conclure

Vous voilà doté de tous les éléments pour rechercher la racine d'une fonction pas trop tordue. Faites donc des essais avec des fonctions de plus en plus complexes pour voir. Vous pouvez aussi réfléchir à la manière de procéder pour traiter les fonctions à zéros multiples...

Contenu et design par Dominique Lefebvre - www.tangenteX.com mars 2013 Licence Creative Commons Contact :