Dérivée partielle en Python

Bonjour,

je cherche à calculer les dérivées partielles d'une fonction en Python. Je sais faire pour une fonction d'une variable grâce au module scipy :
from scipy import misc

def f(x):
    return(x**2)
    
print(misc.derivative(f, 1))

Maintenant, je voudrais savoir faire cela pour une fonction de plusieurs variables, par exemple :
def f(x):
    return(x[0]**2 + x[1]**3 + x[2]**4)

Disons que je veux les dérivées partielles de $f$ en $(3, 1, 4)$.
Si je dérive par rapport à $x[0]$, il faudrait passer la fonction $y \mapsto f(y, x[1], x[2])$ en argument de derivative :
x = [3, 1, 4]
misc.derivative(lambda y: f([y, x[1], x[2]]), x[0])

Ça a l'air de faire ce que je veux. Donc je fais pareil avec les autres variables :
misc.derivative(lambda y: f([x[0], y, x[2]]), x[1])
misc.derivative(lambda y: f([x[0], x[1], y]), x[2])     

Mais je trouve des valeurs aberrantes pour ces commandes. Mes fonctions lambda ont l'air d'être définies correctement, j'ai pu calculer des valeurs de f grâce à celles-ci. Mais la dérivation plante et je ne vois pas pourquoi ?

Ensuite, je voudrai automatiser cela et donc il faudra que je sois capable de générer les listes $[y, x[1], x[2]]$, $[x[0], y, x[2]]$ et $[x[0], x[1], y]$. Comment faire cela sachant que $y$ est symbolique ? Si c'était une valeur, ce serait facile, il me suffirait de faire quelque chose comme :
y = 8
for k in range(len(x)):
    xtemp = x[::]
    xtemp[k] = y
    print(xtemp)
Mais en enlevant la première ligne, le script va planter puisque $y$ ne sera pas défini.

