On rencontre assez souvent des équations implicites à résoudre, c'est à dire des équations qui contiennent l'inconnu à rechercher à plusieurs endroits et qui ne peuvent donc être calculés comme y=f(x), mais comme f(x,y)=0: il s'agit alors de trouver y pour que la fonction s'annule.

Voici le code proposé (il est largement auto-documenté):

#!/usr/bin/python
# -*- coding: utf-8 -*-
 
from __future__ import division
 
def fnimplicite(fn,x1,dx,eps=1.0e-15):
    """fnimplicite(fn,x1,dx,eps=1.0e-15): résolution d'une équation implicite par recherche du zéro d'une fonction"""
    # calcul de la 1ère valeur de la fonction pour le x initial
    f1=fn(x1)
    i=0  # initialisation d'un compteur de boucle
    while True:
        # on vérifie qu'on n'est pas en train de boucler pour rien
        i+=1
        if i>1000:
            raise ValueError ("impossible de trouver le zero de cette fonction")
        # calcul d'une 2ème valeur pour le x précédent + l'incrément
        x2=x1+dx
        f2=fn(x2)
        if abs(f2)<eps:
            # on a atteint la précision demandée: renvoyer x2
            break
        if ((f1>0)and(f2>0))or((f1<0)and(f2<0)):
            #f1 et f2 ont le même signe, ils sont donc du même côté de l'abscisse
            if abs(f2)>abs(f1):
                # on s'éloigne du zéro: revenir en arrière avec un incrément diminué
                dx*=-0.5
        else:
            # on a dépassé le zéro: revenir en arrière avec un incrément diminué 
            dx*=-0.5
        # on fait en sorte que le nouveau point soit désormais le point de base
        x1=x2
        f1=f2
        # et on boucle pour un nouvel essai!
    return x2

• Le 1er argument de cette fonction est le nom de la fonction (sans ses parenthèses!) dont on cherche le zéro.
• Le second argument est la 1ère valeur du paramètre cherché.
• Le 3ème argument est l'incrément de départ du calcul.
• Et le 4ème et dernier argument correspond à la valeur (proche de zéro) de la fonction qui arrêtera la recherche (par défaut: 1.0e-15).

