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 Scilab. Ces programmes sont très facilement portables en n'importe quel langage: C, VB ou autres. En passant, je vous indiquerai les fonctions existantes de Scilab et Maple 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!
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 ε l'erreur qu'on accepte de commettre, alors le nombre k d'itérations minima est donné par k > (log2(b - a)/ε) - 1. 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
// 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 n'appelle pas de commentaire particulier, il 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.0e-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 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, Scilab n'étant pas un logiciel de calcul formel, 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 Maple ou un autre logiciel de calcul formel (Octave par exemple).
Le principe de la méthode de Newton 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. 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 x0 "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...
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 tolérance epsilon égale à 1.0e-6. C'est largement plus convergent que la dichotomie! Vous pouvez essayer avec d'autres valeurs de x0, 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..
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... Qui est dans le vrai? Voyons ce que nous indique Maple!
C'est encore plus simple qu'avec Scilab! 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 (du moins sur mon PC!) la solution : -4,273749 . C'est précisément le résultat obtenu par la méthode de Newton en 4 itérations...
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 - tangenteX.com mars 2013 
Contact : PhysiqueX ou 
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.