Réponses

  • Je me réponds à moi-même pour dire que j'ai écrit n'importe quoi : ma commande misc.derivative(lambda y: f([y, x[1], x[2]]), x[0]) calcule la dérivée de la fonction $y \mapsto f(y, x[1], x[2])$ en $x[0]$, ce qui n'a pas grand chose à voir en général avec la dérivée partielle par rapport à la première variable de la fonction $f$ en $(x[0], x[1], x[2])$. Donc c'est logique que mes résultats ne soient pas ceux attendus.

    Je repose donc ma question mais sans donner ce que j'ai fait car je ne vois pas du tout comment partir : comment calculer les dérivées partielles d'une fonction (disons de trois variables) en Python ?
  • Bonjour,
    Je ne connais pas cela mais une recherche dans moteur monopolistique amène sur cette page.
    EDP
    Bonne soirée !
    Diane Levels
  • Merci, je suis aussi tombé sur cette page et je pensais trouver la réponse à mes questions au paragraphe sur les EDP. Mais, si je comprends bien, ils discrétisent l'EDP pour la résoudre.
    J'ai envisagé cette méthode pour calculer mes dérivées partielles en utilisant $\partial_1 f(3,1,4) = \lim\limits_{t \to 0} \frac{f((3,1,4) + t(0,1,0)) - f(3,1,4)}{t}$, ce qui doit s'implémenter sous la forme :
    def f(x):
        return(x[0]**2 + x[1]**3 + x[2]**4)
        
    t = 10**(-5)
    x = [3, 1, 4]
    for k in range(len(x)):
        L = [0 for i in range(3)]
        L[k] = t
        print(1/t*(f([x[i] + L[i] for i in range(3)]) - f(x)))
    
    (Je précise que je n'ai pas la moindre idée de si c'est une bonne idée de prendre $t = 10^{-5}$ pour $t \to 0$ ; j'imagine que dans mon cas, ce n'est pas un gros problème vu que ma fonction est polynomiale, mais que ce n'est pas forcément terrible si elle est plus compliquée). Ça me donne un résultat correct, mais comme il y a la commande derivative dans le module scipy, j'imagine qu'on doit pouvoir le faire formellement ?
  • Stackoverflow propose différentes solutions. La mieux notée :
    def partial_derivative(func, var=0, point=[]):
        args = point[:]
        def wraps(x):
            args[var] = x
            return func(*args)
        return derivative(wraps, point[var], dx = 1e-6)
    
    Démo, avec la fonction foo définie par :
    def foo(x, y):
      return(x**2 + y**3)
    
    >>> partial_derivative(foo, 0, [3,1])
    6.0000000008386678
    >>> partial_derivative(foo, 1, [3,1])
    2.9999999995311555
    
    La valeur de dx proposée ne me semble pas satisfaisante : si x est de l'ordre de $ 10^{10} $, $ x+dx $ et $ x-dx $ seront arrondis à $ x $ et la dérivée estimée sera nulle. Je modifierais le code proposé en adaptant intelligemment le pas dx :
    from math import frexp
    from scipy.misc import derivative
    def partial_derivative(func, var=0, point=[]):
        args = point[:]
        def wraps(x):
            args[var] = x
            return func(*args)
        n=max(frexp(point[var])[1],0)
        return derivative(wraps, point[var], dx = 2**(n-18))
    
    La fonction frexp est la fonction qui à $ x $ associe fraction et exposant : frexp(x) renvoie le couple (a,n) tel que $ 0.5\leq |a|<1 $ et n est l'entier naturel tel que $ x=a2^{n} $. Cas particulier : frexp(0) renvoie le couple (0,0).

    Je vous laisse en exercice le soin de justifier le choix du pas (pourquoi 18 ?). La démo précédente donne (avec la fonction modifiée) :
    >>> partial_derivative(foo, 0, [3,1])
    6.0
    >>> partial_derivative(foo, 1, [3,1])
    3.0
    
  • Ajoutons que, d'après la documentation, la fonction dérivative de scipy se contente d'estimer une valeur de la dérivée en calculant le taux d'accroissement centré :
    $$ \frac{f(x+dx)-f(x-dx)}{2dx} $$
    En faisant un développement limité d'ordre $ 3 $ de $ f(x\pm dx) $, on peut vérifier que l'erreur théorique est de l'ordre de $ dx^{2} $ lorsque $ dx\to 0 $.

    Pour estimer la pertinence de l'estimation, il faut également tenir compte des éventuelles erreurs d'arrondi dans le calcul de $ x\pm dx $ (nulles pour le choix que je propose) et dans le calcul de l'image par $ f $ de $ x\pm dx $ (qui dépend de l'algorithme de calcul de $ f $...) ce qui permet de justifier le choix de $ dx $ que je propose.
  • OK, donc ça n'a pas de sens de vouloir à tout prix un calcul formel alors que Python fera à peu près le même calcul que je propose (à ceci près que je ne me suis posé aucune question sur le dx pour le moment, mais je vais regarder l'amélioration que tu proposes ).

    En fait, je suis tombé sur la solution proposée sur Stackoverflow, mais je ne comprends pas ce que fais la fonction wraps (enfin, je comprends qu'elle est là pour fournir les arguments à la fonction derivative, mais je ne comprends pas comment). Et comme je veux définir ma fonction $f$ sous la forme :
    def f(x):
        return(x[0]**2 + x[1]**3 + x[2]**4)
    
    et pas :
    def f(x, y, z):
        return(x**2 + y**3 + z**4)
    
    pour éviter de devoir connaître à l'avance le nombre d'arguments de $f$, je ne sais pas l'adapter.

    Si je détaille la fonction partial_derivative de Stackoverflow, elle prend en paramètres :
    - la fonction func dont on veut calculer les dérivées partielles
    - le numéro var de la variable par rapport à laquelle on veut dériver (par défaut choisi à 0)
    - le point sur lequel on veut faire ça (par défaut choisi à vide, ce qui me paraît idiot)
    Puis, on copie dans une liste args la liste point.
    On définit wraps, qui prend en paramètre x, on modifie la liste args en mettant x à la place var (on se retrouve donc avec une liste $[x, 1, 4]$ par exemple, mais je ne comprends pas pourquoi Python ne plante pas puisque $x$ n'est pas défini). Enfin, on retourne la fonction évaluée en *args, mais je ne comprends pas ce *args.

    En relisant mon texte, je crois comprendre pourquoi on met point par défaut à vide : comme le deuxième argument de la fonction est fixé par défaut à $0$, on est obligé de donner une valeur par défaut à tous les arguments suivants ?
  • Bon, c'est vraiment utile d'écrire ce qui gêne : sur la fonction wraps, ce qui m’empêche de comprendre est le *args. Je suis évidemment allé voir ce que ça fait et je lis que ça permet d'empaqueter des variables. Comme je veux passer une liste à ma fonction de test foo, je n'ai pas besoin d'empaqueter. Du coup, le script suivant fonctionne !
    def partial_derivative(func, var=0, point=[]):
        args = point[:]
        def wraps(x):
            args[var] = x
            return func(args)
        return misc.derivative(wraps, point[var], dx = 1e-6)
        
    def foo(x):
      return(x[0]**2 + x[1]**3)
    
Connectez-vous ou Inscrivez-vous pour répondre.