Méthodes de dérivation numérique

Rappels d'analyse

Je pars du principe que vous connaissez la définition lycéenne de la dérivée, celle que l'on apprend en Première :
Soit f une fonction continue dans un intervalle donné. Soit un point a appartenant à I. On dit que la fonction f est dérivable en a si et seulement si la limite lim (1/h)*( f(a+h) -f(a)), quand h tend vers 0, existe.

Cette limite est notée f'(a) dans vos livres de lycée. On la note plutôt f⁽¹⁾(a) pour tenir compte de l'existence possible de dérivées d'ordre supérieur.
A cette définition de la dérivée, il me faut ajouter un élément supplémentaire non connu des lycéens, mais fondamental en calcul numérique : le développement limité d'une fonction ou encore développement en série de Taylor. Je ne démontrerai rien ici, je me contenterai d'énoncer l'existence !
Donc, pour une fonction f analytique, dans un domaine de convergence de rayon R autour d'un point a, pour tout x tel que a-R <= x <= a+R, je peux écrire:
f(x) = f(a) + (x - a)*f⁽¹⁾(a) + (x - a)²/2! * f⁽²⁾(a) + ... + (x-a)n/n! * f⁽ⁿ⁾(a) + ...
Je vous laisse vous rafraichir la mémoire dans votre bouquin d'analyse préféré. Comme vous le savez sans doute, le calcul de la dérivée d'une fonction en un point est un exercice fréquent en physique et en chimie : calcul de vitesse, d'extrema de courbe, etc. Dans la suite on va donc chercher le moyen de calculer la dérivée première f⁽¹⁾ en un point.

Une méthode simple par différences finies

Tout d'abord, remarquons que la définition mathématique d'une dérivée est absolument inopérante à la fois pour les physiciens et pour les informaticiens. Pour la simple raison que faire tendre h vers 0 est impossible, à la fois physiquement et encore plus informatiquement.
Vous savez sans doute que si vous essayez de produire un nombre trop petit dans un programme, votre ordinateur protestera sévèrement ! En effet, un nombre est codé avec un certain nombre de bit (16, 32, 64 ou 128), un nombre toujours limité. Aussi, vous ne pouvez pas choisir un nombre aussi petit que nécessaire !
Il va donc falloir trouver autre chose que la pure définition mathématique de la dérivée. Nous allons utiliser un outil très commun en analyse numérique : les différences finies. C'est pour ça que j'ai introduit le développement d'une fonction en série de Taylor. Pas de panique, si on admet l'existence de ce développement, alors le calcul est de niveau Première ou Terminale....
Je cherche donc la valeur de la dérivée d'une fonction f en un point a. Je vais supposer que f est continue et dérivable, et aussi qu'en a, sa dérivée à gauche est égale à sa dérivée à droite. Cette dernière condition n'est pas obligatoire, mais je n'ai pas envie de vous parler de différence à gauche et de différence à droite... Vous creuserez ce cas à titre d'exercice !
Donc, je vais développer ma fonction f(x) sur un intervalle [-h,h] centré en a, en appliquant la formule de Taylor:
f(a + h) = f(a) + h*f⁽¹⁾(a) + h²/2! * f⁽²⁾(a) + h³/3! * f⁽³⁾(a) + h⁴/4! * f⁽⁴⁾(a) + O(h⁵) (1)
f(a - h) = f(a) - h*f⁽¹⁾(a) + h²/2! * f⁽²⁾(a) - h^>3/3! * f⁽³⁾(a) + h⁴/4! * f⁽⁴⁾(a) + O(h⁵) (2)
O(h⁵) désigne tous les termes d'ordre 5 et supérieurs, que je vais négliger sans vergogne, car je fais une erreur de l'ordre h⁵/5! * f⁽⁵⁾(a), en choisissant bien sur h très petit devant 1, par exemple 0.001 ou 0.0001.
Si je retranche (2) à (1), j'obtiens:
f(a + h) - f(a - h) = 2*h*f⁽¹⁾(a) + h³/3! * f⁽³⁾(a)
et donc:
f⁽¹⁾(a) # (1/2h)*(f(a + h) - f(a - h))
si je néglige les termes d'ordre 3 et les suivants. Le signe # signifie "approximativement égal".
Notons que si j'additionne (1) et (2), j'obtiens de la même manière une approximation de la dérivée seconde de f en a:
f⁽²⁾(a) # (1/h²)*(f(a + h) - 2*f(a) + f(a - h))
Vous comprenez sans doute maintenant pourquoi on appelle cette méthode "différences finies"... Il s'agit en l'occurrence ici de différences centrées.
Voyons ce que cela donne sur un exemple dont on connait déjà le résultat. Par exemple, je veux calculer la dérivée de la fonction sin(x) en a=0. Et je vais essayer avec plusieurs valeurs de h, par exemple 0.01 et 0.001
Nous savons bien sur que la dérivée de sin(x) est cos(x), ce qui nous permet de calculer la valeur théorique de la dérivée de sin(x) en a=0, soit cos(0) = 1.
Le programme de calcul se décompose en plusieurs routines:

