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
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!