Cette fonction converge assez rapidement vers la valeur cherchée, grâce à la division par 2 à chaque boucle, de l'incrément de départ. Par exemple, avec l'incrément de départ = 0.1, on atteint 8.88178419700125e-017 au bout de 50 boucles seulement (c'est facile à calculer: 0.1/2**50).

Attention: si la fonction dont on cherche le zéro ne s'annule jamais, elle va boucler à l'infini! C'est pour cela que j'ai ajouté une condition d'arrêt: si on a bouclé 1000 fois et qu'on a toujours pas trouvé, on ne trouvera plus. Cette condition déclenche alors une exception à récupérer dans un try: .. except:. Par exemple, si la fonction “mafonction(x)” n'a pas de solution x qui l'annule, voilà comment on peut récupérer cette situation (il faut importer le module sys avant si on veut afficher l'erreur renvoyé par raise):

try:
    print fnimplicite(mafonction,1,0.1))
except:
    print "%s" % sys.exc_info()[1]

Ce qui affichera:

impossible de trouver le zero de cette fonction

On peut aussi avoir une fonction qui s'annule pour plusieurs valeurs du paramètre cherché. Il vous appartient alors de donner la bonne valeur initiale pour la recherche, ainsi que le bon incrément avec le bon signe.

Dans certaines fonctions complexes, on pourrait aussi avoir la fonction qui passe à l'infini entre la valeur de départ donnée et la valeur cherchée. La fonction déclenchera ainsi une exception lors du calcul de f2. D'ou l'intérêt d'utiliser toujours try: except: pour la récupération de la valeur. Il faudra pallier ce problème en fixant une valeur de départ plus proche de la valeur cherchée.

Prenons un premier exemple simple

Calculer x tel que a*x + b*x**2 + c/x**3 = 0 avec a=0.509637427336286, b=-0.365419816785096 et c=0.0936145349893529

On construit la fonction:

def mafonction(x):
    a=0.509637427336286
    b=-0.365419816785096
    c=0.0936145349893529
    return a*x+b*x**2+ c/x**3

Appliquons maintenant la fonction de recherche en passant comme paramètre le nom de cette fonction (le “%r” permet d'avoir toutes les décimales):

try:
    print "%r" % fnimplicite(mafonction,1,0.1)
except:
    print "%s" % sys.exc_info()[1]

Ce qui affichera comme résultat:

1.4522568370062716

Vérification

print mafonction(1.4522568370062716)  # affiche 6.34908792208e-016 qui est bien < 1.0e-15

Si on a déjà une fonction qui comporte plusieurs arguments, comme par exemple mafonction(a,b,c,x), il faudra passer par l'intermédiaire de lambda:

try:
    f=lambda x: mafonction(0.509637427336286, -0.365419816785096, 0.0936145349893529, x)
    print "%r" % fnimplicite(f,1,0.1)
except:
    print "%s" % sys.exc_info()[1]

Avec lambda, on peut d'ailleurs utiliser directement l'expression:

try:
    f=lambda x: 0.509637427336286*x-0.365419816785096*x**2+0.0936145349893529/x**3
    print "%r" % fnimplicite(f,1,0.1)
except:
    print "%s" % sys.exc_info()[1]

Et si b devient positif (b=+0.365419816785096), cette expression n'a plus de solution, et la fonction de recherche va boucler et donc va générer une exception:

try:
    f=lambda x: 0.509637427336286*x+0.365419816785096*x**2+0.0936145349893529/x**3
    print "%r" % fnimplicite(f,1,0.1)
except:
    print "%s" % sys.exc_info()[1]

Affichera:

impossible de trouver le zero de cette fonction

La mensualité d'un crédit bancaire (voir la page de ce site credit_bancaire) est donnée par la fonction

M=menscredit(C,IA,N)

avec M = la mensualité, C = le capital emprunté, IA = l'intérêt annuel en % (par exemple 5 pour 5%), et N = le nombre de mois de remboursement.

Par exemple, avec C = 1000, IA = 5 et N = 24, on calcule: M = 43.8713897340686 (à arrondir au centime!)

Mais si on veut calculer à quel intérêt correspond C=1000, M=43 et N=24, cela demande la résolution d'une équation implicite, parce que IA est à plusieurs endroits dans la formule utilisée.

try:
    C=1000
    N=24
    M=43
    f=lambda IA: M-menscredit(C,IA,N)
    print "%r" % fnimplicite(f,5,0.1)    
except:
    print "%s" % sys.exc_info()[1]

Ce qui affiche:

3.0424694829586318

Vérification:

print menscredit(C,3.0424694829586318,N)  # affiche 43.0

La fonction de recherche étudiée jusqu'ici ne présuppose pas de rechercher dans une direction donnée. Par exemple:

print fnimplicite(sin,pi/4,0.1) trouveras 0 même si on a démarré avec pi/4.

Dans certain cas, on veut pouvoir chercher dans une direction donnée. Voilà le code adapté à cet objectif:

#!/usr/bin/python
# -*- coding: utf-8 -*-
 
from __future__ import division
 
def fnimplicite(fn,x1,dx,eps=1.0e-15):
    """fnimplicite(fn,x1,dx,eps=1.0e-15): résolution d'une équation implicite par recherche du zéro d'une fonction"""
    # calcul de la 1ère valeur de la fonction pour le x initial
    f1=fn(x1)
    i=0  # initialisation d'un compteur de boucle
    while True:
        i+=1
        if i>1000:
            raise ValueError ("impossible de trouver le zero de cette fonction")
        # calcul d'une 2ème valeur pour le x précédent + l'incrément
        x2=x1+dx
        f2=fn(x2)
        if abs(f2)<eps:
            # on a atteint la précision demandée: renvoyer x2
            break
        if ((f1>0)and(f2>0))or((f1<0)and(f2<0)):
            x1=x2
            f1=f2
        else:
            #f1 et f2 sont de part et d'autre de l'abscisse
            dx*=0.5
    return x2

Avec ce code:

print fnimplicite(sin,pi/4,0.1) trouva bien que le prochain passage de sin(x) à zéro à partir de pi/4 se produira pour x=pi.

L'incrément de départ doit donner la direction de recherche. Ainsi:

print fnimplicite(sin,pi/4,-0.1) trouva bien que le prochain passage de sin(x) à zéro à partir de pi/4 mais dans le sens des aiguilles d'une montre se produira bien pour x=0.

Et, bien sûr, si dans la direction cherchée la fonction ne passe jamais à zéro, on se retrouvera avec l'exception correspondant à la condition de bouclage.

Amusez-vous bien!