la routine de définition de la fonction:

 double Fonction(double x)

        { 

           return sin(x);

        }

la routine de calcul de la dérivée:

 double DeriveOrdre1P3(double a, double h, double(*fonction)(double
        x))

        { 

           double x; 

           x = (fonction(a+h) - fonction(a-h))/(2*h);

           return x; 

        }

le programme principal:

 double x,a,h; 

        // saisie des données 

        cout << "Valeur de a: "; cin >> a; 

        cout << "Valeur de h: "; cin >> h; 

        x = DeriveOrdre1P3(a,h,Fonction);

        cout << "Derivee de f en " << a << ": " << x
        << endl;

Voyons les résultats obtenus:

pour a = 0 et h = 0.01, j'obtiens 0.999983
pour a = 0 et h = 0.001, j'obtiens 1
pour a = 3.141592 et h = 0.001, j'obtiens -1

Plutôt corrects, ces résultats. Il se peut cependant que pour certaines fonctions un peu "raides", les résultats ne soient pas très précis en utilisant une différence centrée sur 3 points. On peut améliorer un peu les choses en calculant la dérivée sur un voisinage plus étendu, en utilisant une différence à 5 points. Dans l'algorithme ci-dessus, j'ai négligé les termes en O(3). En travaillant un peu le calcul à partir des séries de Taylor (1) et (2), on peut obtenir une équation en O(5), ce qui améliore la précision de deux ordres de grandeur.
Au final, j'obtiens:
f⁽¹⁾(a) # (1/12h)*(f(a - 2*h) - 8*f(a - h) + 8*f(a + h) - f(a + 2*h))
ce qui nous donne la routine:

 double DeriveOrdre1P5(double a, double h, double(*fonction)(double
        x))

        { 

          double x; 

          x = (fonction(a-2*h) - 8*fonction(a-h) + 8*fonction(a+h) -
        fonction(a+2*h))/(12*h);

          return x;

        }

Vous pourrez contrôler qu'on obtient les mêmes résultats pour le calcul précédent.
Notons cependant que ces méthodes souffrent de quelques défauts, en particulier des erreurs d'arrondis et de troncatures. Elles peuvent être très gênantes sur des calculs précis.

Une méthode plus précise de dérivation

Il existe des méthodes qui améliorent grandement la technique des différences finies, en affinant le calcul des ordres supérieurs. Elles utilisent des algorithmes particuliers comme ceux de Neville, qu'il n'est pas question de décrire ici. Les "Numerical Recipes" nous fournissent une implémentation d'un algorithme de Neville revu et adapté par Ridders, qui est très efficace.
Vous en trouverez le code dans le programme Derivee.cpp.

Un programme C pour tester les différentes méthodes de dérivation

Le programme Derivee.cpp, téléchargeable, permet, sur une fonction déclarée dans la routine Fonction, de comparer la précision des différentes méthodes.

La dérivation avec Python

L’outil

Le package scipy de Python propose une routine de calcul de dérivée d'ordre n en un point, construite sur l'algorithme des différences centrées. Il s'agit de la routine derivative() que vous trouverez dans scipy.misc.

Les principaux arguments sont : le nom de la fonction à dériver, le point auquel vous voulez la dériver, le pas de calcul et l'ordre de la dérivation. Sur ce dernier point, assurez-vous que votre fonction est bien de classe C^_n, sinon vous aurez droit à un beau message d'erreur. Attention au pas de calcul dx, qui doit être assez petit, le "assez" dépendant de la fonction. Faites plusieurs essais...

Vous trouverez la documentation de la fonction derivative() sur cette page .

Exemples

Dériver la fonction sin(x)

Commençons par un exemple très simple : imaginons que je veuille calculer la dérivée de la fonction f(x) = sin(x) sur l'intervalle de 0 à 2*PI. Cette fonction est continue et dérivable sur le domaine, donc pas de problème.

Comme nous l'avons vu ci-dessus, la fonction derivative() calcule la dérivée de f(x) en un point. Il faut donc découper le domaine en un grand nombre de points, calculer la dérivée en chaque point et construire la courbe dérivée à l'aide de ces calculs. Pas très pratique...
Heureusement, Python qui sait travailler avec des vecteurs nous permet de réaliser ce calcul très simplement. Voyons comment dans le script BODerivationSinus.py.

Pour commencer les habituels appels de librairie :

 from scipy import sin, pi, arange

        from scipy.misc import derivative

        from matplotlib.pyplot import *

La fonction derivative() est dans le package scipy.misc.

Je construit maintenant le vecteur x, à l'aide de la fonction arange(). Puis je calcule le vecteur y = f(x):

 x = arange(0,2*pi,0.01)

        y = f(x)

Et maintenant l'objet de la manoeuvre : je vais dériver la fonction f(x) pour tous les points du domaine et construire le vecteur yprime qui contiendra ces dérivées. J'appelle la fonction derivative(), en lui indiquant que je souhaite une dérivée première (le n = 1) et que le pas de dérivation est égal à 0.01 (le dx = 0.01). C'est bien sur le même que pour arange(). Cela donne :

yprime = derivative(f,x,n=1,dx=0.01)

Ensuite viennent les instructions de tracé de courbe.

Et voilà ce que cela donne :

Vous reconnaissez bien sur la courbe représentative de cos(x) en la courbe de la fonction dérivée !

Un autre exemple de l'utilisation de la fonction derivative()

Le calcul des dérivées intervient en optimisation lorsqu’il s’agit de rechercher le minimum d’une fonction. Il est alors associé au calcul des racines d’une équation que nous avons vu plus haut. Voyons un exemple simple et très classique. Une usine fabrique des boites de conserve de 1 l en métal et voudrait économiser le plus possible sur ses coûts de fabrication. On va supposer que les coûts de procédé ont été optimisés et qu’il ne reste plus qu’à travailler sur les coûts de matières premières. L’objectif est donc de réduire au mieux la surface de métal nécessaire pour englober un volume d’un litre, en supposant l’épaisseur du métal constante.

Appelons h la hauteur de la boite et r son rayon. La surface du cylindre que constitue la boite s’exprime donc par \(S = 2 \pi r h + 2\pi r^2 \) et son volume \( \pi r^2h \), si j’exprime les grandeurs en décimètres. Rappelons que 1 litre = 1dm^3 et donc que \( \pi r^2 h = 1\). Je me sers de cette dernière égalité pour exprimer S en fonction de r uniquement, ce qui me donne \(S = \dfrac{2}{r} + 2\pi r^2 \).

Il me reste à calculer la dérivée de cette fonction, et à chercher sa racine. Le script BODeriveeOpt.py va le faire pour moi !

Après avoir défini la fonction de calcul de la surface du cylindre et le domaine de variation du rayon (attention à éliminer le rayon nul!), je calcule la dérivée de la fonction Surface et je cherche le rayon ropt pour lequel elle s'annule.

Le programme trace les deux courbes S et S', ce qui nous permet de lire le rayon optimum, qui est le même que celui qu'affiche le programme, soit 0,54 dm.
Voici les deux courbes :